基于Hydrus-2D模型的膜下滴灌对降雨入渗影响研究

(辽宁水利土木工程咨询有限公司，辽宁 沈阳 110003)

入渗(又称为下渗)，是水从地表垂直或水平渗入到土壤和地下的运动过程，这是自然界中水文循环必不可少的一部分。Hydrus是一个可用来模拟土壤水流及溶质三维运动的有限元计算机模型，它可以灵活处理各类水流边界，包括定水头和变水头边界、给定流量边界、渗水边界、大气边界以及排水沟等。

入渗是半干旱地区水分循环的重要环节，以内蒙古东部平原区这一半干旱地区土壤在膜下滴灌方式下的降雨入渗过程为例，本文研究了Hydrus-2D模型在土壤水分入渗过程模拟中的应用，以期为模型在该地区的应用及推广提供借鉴。

1 试验区概况

试验区位于内蒙古东部平原区，在通辽市开鲁县，为西辽河流域，属大陆性温带半干旱季风性气候，年平均气温6.80℃，年平均降雨量340.5mm，无霜期148d，试验区属于典型的北方半干旱区，水资源匮乏，供水主要以地下水为主，而用水主要以农牧业为主，随着农业生产规模的不断扩大和气候的变化，致使该地区的水资源短缺现象日益严重。

2 试验区土壤条件

试验点位于东经121°23′，北纬43°37′，试验点2m 深度范围内土壤分为3层，其中第1层0～30cm，第2层30～60cm，第3层60～200cm。用环刀在各层取土样，同时，用自封袋和铝盒在对应各土层取土样，每层3个土样重复操作。将用环刀取的土样带回实验室测定土壤容重，将用铝盒取的土样带回实验室测定土壤饱和含水率，将自封袋中的土样风干、碾压、筛分、均匀混合后用筛分法测定土壤粒径(各测定结果见表1)。

表1 实验区土壤基本物理性质

3 基于Hydrus-2D模型的膜下滴灌方式下的降雨入渗规律研究

Hydrus是用于模拟饱和与非饱和多孔介质中水分、能量、溶质运移的一维、二维和三维有限元计算模型软件。该模型的水流状态为饱和-非饱和达西水流，水流控制方程采用Richards方程，溶质运移方程采用对流-弥散(CED)方程，模型方程求解采用Galerkin线性有限元法。软件可以灵活处理各类水流边界，包括定水头和变水头边界、给定流边界、渗水边界、自由排水边界、大气边界等。另外，该软件也包括一个参数优化算法，可用于各种土壤的水压和溶质运移参数的逆向估计，为使用者提供了反推土壤特性经验参数的途径，从而提高了拟合效率和拟合效果。

3.1 建立模型

3.1.1 水分运动基本方程及模型

(1) 降雨入渗模型方程

因为膜下滴灌和地面灌的降雨入渗过程可以看成土壤水分在水平方向和垂直方向的运移过程，因此，本文基于Hydrus-2D软件，建立膜下滴灌和地面灌二维降雨入渗模型。根据Richards方程，并考虑植物根系吸水的二维土壤水分运动方程为：

(1)

式中θ——土壤体积含水率，cm3/cm3；

t——时间，h；

k(θ)——土壤非饱和导水率，cm/h；

φ——土壤总水势，cm，对于非饱和土，φ=φm±z，其中φm为土壤基质势；

z——观测点垂直坐标，坐标向上为正，向下为负；

S——汇项，表示根系吸水速率，cm/h。

(2) 土壤水分运动模型

Hydrus模拟涉及土壤水分特征曲线和土壤导水率曲线两类参数，本文对这两类参数的刻画釆用vanGenuchten-Mualem模型(Van-Genuchten 1980)：

(2)

(3)

式中θ(x)——土壤物理特征曲线；

θs——饱和含水率，cm3/cm3；

θr——剩余饱和率，cm3/cm3；

h——压力水头，cm；

K(h)——压力水头下的导水率；

Ks——饱和导水率；

α、m、n和l——经验常数。其中m和n之间有一定的相关关系，一般认为m=1-1/n时模型的简化效果最好(Van-Genuchten 1980)。对于大多数土壤，l参数一般取0.5；

Se——土壤饱和度。

(3) 根系吸水模型

在根系吸水模块，Hydrus软件共提供两种模型，Feddes模型和S-Shape模型，本次模拟采用Feddes模型，计算公式为:

S(h)=α(h)b(x)Tp

(4)

式中S(h)——根系吸水速率，即从单位体积土体单位时间内吸取的水量；

α(h)——水分胁迫响应函数；

b(x)——根系吸水分配密度函数；

Tp——潜在蒸腾量。

3.1.2 模型概化

(1) 时间和空间离散化

对试验区进行原位观测的时间内，7月份降雨量最大，为81.10mm，其中21日及25日两天降雨量分别为35.90mm、25.40mm，基于本文的研究目标，选取这两天的降雨量进行模拟计算。土壤模拟时期分为两段：第一段自2016年7月20日至2016年7月24日，共120个小时；第二段从 2016年7月24日至 2016年7月29日，共120个小时，各组模拟的时间都以h为单位，每隔2h输出一个计算结果，用以与实测值进行比较。模拟土层为0～2m的垂直剖面。将一维的土壤剖面平均剖分为201个节点，节点间距为1cm。在上述土层剖分的基础上，根据土壤剖面特征和实测的土壤质地、颗粒组成、容重等数据，将土层分成若干介质层，以便赋予不同的土壤水分运动参数(同一介质层的土壤水分运动参数相同，各层之间的土壤水分运动参数不同)。进行土壤水分动态模拟时将土壤剖面分为0～30cm、30～60cm、60～200cm3层。

(2) 作物根系概化

玉米根系在土壤中的分布与生长方式可分为水平根和垂直根。根据一些学者研究发现玉米根系在水平方向延展距离一般为0.8m，在竖直方向延展距离一般为1.25m。玉米在不同质地的土壤中根系分布各不相同，但总体上玉米大部分根系主要分布在0～40cm土层内并在生育后期深层土壤根系量有所增加。在本次模拟研究中设定地表以下0～40cm深度范围内为主要根系吸水区。

(3) 观测点的设置

因为膜下滴灌存在着覆膜和起垄，而它们的存在影响了农田土壤水的入渗、蒸发、运移等过程，因此为了更好地研究该灌溉方式下降雨入渗过程，可以把膜下滴灌分为3个土壤剖面进行研究，分别为膜中剖面、膜边剖面和沟中剖面，同时，在各个剖面设置观测点，共设置30个观测点，平均分布在膜中、膜边和沟中3个剖面上，即每个断面分配10个观测点，其中在土壤深度0～40cm范围内每隔10cm设置1个观测点，在40cm深度一下每隔20cm设置1个观测点。

3.1.3 初始及边界条件的设定

a.初始条件：Hydrus软件提供了负压水头和土壤含水率两种初始条件形式，本次模拟选取了土壤含水率初始条件形式。因为模拟时期分为两段，所以在进行模拟时以降雨前一天作为每一次模拟的初始时间，即 2016年7月 20日和 2016年7月24日，给定当天的整个剖面各个位置的实测土壤含水率作为初始条件，相邻测点间的含水率按线性插值程序自动给出。每一次模拟又分成膜下滴灌组和地面灌组并分别给出初始条件。

b.边界条件：根据实际情况边界条件设定如下：ⓐ上边界，将膜下滴灌覆膜的部分设定为无水流边界，而无膜部分土壤与大气直接相通，既有降雨的入渗又受到蒸发作用，所以将其设定为大气边界。地面灌方式下的土壤与大气直接相通，所以其上边界也设定为大气边界；ⓑ下边界，研究区域内年平均地下水埋深为9m，大于所模拟的2m深度，故膜下滴灌和地面灌下边界都设定为自由排水边界。

c.大气边界为时变边界条件，主要包括降雨量、蒸发量、蒸腾量等信息，其中降雨数据根据实测降雨资料输入，蒸发状况则由模型根据输入的气象条件及植被情况自行计算。

3.2 模型参数的确定

在此模型模拟中将土层划分为3层进行模拟。根据土壤的颗粒组成和容重，利用Hydrus软件自带的神经网络预测模块，首先对土壤水分运动参数设定初始值，然后求解正问题，并将结果与实测值进行比较，反馈并修改土壤水分运动参数，直至模拟值与实测值差距最小，这时得到的土壤水分运动参数认为是最优的结果(各参数最终结果见表2)。

表2 土壤水分运动参数

3.3 运行模型及拟合效果分析

为直观评价模型的模拟效果，下面给出了模拟时间段内各层土壤含水率实际观测值和对应模拟值的对比关系图(如图1所示)。

由图1可看出：ⓐ在观测深度下土壤含水率的实际值和模拟值的变化趋势大体相同；ⓑ在土壤含水率较高时的模拟值与实测值偏差相对较小；ⓒ膜下滴灌三个剖面的模拟效果比膜下滴灌整体层面上的模拟效果相对较好。总的来说，模拟值能很好地反映出土壤含水率随土层深度变化时的变化趋势，所建模型对土壤含水率的模拟效果较好，可用于试验区田间土壤水分运动的模拟。

3.4 模型应用

降雨入渗过程主要受到降雨强度、初始含水率、下垫面方式、土壤特性等因素的影响。对于相同的研究对象其土壤特性一般保持不变，所以在研究同一区域降雨入渗时主要考虑另外三个影响因素。

3.4.1 模拟试验方案设置

模拟试验设置为两因素多水平试验，两因素分别为：降雨强度、初始含水率。根据试验区多年降雨资料，降雨强度分别设置为2mm/h、4mm/h，8mm/h3个水平；所观测的土壤深度范围内共有3种土质，共分3层，结合2016年试验所测得的含水率数据，3层土壤初始含水率分别设置为对应层土壤饱和含水率的40%、60%和80%。试验因素水平见表3。

表3 试验因素水平

3.4.2 模型模拟与求解

本次模拟试验是在已建立好的膜下滴灌和地面灌模型基础上进行的，在模拟时仅对降雨条件和初始含水率条件进行改变，其他信息不做改变。所模拟的一次累积降雨量为32mm，当雨强为2mm/h时设置降雨历时为16h，连续降雨，中间无间隔；当雨强为4mm/h时设置降雨历时为8h，连续降雨，中间无间隔；当雨强为8mm/h时设置降雨历时为4h，连续降雨，中间无间隔。在进行各次模拟试验时设置模拟总时长都为120h，初始时刻设置为降雨开始前4h。

模型计算中，初始步长为0.5h，最小步长为0.001h，最大步长为1h，每隔2h输出一个计算结果。

3.4.3 模拟结果分析

因为降雨入渗基本上在雨后1—2d内就能完成，所以在分析时取雨后48h输出的模拟结果，通过对模拟结果进行整理，绘制出不同降雨强度、不同初始含水率条件下，膜下滴灌在雨后48h时的土壤各层含水率状况图(如图2所示)。

图2 不同初始含水率条件下雨后48h时的土壤含水率分布图

由图2可知：ⓐ在降雨强度不变的情况下，初始含水率越低，雨后经过一段时间其土壤含水率就越大，这种现象在浅层尤其明显，原因是初始含水率越低其土壤入渗能力越强，在同等的降雨强度下，入渗量就越多，从而含水率也就越大；ⓑ在初始含水率不变的情况下降雨强度越低，降雨后土壤含水率在整体上就越大，并且这种现象在浅层最为明显，由此可见土壤入渗能力随着降雨强度的增加而降低。

为进一步了解初始含水率和降雨强度对土壤入渗的影响，做出了模拟试验雨后48h的累积入渗量分布图(如图3所示)。

图3 不同雨强、不同初始含水率条件下雨后48h累积入渗量分布图

由图3可知：初始含水率和降雨强度对降雨累积入渗量都有影响，具体表现为：在降雨强度一定的情况下，初始含水率越高降雨累积入渗量越低；在初始含水率一定的情况下，降雨强度越大降雨累积入渗量相对来说就越低。

4 结 语

本文以开鲁县田间降雨入渗过程为例，介绍了 Hydrus-2D 软件组成，基于室内试验取得的各项土壤物理特征参数，模拟了膜下滴灌方式下土壤含水量变化过程，并以部分实测数据对模拟的入渗过程进行了验证，从而实践了Hydrus-2D 模型在土壤水入渗过程中的应用，模型能很好地模拟出膜下滴灌方式下降雨入渗过程，可用于实验区田间土壤水分运动的模拟计算。

试验结果表明：在降雨强度不变情况下，初始含水率越低，雨后经过一段时间土壤含水率就越大；在初始含水率不变的情况下降雨强度越低，降雨后土壤含水率在整体上就越大；在降雨强度一定的情况下初始含水率越高，降雨累积入渗量越低；在初始含水率一定的情况下降雨强度越大，降雨累积入渗量相对来说就越低。