CN114004169A - 一维饱和-非饱和水分运动参数计算方法及装置 - Google Patents

Abstract

本发明提供基于水均衡的一维饱和‑非饱和水分运动参数计算方法及装置，能够准确地模拟实际非均值土的完整的饱和‑非饱和带水分运动，计算速度快，计算效率高。方法包括：步骤1.基础资料收集；步骤2.垂向计算厚度确定和静态网格划分；步骤3.动态网格划分；步骤4.运行非饱和带模块；步骤4.1计算入渗水分配；步骤4.2计算重力作用下的土壤水分运动；步骤4.3计算源汇项作用下的土壤水分运动；步骤4.4计算扩散作用下的土壤水分运动；步骤5.运行饱和带模块；步骤6.分配饱和带与非饱和带之间的净通量；步骤7.采用上述步骤3～6对各时间步下的饱和‑非饱和水分运动参数均进行计算，得到计算区各时间步下各网格的运动参数值。

Description

一维饱和-非饱和水分运动参数计算方法及装置

技术领域

本发明属于地下水分及溶质计算领域，具体涉及基于水均衡的一维饱和-非饱和水分运动参数计算方法及装置。

技术背景

在干旱半干旱地区，农业用水占经济社会用水总量的70％以上。在地表水资源稀缺的地区，地下水对于国家水安全、粮食生产、生态文明建设具有重要意义。对于存在大量引水灌溉的农业灌溉，地下水埋深较浅，地下水和土壤水联系密切。准确模拟土壤含水率和地下水埋深及饱和带与非饱和带之间的交换通量是制定合理的灌溉策略和合理利用地下水的重要依据。制定合理的灌溉措施与合理利用地下水对于生态文明建设，国家水安全和粮食安全都具有重要意义。

目前常用于模拟饱和-非饱和水分运动的是基于水头型Richards方程的模拟方法。该技术至少存在以下问题：由于水头型Richards方程为非线性方程，求解过程复杂且耗时，计算过程中不能保证水均衡，因此在实际应用时存在困难。而基于含水量型Richards方程计算，则在土质变化时含水率梯度不能完成反映土壤水势梯度的变化，不能刻画非均质土的水分运动。同时由于含水率不能反映饱和带的水势梯度，因此该类技术方案不能用于模拟饱和带。

发明内容

本发明是为了解决上述问题而进行的，目的在于提供一种基于水均衡的一维饱和-非饱和水分运动参数计算方法及装置，能够准确地模拟实际非均值土的完整的饱和-非饱和带水分运动，计算得到土壤含水率和地下水埋深，并能得到饱和带与非饱和带之间的交换通量，计算速度快，计算效率高。

本发明为了实现上述目的，采用了以下方案：

<方法>

本发明提供了一维饱和-非饱和水分运动参数计算方法，其特征在于，包括以下步骤：

步骤1.基础资料收集：

收集研究区的气象资料、灌溉资料、土质资料、土壤含水率资料、地下水埋深资料；

步骤2.垂向计算厚度确定和静态网格划分：

根据计算时段内的最大埋深确定垂向计算深度(推荐比最大埋深更深的范围作为垂向计算厚度)，确定垂向网格划分厚度，根据这些基础信息进行静态垂向网格划分；

步骤3.动态网格划分：

根据地下水埋深位置在静态网格基础上增加一层；设当前正处于第t个时间步，上一个时间步计算得到的地下水埋深为D^t-1，找到地下水埋深所在的第j层，将第j层从地下水埋深处分开，上一层为最后一层非饱和带计算层，下一层为第一层饱和带计算层，这2层的底部距地表的距离为：

z^j，t＝D^t-1 (3-1)

z^j+1，t＝z^j，t-1 (3-2)

式中，z^j,t和z^j,t-1为t和t-1时间步第j层的底部距地表的距离；z^j+1,t为t时间步第j+1层距地表的距离；D^t-1为t-1时间步的地下水埋深；

步骤4.运行非饱和带模块：

非饱和带模块基于含水量型Richards方程，在一个时间步内非饱和带模块包括4个分步骤，更新非饱和带的土壤含水率，并得到饱和带与非饱和带之间的各项交换通量，具体计算方法见以下分步骤；

步骤4.1计算入渗水分配；

入渗水的分配采用“tipping-bucket”方法，入渗水自上而下按饱和含水率进行填充，入渗水分配后，第j层的填充水量为：

式中，q^j为单位面积的第j层土壤所分配的水量；M^j为第j层土壤的厚度；

是第j层土壤的饱和含水量；θ^j是第j层土壤的含水率；I是单位面积的入渗水的总量；

是第j层以上土壤入渗过程中已经分配过的入渗水量；本发明针对一维问题，如果入渗水将整个计算区填至饱和，则多余的水分为地表径流，忽略地表径流；

步骤4.2计算重力作用下的土壤水分运动；

当土壤含水率超过田间持水率时，土壤水分在重力作用下向下层排水，其计算公式为，

式中，K(θ)为非饱和水力传导度，采用以下公式计算：

式中，θ为土壤含水率，K_s为饱和水力传导度，θ_f为田间持水率；田间持水率θ_f受地下水埋深的影响，采用以下公式计算：

θ′_f＝ωθ_f (4-4)

式中，θ′_f为修正后田间持水率；w为田间持水率修正系数；h_d为田间持水率受地下水埋深影响的高度；D为地下水埋深；

求解重力作用下的土壤水分运动方程后，可得到通过重力作用补给饱和带的水通量，记为DQ；

步骤4.3计算源汇项作用下的土壤水分运动；

在源汇项作用下，土壤水分运动方程如下：

式中，W是源汇项，主要包括土壤蒸发和作物蒸腾；当实际蒸发蒸腾的消耗深度比地下水埋深更深时，将会消耗地下水，产生向上的水通量，记为UQ；

步骤4.4计算扩散作用下的土壤水分运动；

采用含水率型方程描述的扩散作用下(基质势梯度作用)的土壤水分运动，其计算公式：

式中，D(θ)为扩散度，p_l(l＝1,…,L)是土壤水分特征曲线的参数，L为参数个数；

为描述土质垂向变异性的修正项；在扩散作用下，饱和带的水分将会补给至非饱和带，通过的水量记为DF；

步骤5.运行饱和带模块；

一个时间步内，经过入渗水的分配及在重力势、源汇项和扩散项作用下的土壤水分运动，计算得到在地下水面的交换通量，计算公式如下：

FQ＝DQ-UQ-DF (5-1)

式中，FQ为饱和带与非饱和之间的净通量；DQ为重力势作用下向饱和带补给的通量；UQ为源汇项作用下饱和带向非饱和带补给的通量；DF为扩散作用下饱和带向非饱和带补给的通量；

根据饱和非饱和界面通量，采用地下水埋深波动法更新地下水埋深：

式中，D^t和D^t-1为第t和t-1时间步的地下水埋深；μ为给水度；

步骤6.分配饱和带与非饱和带之间的净通量；

当地下水埋深更新后，将饱和带与非饱和带之间的净通量分配到地下水埋深变动带，计算公式：

式中，θ^t为t时间步更新后的土壤含水率；θ_k ^t为t时间步内经过扩散作用计算后的土壤含水率；θ_s为饱和含水率；

步骤7.采用上述步骤3～6对各时间步下的饱和-非饱和水分运动参数均进行计算，得到计算区各时间步下各网格的运动参数值。

优选地，本发明提供的一维饱和-非饱和水分运动参数计算方法，还可以具有以下特征：运动参数包括：土壤水分、地下水埋深、非饱和带向饱和带的补给通量、饱和带向非饱和带的补给通量以及饱和带与非饱和之间的净通量。

优选地，本发明提供的一维饱和-非饱和水分运动参数计算方法，还可以具有以下特征：在步骤2中，采用最大埋深+1m的深度作为垂向计算厚度。

优选地，本发明提供的一维饱和-非饱和水分运动参数计算方法，还可以具有以下特征：在步骤2中，垂向网格划分厚度设置为0.05～0.20m之间。

优选地，本发明提供的一维饱和-非饱和水分运动参数计算方法，还可以具有以下特征：在步骤3中根据地下水埋深位置将地下水位所在分层分为2层，为了后续计算，在一个时间步结束后，需要将网格调整回静态网格，即合并计算过程中地下水埋深以上的第一层(j层)和地下水埋深以下的第一层(j+1层)为新的第j层，合并层的含水率为：

式中：

为合并后新的第j层的含水率；θ^j和θ^j+1为第j层和j+1层的含水率；M^j和M^{j +1}为第j层和第j+1层的厚度。

<装置>

进一步，本发明还提供了一种自动实现上述<方法>的区域饱和-非饱和水分及溶质运移计算装置，其特征在于，包括：

资料获取部，获取研究区的气象资料、灌溉资料、土质资料、土壤含水率资料、地下水埋深资料；

静态网格划分部，根据计算时段内的最大埋深确定垂向计算深度，确定垂向网格划分厚度，并根据这些基础信息进行静态垂向网格划分；

动态网格划分部，根据地下水埋深位置在静态网格基础上增加一层；设当前正处于第t个时间步，上一个时间步计算得到的地下水埋深为D^t-1，找到地下水埋深所在的第j层，将第j层从地下水埋深处分开，上一层为最后一层非饱和带计算层，下一层为第一层饱和带计算层，这2层的底部距地表的距离为：

z^j，t＝D^t-1 (3-1)

z^j+1，t＝z^j，t-1 (3-2)

式中，z^j,t和z^j,t-1为t和t-1时间步第j层的底部距地表的距离；z^j+1,t为t时间步第j+1层距地表的距离；D^t-1为t-1时间步的地下水埋深；

非饱和带模块运行部，在一个时间步内执行步骤4.1～4.4四个分步骤，更新非饱和带的土壤含水率，并得到饱和带与非饱和带之间的各项交换通量；具体计算方法见以下分步骤；

步骤4.1计算入渗水分配；

入渗水自上而下按饱和含水率进行填充，入渗水分配后，第j层的填充水量为：

式中，q^j为单位面积的第j层土壤所分配的水量；M^j为第j层土壤的厚度；

是第j层土壤的饱和含水量；θ^j是第j层土壤的含水率；I是单位面积的入渗水的总量；

是第j层以上土壤入渗过程中已经分配过的入渗水量；

步骤4.2计算重力作用下的土壤水分运动；

当土壤含水率超过田间持水率时，土壤水分在重力作用下向下层排水，其计算公式为，

式中，K(θ)为非饱和水力传导度，采用以下公式计算：

式中，θ为土壤含水率，K_s为饱和水力传导度，θ_f为田间持水率；田间持水率θ_f受地下水埋深的影响，采用以下公式计算：

θ′_f＝ωθ_f (4-4)

式中，θ′_f为修正后田间持水率；w为田间持水率修正系数；h_d为田间持水率受地下水埋深影响的高度；D为地下水埋深；

求解重力作用下的土壤水分运动方程后，可得到通过重力作用补给饱和带的水通量，记为DQ；

步骤4.3计算源汇项作用下的土壤水分运动；

在源汇项作用下，土壤水分运动方程如下：

式中，W是源汇项，包括土壤蒸发和作物蒸腾；当实际蒸发蒸腾的消耗深度比地下水埋深更深时，将会消耗地下水，产生向上的水通量，记为UQ；

步骤4.4计算扩散作用下的土壤水分运动；

采用含水率型方程描述的扩散作用下的土壤水分运动，其计算公式：

式中，D(θ)为扩散度，p_l(l＝1,…,L)是土壤水分特征曲线的参数，L为参数个数；

为描述土质垂向变异性的修正项；在扩散作用下，饱和带的水分将会补给至非饱和带，通过的水量记为DF；

饱和带模块运行部，在一个时间步内，经过入渗水的分配及在重力势、源汇项和扩散项作用下的土壤水分运动，计算得到在地下水面的交换通量，计算公式如下：

FQ＝DQ-UQ-DF (5-1)

式中，FQ为饱和带与非饱和之间的净通量；DQ为重力势作用下向饱和带补给的通量；UQ为源汇项作用下饱和带向非饱和带补给的通量；DF为扩散作用下饱和带向非饱和带补给的通量；

埋深更新部，根据饱和非饱和界面通量，采用地下水埋深波动法更新地下水埋深：

式中，D^t和D^t-1为第t和t-1时间步的地下水埋深；μ为给水度；

净通量分配部，当地下水埋深更新后，将饱和带与非饱和带之间的净通量分配到地下水埋深变动带，计算公式：

式中，θ^t为t时间步更新后的土壤含水率；θ_k ^t为t时间步内经过扩散作用计算后的土壤含水率；θ_s为饱和含水率；

参数值获取部，采用动态网格划分部、非饱和带模块运行部、饱和带模块运行部、净通量分配部对每一时间步下的饱和-非饱和水分运动参数均进行计算，得到计算区各时间步下各网格的运动参数值；以及

控制部，与资料获取部、静态网格划分部、动态网格划分部、非饱和带模块运行部、饱和带模块运行部、净通量分配部、参数值获取部通信相连，控制它们的运行。

优选地，本发明提供的区域饱和-非饱和水分及溶质运移计算装置，还可以包括：输入显示部，与控制部通信相连，用于让用户输入操作指令，并进行相应显示。

优选地，本发明提供的区域饱和-非饱和水分及溶质运移计算装置，还可以具有这样的特征：输入显示部还能够根据操作指令对计算得到的饱和-非饱和水分运动参数按照时间顺序以动态演化的方式显示在相应的研究区模型图上。

优选地，本发明提供的区域饱和-非饱和水分及溶质运移计算装置，还可以具有这样的特征：运动参数包括：土壤水分、地下水埋深、非饱和带向饱和带的补给通量、饱和带向非饱和带的补给通量以及饱和带与非饱和之间的净通量。

发明的作用与效果

本发明基于水均衡方法，并分别计算了入渗水分配、重力作用下的土壤水分运动、源汇项作用下的土壤水分运动、扩散作用下的土壤水分运动情况，在计算过程中充分考虑了田间持水率受地下水埋深的影响，采用地下水埋深对田间持水率进行修正处理，同时本发明还引进了土质垂向变异性的修正项，可以有效刻画非均质情况，更能适应于实际情况，进一步本发明中考虑到扩散作用下的土壤水分通量，不仅可以模拟饱和带向非饱和带的通量，而且可以模拟非饱和带向饱和带的通量，对于地下水浅埋区和深埋区都同样适用，前述这些处理方法都极大地提高了本发明的模拟计算精度，使得本发明能够计算获得符合真实情况、准确可靠的土壤含水率、地下水埋深、非饱和带向饱和带的补给通量(水分通量)、饱和带向非饱和带的补给通量以及饱和带与非饱和之间的净通量，根据这些数据可以科学合理地确定灌溉措施、实施灌溉、合理利用水资源，例如，确定相应时期的灌溉量，将模拟计算结果与不同灌溉方案进行比较，选择满足模拟计算结果且用水量最少的灌溉方案作为最优灌溉方案并执行。

附图说明

图1为本发明实施例中涉及的基于水均衡的一维饱和-非饱和水分运动参数计算方法的流程图；

图2为本发明实施例中涉及的永联试验区地理位置图；

图3为本发明实施例中涉及的永联试验区2007年气象数据和灌溉数据图；

图4为本发明实施例中涉及的永联试验区6个监测点的垂向土质分布图；

图5为本发明实施例计算得到的永联试验区不同时间的地下水埋深和剖面含水率与实测值的对比图；

图6为本发明实施例计算得到的各监测点的地下水埋深和含水率与实测值的对比图。

具体实施方式

以下结合附图对本发明涉及的一维饱和-非饱和水分运动参数计算方法及装置进行详细地说明。

<实施例>

本实施例中，以中国内蒙古河套灌区永联试验区的6处农田监测点的饱和-非饱和水分运动模拟为例进行说明。本实施例目的是利用研究区2007年的基础数据，模拟2007年6处监测点的饱和-非饱和带的水分运动过程。如图1所示，本实施例所提供的基于水均衡的一维饱和-非饱和水分运动模拟方法，包括如下步骤：

步骤1.基础资料收集：

收集研究区的气象资料、灌溉资料、土质资料、土壤含水率资料、地下水埋深资料。

河套灌区永联试验区的分布图如图2所示，气象资料(降雨量、参考作物蒸发蒸腾量、实际蒸发蒸腾量)和灌溉资料如图3所示，垂向土质分布如图4所示，相应的土壤水力参数如表1所示。土壤含水率资料和地下水埋深资料将用于与模拟值对比。

表1不同土质水力参数取值表

步骤2.垂向计算厚度确定和静态网格划分：

永联试验区地下水埋深在0～3m之间，根据计算时段内的最大埋深确定垂向计算深度，本实施例以4m为垂向计算厚度。垂向网格划分厚度取为0.05m，垂向共分为80层。计算时间段为5月1日至11月30日。之后采用步骤3至6计算各时间步下的饱和-非饱和水分运动，得到计算区各时间步下各网格的土壤水分和地下水埋深。

步骤3.动态网格划分：

根据地下水埋深位置在静态网格基础之上增加一层，由80层增加至81层。设当前正处于第t个时间步，上一个时间步计算得到的地下水埋深为D^t-1，找到地下水埋深所在的第j层，将第j层从地下水埋深处分开，上一层为最后一层非饱和带计算层，下一层为第一层饱和带计算层，这2层的底部距地表的距离为，

z^j，t＝D^t-1 (3-1)

z^j+1，t＝z^j，t-1 (3-2)

式中，z^j,t和z^j,t-1为t和t-1时间步第j层的底部距地表的距离；z^j+1,t为t时间步第j+1层距地表的距离；D^t-1为t-1时间步的地下水埋深。

步骤4.运行非饱和带模块：

非饱和带模块基于含水量型Richards方程，在一个时间步内的计算包括：入渗水分配计算、重力作用下的土壤水分运动、源汇项作用下的土壤水分运动和扩散作用下的土壤水分运动。通过运行非饱和带模块更新非饱和带的土壤含水率，并将地下水面的各项交换通量传递给饱和带模块，用于更新地下水埋深。

步骤4.1计算入渗水分配；

入渗水自上而下按饱和含水率进行填充，入渗水分配后，第j层的填充水量为：

式中，q^j为单位面积的第j层土壤所分配的水量；M^j为第j层土壤的厚度；

是第j层土壤的饱和含水量；θ^j是第j层土壤的含水率；I是单位面积的入渗水的总量；

是第j层以上土壤入渗过程中已经分配过的入渗水量；

步骤4.2计算重力作用下的土壤水分运动；

当土壤含水率超过田间持水率时，土壤水分在重力作用下向下层排水，其计算公式为，

式中，K(θ)为非饱和水力传导度，采用以下公式计算：

式中，θ为土壤含水率，K_s为饱和水力传导度，θ_f为田间持水率；田间持水率θ_f受地下水埋深的影响，采用以下公式计算：

θ′_f＝ωθ_f (4-4)

式中，θ′_f为修正后田间持水率；w为田间持水率修正系数；h_d为田间持水率受地下水埋深影响的高度；D为地下水埋深；

求解重力作用下的土壤水分运动方程后，可得到通过重力作用补给饱和带的水通量，记为DQ；

步骤4.3计算源汇项作用下的土壤水分运动；

在源汇项作用下，土壤水分运动方程如下：

式中，W是源汇项，包括土壤蒸发和作物蒸腾；当实际蒸发蒸腾的消耗深度比地下水埋深更深时，将会消耗地下水，产生向上的水通量，记为UQ；

步骤4.4计算扩散作用下的土壤水分运动；

采用含水率型方程描述的扩散作用下的土壤水分运动，其计算公式：

式中，D(θ)为扩散度，p_l(l＝1,…,L)是土壤水分特征曲线的参数，L为参数个数；

为描述土质垂向变异性的修正项；在扩散作用下，饱和带的水分将会补给至非饱和带，通过的水量记为DF；

步骤5.运行饱和带模块：

一个时间步内，经过入渗水的分配及在重力势、源汇项和扩散项作用下的土壤水分运动，计算得到了在地下水面的交换通量，计算公式如下，

FQ＝DQ-UQ-DF (5-1)

式中，FQ为饱和带与非饱和之间的净通量；DQ为重力势作用下向饱和带补给的通量；UQ为源汇项作用下饱和带向非饱和带补给的通量；DF为扩散项作用下饱和带向非饱和带补给的通量。

接下来根据饱和非饱和界面通量，采用地下水埋深波动法更新地下水埋深，

式中，D^t和D^t-1为第t和t-1时间步的地下水埋深；μ为给水度。

步骤6.分配饱和带与非饱和带之间的净通量：

当地下水埋深更新后，将饱和带与非饱和带之间的净通量分配到地下水埋深变动带，计算公式为，

式中，θ^t为t时间步更新后的含水率；θ_k ^t为t时间步内经过扩散作用计算后的含水率；θ_s为饱和含水率。

随后将动态网格调整回静态网格。经过步骤3～6，一个时间步的土壤含水率和地下水埋深已全部更新。

以上计算是针对一个监测点，6个监测点的计算方法全部相同。根据以上计算得到2007年永联试验区6个监测点的饱和-非饱和水分运动过程，并将计算结果与实测值对比：

图5展示了永联试验区不同时间的地下水埋深和剖面含水率与实测值的对比，图5中的值为6个监测点的平均值。以纳什效率系数(NSE)、均方根误差(RMSE)和平均绝对百分比误差(MAPE)作为评价指标。整体的地下水埋深模拟精度高，其中NSE为0.923，RMSE为0.180m。图5种展示了模拟期内15次剖面含水率模拟值和监测值的对比结果，在大部分模拟日期NSE均大于0.60，RMSE和MAPE小于0.040cm³cm^-3和10.5％。图6展示了各监测点的地下水埋深和含水率与实测值的对比。各点的土壤含水率的NSE几乎均大于0.50，RMSE和MAPE小于0.060cm³cm^-3和15.0％。各点的NSE均大于0.70，RMSE均小于0.357m。

根据以上模拟结果和模拟指标可知，本发明可以较好地模拟实际情况下的饱和-非饱和水分运动，计算得到的参数值与实测值非常接近，具有很高的模拟和计算精度。

另外，由于本发明的计算方法无需迭代，本实施例中每个监测点的计算时长在2s内，累计计算时长在12s内；同时本发明基于水均衡模型，计算过程中水均衡误差在10^-5内。

此外，将本发明用于实际灌溉过程可以：1)评价现有灌溉措施(灌溉水量和灌溉日期)是否合理。如果计算出来的埋深过深，含水率过小，则说明水量偏少。如果计算出来某些日子含水率过大，有些日值含水率过小，则说明灌溉日期设置不合适。如果计算出来的含水率整体偏大，埋深整体偏浅，则说明灌溉水量不足。例如，根据以上方法计算得到的土壤含水率和地下水埋深等参数值可以判断土体水分是否满足作物生长需要，从而确定灌溉措施：假如取土体计划湿润层为80cm，要求的0-80cm平均体积含水率不得小于0.20，采用本发明模拟计算后发现小于了0.20，则表明此时需要灌溉；若加入要求的地下水埋深不得小于3m，计算过程中有出现地下水埋深超过3m的，则在当日需施加灌溉。2)设置不同的灌溉序列(包括灌溉水量和灌溉日期)，通过本发明进行模拟，选择出最优的灌溉措施(即用水量最少)。

进一步，本实施例还提供能够自动实现上述方法的一维饱和-非饱和水分运动参数计算装置，该装置包括资料获取部、静态网格划分部、动态网格划分部、非饱和带模块运行部、饱和带模块运行部、净通量分配部、参数值获取部、输入显示部以及控制部。

资料获取部用于获取研究区的气象资料、灌溉资料、土质资料、土壤含水率资料、地下水埋深资料。

静态网格划分部根据计算时段内的最大埋深确定垂向计算深度，确定垂向网格划分厚度，并根据这些基础信息进行静态垂向网格划分；

动态网格划分部根据地下水埋深位置在静态网格基础上增加一层；设当前正处于第t个时间步，上一个时间步计算得到的地下水埋深为D^t-1，找到地下水埋深所在的第j层，将第j层从地下水埋深处分开，上一层为最后一层非饱和带计算层，下一层为第一层饱和带计算层，这2层的底部距地表的距离为：

z^j，t＝D^t-1 (3-1)

z^j+1，t＝z^j，t-1 (3-2)

式中，z^j,t和z^j,t-1为t和t-1时间步第j层的底部距地表的距离；z^j+1,t为t时间步第j+1层距地表的距离；D^t-1为t-1时间步的地下水埋深；

非饱和带模块运行部在一个时间步内执行步骤4.1～4.4四个分步骤，更新非饱和带的土壤含水率，并得到饱和带与非饱和带之间的各项交换通量；具体计算方法见以下分步骤；

步骤4.1计算入渗水分配；

入渗水自上而下按饱和含水率进行填充，入渗水分配后，第j层的填充水量为：

式中，q^j为单位面积的第j层土壤所分配的水量；M^j为第j层土壤的厚度；

是第j层土壤的饱和含水量；θ^j是第j层土壤的含水率；I是单位面积的入渗水的总量；

是第j层以上土壤入渗过程中已经分配过的入渗水量；

步骤4.2计算重力作用下的土壤水分运动；

当土壤含水率超过田间持水率时，土壤水分在重力作用下向下层排水，其计算公式为，

式中，K(θ)为非饱和水力传导度，采用以下公式计算：

式中，θ为土壤含水率，K_s为饱和水力传导度，θ_f为田间持水率；田间持水率θ_f受地下水埋深的影响，采用以下公式计算：

θ′_f＝ωθ_f (4-4)

式中，θ′_f为修正后田间持水率；w为田间持水率修正系数；h_d为田间持水率受地下水埋深影响的高度；D为地下水埋深；

求解重力作用下的土壤水分运动方程后，可得到通过重力作用补给饱和带的水通量，记为DQ；

步骤4.3计算源汇项作用下的土壤水分运动；

在源汇项作用下，土壤水分运动方程如下：

式中，W是源汇项，包括土壤蒸发和作物蒸腾；当实际蒸发蒸腾的消耗深度比地下水埋深更深时，将会消耗地下水，产生向上的水通量，记为UQ；

步骤4.4计算扩散作用下的土壤水分运动；

采用含水率型方程描述的扩散作用下的土壤水分运动，其计算公式：

式中，D(θ)为扩散度，p_l(l＝1,…,L)是土壤水分特征曲线的参数，L为参数个数；

为描述土质垂向变异性的修正项；在扩散作用下，饱和带的水分将会补给至非饱和带，通过的水量记为DF；

饱和带模块运行部，在一个时间步内，经过入渗水的分配及在重力势、源汇项和扩散项作用下的土壤水分运动，计算得到在地下水面的交换通量，计算公式如下：

FQ＝DQ-UQ-DF (5-1)

式中，FQ为饱和带与非饱和之间的净通量；DQ为重力势作用下向饱和带补给的通量；UQ为源汇项作用下饱和带向非饱和带补给的通量；DF为扩散作用下饱和带向非饱和带补给的通量；

埋深更新部根据饱和非饱和界面通量，采用地下水埋深波动法更新地下水埋深：

式中，D^t和D^t-1为第t和t-1时间步的地下水埋深；μ为给水度；

净通量分配部在地下水埋深更新后，将饱和带与非饱和带之间的净通量分配到地下水埋深变动带，计算公式：

式中，θ^t为t时间步更新后的土壤含水率；θ_k ^t为t时间步内经过扩散作用计算后的土壤含水率；θ_s为饱和含水率；

参数值获取部采用动态网格划分部、非饱和带模块运行部、饱和带模块运行部、净通量分配部对每一时间步下的饱和-非饱和水分运动参数均进行计算，得到计算区各时间步下各网格的运动参数值；运动参数包括：土壤水分、地下水埋深、非饱和带向饱和带的补给通量、饱和带向非饱和带的补给通量以及饱和带与非饱和之间的净通量。

输入显示部，用于让用户输入操作指令，并进行相应显示。例如，输入显示部显示提示信息让操作员输入相应的研究区资料信息，输入显示部还能够根据操作指令对计算得到的饱和-非饱和水分运动参数按照列表或者图示方式进行显示，并能够根据操作指令对计算得到的饱和-非饱和水分运动参数按照时间顺序以动态演化的方式显示在相应的研究区模型图上。

控制部与资料获取部、静态网格划分部、动态网格划分部、非饱和带模块运行部、饱和带模块运行部、净通量分配部、参数值获取部、输入显示部均通信相连，控制它们的运行。

以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的一维饱和-非饱和水分运动参数计算方法及装置并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。

Claims (10)

1.一维饱和-非饱和水分运动参数计算方法，其特征在于，包括以下步骤：

步骤1.基础资料收集：

收集研究区的气象资料、灌溉资料、土质资料、土壤含水率资料、地下水埋深资料；

步骤2.垂向计算厚度确定和静态网格划分：

根据计算时段内的最大埋深确定垂向计算深度，确定垂向网格划分厚度，根据这些基础信息进行静态垂向网格划分；

步骤3.动态网格划分：

根据地下水埋深位置在静态网格基础上增加一层；设当前正处于第t个时间步，上一个时间步计算得到的地下水埋深为D^t-1，找到地下水埋深所在的第j层，将第j层从地下水埋深处分开，上一层为最后一层非饱和带计算层，下一层为第一层饱和带计算层，这2层的底部距地表的距离为：

z^j，t＝D^t-1 (3-1)

z^j+1，t＝z^j，t-1 (3-2)

式中，z^j,t和z^j,t-1为t和t-1时间步第j层的底部距地表的距离；z^j+1,t为t时间步第j+1层距地表的距离；D^t-1为t-1时间步的地下水埋深；

步骤4.运行非饱和带模块：

在一个时间步内非饱和带模块包括4个分步骤，更新非饱和带的土壤含水率，并得到饱和带与非饱和带之间的各项交换通量，具体计算方法见以下分步骤；

步骤4.1计算入渗水分配；

入渗水自上而下按饱和含水率进行填充，入渗水分配后，第j层的填充水量为：

式中，q^j为单位面积的第j层土壤所分配的水量；M^j为第j层土壤的厚度；

是第j层土壤的饱和含水量；θ^j是第j层土壤的含水率；I是单位面积的入渗水的总量；

是第j层以上土壤入渗过程中已经分配过的入渗水量；

步骤4.2计算重力作用下的土壤水分运动；

当土壤含水率超过田间持水率时，土壤水分在重力作用下向下层排水，其计算公式为，

式中，K(θ)为非饱和水力传导度，采用以下公式计算：

式中，θ为土壤含水率，K_s为饱和水力传导度，θ_f为田间持水率；田间持水率θ_f受地下水埋深的影响，采用以下公式计算：

θ′_f＝ωθ_f (4-4)

式中，θ′_f为修正后田间持水率；w为田间持水率修正系数；h_d为田间持水率受地下水埋深影响的高度；D为地下水埋深；

求解重力作用下的土壤水分运动方程后，可得到通过重力作用补给饱和带的水通量，记为DQ；

步骤4.3计算源汇项作用下的土壤水分运动；

在源汇项作用下，土壤水分运动方程如下：

式中，W是源汇项，包括土壤蒸发和作物蒸腾；当实际蒸发蒸腾的消耗深度比地下水埋深更深时，将会消耗地下水，产生向上的水通量，记为UQ；

步骤4.4计算扩散作用下的土壤水分运动；

采用含水率型方程描述的扩散作用下的土壤水分运动，其计算公式：

式中，D(θ)为扩散度，p_l(l＝1,…,L)是土壤水分特征曲线的参数，L为参数个数；

为描述土质垂向变异性的修正项；在扩散作用下，饱和带的水分将会补给至非饱和带，通过的水量记为DF；

步骤5.运行饱和带模块；

一个时间步内，经过入渗水的分配及在重力势、源汇项和扩散项作用下的土壤水分运动，计算得到在地下水面的交换通量，计算公式如下：

FQ＝DQ-UQ-DF (5-1)

式中，FQ为饱和带与非饱和之间的净通量；DQ为重力势作用下向饱和带补给的通量；UQ为源汇项作用下饱和带向非饱和带补给的通量；DF为扩散作用下饱和带向非饱和带补给的通量；

根据饱和非饱和界面通量，采用地下水埋深波动法更新地下水埋深：

式中，D^t和D^t-1为第t和t-1时间步的地下水埋深；μ为给水度；

步骤6.分配饱和带与非饱和带之间的净通量；

当地下水埋深更新后，将饱和带与非饱和带之间的净通量分配到地下水埋深变动带，计算公式：

式中，θ^t为t时间步更新后的土壤含水率；θ_k ^t为t时间步内经过扩散作用计算后的土壤含水率；θ_s为饱和含水率；

步骤7.采用上述步骤3～6对各时间步下的饱和-非饱和水分运动参数均进行计算，得到计算区各时间步下各网格的运动参数值。

2.根据权利要求1所述的一维饱和-非饱和水分运动参数计算方法，其特征在于：

其中，运动参数包括：土壤水分、地下水埋深、非饱和带向饱和带的补给通量、饱和带向非饱和带的补给通量以及饱和带与非饱和之间的净通量。

3.根据权利要求1所述的一维饱和-非饱和水分运动参数计算方法，其特征在于：

其中，在步骤2中，采用最大埋深+1m的深度作为垂向计算厚度。

4.根据权利要求1所述的一维饱和-非饱和水分运动参数计算方法，其特征在于：

其中，在步骤2中，垂向网格划分厚度设置为0.05～0.20m之间。

5.根据权利要求1所述的一维饱和-非饱和水分运动参数计算方法，其特征在于：

其中，在步骤3中根据地下水埋深位置将地下水位所在分层分为2层，在一个时间步结束后，将网格调整回静态网格，合并计算过程中地下水埋深以上的第一层j和地下水埋深以下的第一层j+1为新的第j层，合并层的含水率为：

式中：

为合并后新的第j层的含水率；θ^j和θ^j+1为第j层和j+1层的含水率；M^j和M^j+1为第j层和第j+1层的厚度。

6.一维饱和-非饱和水分运动参数计算装置，其特征在于，包括：

资料获取部，获取研究区的气象资料、灌溉资料、土质资料、土壤含水率资料、地下水埋深资料；

静态网格划分部，根据计算时段内的最大埋深确定垂向计算深度，确定垂向网格划分厚度，并根据这些基础信息进行静态垂向网格划分；

动态网格划分部，根据地下水埋深位置在静态网格基础上增加一层；设当前正处于第t个时间步，上一个时间步计算得到的地下水埋深为D^t-1，找到地下水埋深所在的第j层，将第j层从地下水埋深处分开，上一层为最后一层非饱和带计算层，下一层为第一层饱和带计算层，这2层的底部距地表的距离为：

z^j，t＝D^t-1 (3-1)

Z^j+1，i＝z^j，t-1 (3-2)

式中，z^j,t和z^j,t-1为t和t-1时间步第j层的底部距地表的距离；z^j+1,t为t时间步第j+1层距地表的距离；D^t-1为t-1时间步的地下水埋深；

非饱和带模块运行部，在一个时间步内执行步骤4.1～4.4四个分步骤，更新非饱和带的土壤含水率，并得到饱和带与非饱和带之间的各项交换通量；具体计算方法见以下分步骤；

步骤4.1计算入渗水分配；

入渗水自上而下按饱和含水率进行填充，入渗水分配后，第j层的填充水量为：

式中，q^j为单位面积的第j层土壤所分配的水量；M^j为第j层土壤的厚度；

是第j层土壤的饱和含水量；θ^j是第j层土壤的含水率；I是单位面积的入渗水的总量；

是第j层以上土壤入渗过程中已经分配过的入渗水量；

步骤4.2计算重力作用下的土壤水分运动；

当土壤含水率超过田间持水率时，土壤水分在重力作用下向下层排水，其计算公式为，

式中，K(θ)为非饱和水力传导度，采用以下公式计算：

式中，θ为土壤含水率，K_s为饱和水力传导度，θ_f为田间持水率；田间持水率θ_f受地下水埋深的影响，采用以下公式计算：

θ′_f＝ωθ_f (4-4)

式中，θ′_f为修正后田间持水率；w为田间持水率修正系数；h_d为田间持水率受地下水埋深影响的高度；D为地下水埋深；

求解重力作用下的土壤水分运动方程后，可得到通过重力作用补给饱和带的水通量，记为DQ；

步骤4.3计算源汇项作用下的土壤水分运动；

在源汇项作用下，土壤水分运动方程如下：

式中，W是源汇项，包括土壤蒸发和作物蒸腾；当实际蒸发蒸腾的消耗深度比地下水埋深更深时，将会消耗地下水，产生向上的水通量，记为UQ；

步骤4.4计算扩散作用下的土壤水分运动；

采用含水率型方程描述的扩散作用下的土壤水分运动，其计算公式：

式中，D(θ)为扩散度，p_l(l＝1,…,L)是土壤水分特征曲线的参数，L为参数个数；

为描述土质垂向变异性的修正项；在扩散作用下，饱和带的水分将会补给至非饱和带，通过的水量记为DF；

饱和带模块运行部，在一个时间步内，经过入渗水的分配及在重力势、源汇项和扩散项作用下的土壤水分运动，计算得到在地下水面的交换通量，计算公式如下：

FQ＝DQ-UQ-DF (5-1)

式中，FQ为饱和带与非饱和之间的净通量；DQ为重力势作用下向饱和带补给的通量；UQ为源汇项作用下饱和带向非饱和带补给的通量；DF为扩散作用下饱和带向非饱和带补给的通量；

埋深更新部，根据饱和非饱和界面通量，采用地下水埋深波动法更新地下水埋深：

式中，D^t和D^t-1为第t和t-1时间步的地下水埋深；μ为给水度；

净通量分配部，当地下水埋深更新后，将饱和带与非饱和带之间的净通量分配到地下水埋深变动带，计算公式：

式中，θ^t为t时间步更新后的土壤含水率；θ_k ^t为t时间步内经过扩散作用计算后的土壤含水率；θ_s为饱和含水率；

参数值获取部，采用所述动态网格划分部、所述非饱和带模块运行部、所述饱和带模块运行部、所述净通量分配部对每一时间步下的饱和-非饱和水分运动参数均进行计算，得到计算区各时间步下各网格的运动参数值；以及

控制部，与所述资料获取部、所述静态网格划分部、所述动态网格划分部、所述非饱和带模块运行部、所述饱和带模块运行部、所述净通量分配部、所述参数值获取部通信相连，控制它们的运行。

7.根据权利要求6所述的一维饱和-非饱和水分运动参数计算装置，其特征在于，还包括：

输入显示部，与所述控制部通信相连，用于让用户输入操作指令，并进行相应显示。

8.根据权利要求7所述的一维饱和-非饱和水分运动参数计算装置，其特征在于：

其中，所述输入显示部能够根据操作指令对计算得到的饱和-非饱和水分运动参数按照列表或者图示方式进行显示。

9.根据权利要求7所述的一维饱和-非饱和水分运动参数计算装置，其特征在于：

其中，所述输入显示部还能够根据操作指令对计算得到的饱和-非饱和水分运动参数按照时间顺序以动态演化的方式显示在相应的研究区模型图上。

10.根据权利要求6所述的一维饱和-非饱和水分运动参数计算装置，其特征在于：

其中，所述运动参数包括：土壤水分、地下水埋深、非饱和带向饱和带的补给通量、饱和带向非饱和带的补给通量以及饱和带与非饱和之间的净通量。

Images