滴灌双点源土壤水分入渗特性数值模拟研究

李耀刚，蒋存仁

（1.甘肃天水市水务局，甘肃天水741000；2.甘肃省水利水电勘测设计研究院，兰州730000）

滴灌作为一种现代高效节水灌溉技术在农田节水灌溉中被广泛使用，它可以根据田间不同作物需水要求进行灌溉。滴灌属于局部微灌溉，通过频繁、缓慢地将滴头处消能后的少量水浸润作物根系层土壤，在时间与空间上能够比较精准地调控土壤含水量，同时，土壤地表水分无效蒸发和深层渗漏损失大幅减小，水分利用效率高达90%以上，相较其他农田灌溉方式节水效益明显[1]。作物根系经由滴灌入渗形成的湿润体获取水分，根系分布状况及生长状态一定程度上受土壤湿润体水分分布的影响。滴灌条件下土壤水分运动规律的研究可为滴灌系统设计中合理确定灌水持续时间、灌水流量及滴头布置间距提供重要的理论依据。点源滴灌入渗土壤水分分布规律及湿润体特性被国内外学者进行了大量试验研究和数模模拟，得出了一些比较成熟的理论和方法[2-4]。张志刚[5]等研究了不同滴头流量对地上滴灌土壤水分分布的影响，得出土壤湿润锋的入渗距离比值、运移速率与灌水时间呈幂函数关系。张松[6]等通过室内试验研究了地下点源滴灌入渗土壤水分入渗规律，初步研究了地下点源滴灌在沙土中的适应性。数值计算方面以Richards方程为基础的点源入渗数值模型被提出，裴青宝[7]等对多点源入渗条件下土壤水分溶质运移进行了模拟；李斯[8]运用Hydrus模型模拟了土壤有效水边界对滴灌入渗土壤湿润体特性的影响；许迪等[9]通过建立的非饱和土壤水分运动和溶质运移数值模型模拟了地下单点源滴灌入渗中水、肥分布规律；Skaggs[10]等对地上、地下滴灌土壤水分运动过程应用Hydrus-2d软件进行了模拟。前人对点源滴灌入渗试验及数值模拟进行了许多有益的探索，然而在农田灌溉中，对于株间距较小的农作物如棉花、玉米等，为使每株植物都能从土壤中获取水分，滴灌入渗就需要在湿润区形成湿润带，而多点源入渗在滴头连线方向土壤湿润体呈带状分布，因此多点源滴灌入渗湿润体交汇情况下土壤水分运动的研究具有重要的现实意义。本文在上述研究的基础上，以双点源滴灌条件下的土壤水分分布为研究对象，根据非饱和土壤水动力学理论，利用软件HYDRUS-3D建立了物理和数值模型，利用模型计算输出结果分析了双点源滴灌滴头流量和灌水时间对土壤水分运动特征的影响，以期为进一步研究多点源滴灌条件下土壤水分运移规律提供理论依据。

1 数学模型

1.1 基本方程

忽略气象因素、容质势、温度势对土壤水分运动的影响，假定土壤各项同性且为刚性介质，土壤中的水分不可压缩，不考虑土壤水分滞后效应。三维等温均匀达西水流水分运移方程用修改过的Richards方程：

式中：ψ 为土壤总水势， cm；θ 为土壤体积含水率，cm3/cm3； k(θ)为土壤非饱和导水率， cm/min；t 为计算时间，min，x、y、z分别为横纵轴向坐标，cm。

1.2 定解条件

1.2.1 初始条件

滴灌模拟区域几何边界如图1 所示。尺寸为：100 cm×100 cm×60 cm（长×宽×高)，模拟土壤介质选择沙壤土，容重为1.35 g/cm3，假定土壤初始剖面含水率（负压水头）均匀分布，为方便计算，本文将土壤含水率通过土壤水分特征曲线转换为负压水头，其值为0.085 cm3/cm3(负压水头为-1 200 cm水柱）。求解方程的初始条件为：

式中：h0为土壤初始负压水头，cm。X、Y、Z分别为模拟计算区域的横向、纵向及垂向最大距离。本文中，X=100 cm，Y=100 cm，Z=60 cm。

1.2.2 边界条件

上边界：

（1）滴头处：

式中：q为滴头流量，L/h；E为土壤表面蒸发强度，cm/min。

（2）滴头之外：

（3）地表积水，饱和区h=0；非饱和区：

下边界：h=h0，z=Z，0≤t≤T。

左右边界：假定边界对称，考虑边界上水分通量为0：

1.3 模型求解

采用隐式差分法离散求解方程式（1）：

式中：i、j 分别为距离、时间节点编号；j+1、j 为当前和前一计算时刻；k为迭代次数。

应用HYDRUS-3D对建立的数值模型进行求解运算。计算区域采用有限元网格进行剖分，通过GALERKIN 有限元法求解。滴头处水势梯度大，通过增加计算节点加密网格，满足计算精度要求。初始、最小和最大时间步长分别按0.01、0.001、10 min 计；运算结果每5 min间隔输出一次；最大迭代次数为10 次；压力水头、含水率容差分别为1、0.001。初始条件采用土壤初始负压水头，假定计算区域土壤初始含水率均匀分布，计算区内负压水头均为-1 200 cm。边界条件采用连续流输入，计算区下边界设为自由排水面。为分析湿润体内土壤水分变化情况，在计算区距离滴头水平、竖直间隔5 cm处设置20个观测点。

2 试验验证

2.1 土壤水动力学特性

采用沙壤土作为供试土样，土壤容重为1.35 g/cm3，土壤平均含水率为0.08 cm3/cm3。供试土样物理性状：粒径大于0.02 mm 的沙粒含量占比55.29%，粒径介于0.02～0.002 mm 的粉粒含量占比33.87%，粒径小于0.002 mm 的黏粒含量占比10.84%。土壤水分特征曲线通过压力膜仪测定拟合，土壤饱和导水率通过定水头法测定。土壤导水率和土壤水分特征曲线通过van Genuchten模型进行拟合：

式中：5 个参数θs，θr，α，n，Ks独立；θr为土壤残余含水率，cm3/cm3；θs为土壤饱和含水率，cm3/cm3；α、m、n 为拟合参数，其中α 的大小与土壤物理性质密切相关，cm-1；L 为孔隙连通性参数；Se为土壤有效含水量，cm3/cm3；Ks为土壤渗透系数，cm/min。

本文中，土壤水分特性参数拟合值分别为：θs=0.45 cm3/cm3，θr=0.018 cm3/cm3，α=0.009 cm-1，n=1.6，Ks=0.006 cm/min。

2.2 结果验证

图2 给出了沙壤土在滴头流量为1.2 L/h，灌水时间为300 min时距滴头水平和垂直方向不同位置处土壤含水率的实测与模拟结果对比。从含水率柱状图可以看出，实测与模拟结果比较接近，水平距离土壤含水率相对误差分别为5.74%、3.38%、4.42%、4.21%；垂直距离土壤含水率相对误差分别为0.68%、3.36%、5.09%、8.45%、4.94%，均在10%以内。因此，湿润体内土壤水分分布的数值模拟能较好地反映土壤水分分布的实际情况。

3 双点源滴灌湿润体内土壤水分分布模拟结果与分析

利用建立的双点源入渗模型模拟了灌水时间、滴头流量对滴灌土壤水分分布规律的影响，对湿润体内土壤水分分布规律进行了初步探索研究。

3.1 灌水时间对湿润体内土壤水分分布的影响

土壤水分入渗湿润体大小的变化规律，是正确进行双点源滴灌系统田间布置和灌溉水分管理的重要依据。沿滴头布置方向的横断面可以反映双点源灌水条件下的土壤水分交汇过程，因此本文主要对该断面上的土壤水分分布规律进行分析。图3反映了滴头流量为1.2 L/s，灌水持续时间分别为250、300、350、400、500、600 min 时，沙壤土土壤水分分布的模拟过程。坐标（0、0）、（30、0）为滴头所在位置。从图3 可以看出，入渗开始后，湿润锋未交汇前，土壤水分入渗遵循点源入渗，湿润体内土壤含水率等值线剖面形状近似呈椭圆形。土壤水分随时间推移由滴头向外不断扩散，土壤含水率等值线图离滴头由近及远的位置上呈由疏到密分布。灌水时间为250 min 时，点源入渗湿润锋水平距离小于15 cm，垂直入渗深度大于15 cm，湿润锋没有交汇，滴头处土壤含水量达到土壤饱和含水量。当灌水时间超过268 min 后，湿润体开始融合交汇。随着入渗时间推移，湿润体形状逐渐由2个分离的半椭圆变成花生壳形状，除靠近滴头处土壤含水率较高的区域土壤含水率等值线图近似于1/4 椭圆形外，远离滴头的湿润区内土壤等值线图融合为一条曲线，并且随着时间的推移，曲线向下移动，等值线形状从屋脊形过渡为水平。交汇区内下部土壤含水率等值线分布平缓，上部分布弯曲。灌水时间增加，滴头处土壤水分饱和区范围逐渐增大，交汇面土壤含水率增加。

分析双点源滴灌入渗条件下湿润体内的土壤含水率变化对于科学合理地施行作物水分管理，制定合理的滴灌灌溉制度有一定的指导意义。图4 所示滴头流量为1.2 L/h 灌水600 min 情况下的距离滴头不同位置观测点（测点布置见图5）的模拟土壤含水率变化情况，由于滴头对称布置，湿润体对称分布，故选取一侧的观测点模拟数据进行分析。由图4 可和，灌水开始后，滴头处（N1）土壤含水率迅速升高，10 min 内由0.08 cm3/cm3增加到0.42 cm3/cm3，随后增长速率趋于平缓。由测点N4、N3、N2可知，当时间分别为25、110、268 min时水分运移至距滴头水平方向5、10、15 cm 处；当时间分别为25、90、250、480 min 时水分运移至距滴头正下方竖直方向5、10、15、20 cm 处。土壤水分平均水平运移速率分别为0.2、0.09、0.05 cm/min，平均垂直运移速率分别为0.2、0.11、0.06、0.04 cm/min。垂直方向土壤入渗速率高于水平方向，这是由于土壤水分在垂直方向入渗受到基质势和重力势共同作用。随着时间推移水势梯度减小，土壤水分在水平和垂直方向平均入渗速率均减小。当湿润锋开始交汇时，N2 的土壤含水率迅速升高，灌水结束后土壤含水率由0.08 cm3/cm3提高到0.3 cm3/cm3。灌水结束时地表N4点的含水率小于0.35 cm3/cm3，地表以下N5的含水率大于0.35 cm3/cm3，N6的含水率大于N3，可以看出，滴灌点源入渗过程中，距滴头水平和垂直相同距离测点含水率有差距，垂直方向含水率高于水平方向，随着水分扩散，水平和垂直方向土壤含水率的差距逐渐减小。这是因为在滴灌入渗过程中基质势是土壤水分水平运动的主要动力，而土壤水分向下运动的动力除基质势外还有重力势。N7 的含水率小于N2，因为湿润锋的交汇先从地表开始，接着逐渐向下进行交汇，交汇时，交汇面处土壤含水率明显大于下层土壤。由测点N9、N11、N14可以看出，土壤水分运移到该点的时间分别为60、225、480 min。

3.2 流量对土壤湿润体内水分分布的影响

对双点源土壤入渗不同滴头流量供水情况下的土壤湿润体内含水率分布状况，采用试验实测的方法来获得土壤含水率分布剖面相对比较困难，而采用数值模拟的方法可以得到可靠的计算结果。图6 为2 种滴头流量（1.2、2.4 L/h）、灌水时间（250、600 min）组合情况下的湿润体内土壤含水率等值线剖面。由图6 可知，土壤含水率等值线图类似于湿润体形状，土壤水分由滴头向外扩散，水势梯度减小，等值线图由疏变密。灌水时间相同，滴头附近各节点处土壤含水量，土壤湿润体剖面面积，高含水区面积及同一位置湿润区的土壤平均含水率均随灌水流量的增大而变大。其它条件不变，流量增加，湿润锋运移速率加快，水平和垂直距离增加，交汇时间变短，反之亦然。地表因滴头流量增加而形成积水区，湿润锋水平运移速率随着积水区半径增大而加快，滴头流量增加使得湿润锋与积水区边界间水分梯度减小。说明滴头流量对湿润体形状和土壤含水率的分布均有影响。

交汇面含水率分布状况的研究对合理制定灌溉制度有一定指导意义。图7所示为灌水结束后，地表滴头连线中心交汇面纵剖面含水率变化情况。q=2.4 L/h 时，距地表面竖直向下5 cm 处含水率为0.38 cm3/cm3；0～15 cm 范围内剖面含水率较高，高于0.35 cm3/cm3，未达到土壤饱和含水率；土壤表层向下20 cm 以内土壤含水率缓慢减小，由0.38 cm3/cm3减小到0.35 cm3/cm3；距土壤表层向下20 cm 以下，含水率减小速率增大，竖直向下20 ～30 cm 范围内，含水率由0.35 cm3/cm3减小到0.08 cm3/cm3。q=1.2 L/h 时，距地表面竖直向下3 cm 处土壤含水率最大，为0.31 cm3/cm3；0～7 cm 范围内含水率超过0.3 cm3/cm3，距土壤表层向下7 cm 范围内土壤含水率缓慢减小，由0.31 cm3/cm3减小到0.30 cm3/cm3；距土壤表层向下7 cm 以下，含水率减小速率增加，竖直向下7～18 cm 范围内，含水率由0.3 cm3/cm3减小到0.08 cm3/cm3。滴头流量对湿润体内含水率的变化有一定影响。

4 结 语

（1）本文在非饱和土壤水动力学理论的基础上试图建立了双点源滴灌土壤水分运动数值模型，应用流体建模有限元分析软件HYDRUS-3D 进行模拟，其结果与实测数据进行对比，一致性较好。在理想状态下，可以通过数值计算方法来探索研究双点源土壤水分入渗影响因素。

（2）滴灌双点源入渗土壤含水率的分布同时受灌水时间和滴头流量的影响。当滴头流量q=1.2 L/h，灌水时间为268 min时，湿润锋开始由点源入渗过渡为双点源交汇入渗，湿润体形状逐渐由2个分离的半椭圆变成花生壳形状，除靠近滴头处土壤含水率较高的区域土壤含水率等值线剖面近似于1/4 椭圆形外，远离滴头的湿润区内土壤等值线图融合为一条曲线，并且随着时间的推移，曲线向下移动，形状从屋脊形过渡为水平。随着灌水时间延长，滴头处土壤水分饱和区范围逐渐增大，交汇处土壤含水率增加。灌水时间相同，湿润锋运移速率、滴头附近各节点处土壤含水量，土壤湿润体剖面面积，高含水区面积及同一位置湿润区的土壤平均含水率均随灌水流量的增大而变大。滴头流量、灌水时间对双点源滴灌土壤水分入渗均有一定影响。