涌灌条件下非饱和土壤水分运动的数值模拟

张巧利，张 亮

(1.河南广播电视大学机电工程学院,郑州 450046；2.郑州职业技术学院，郑州 450000)

0 引 言

涌流灌溉是一种新型的适用于果树及植树造林的节水灌溉技术，它由滴灌演变而来，通过土壤根系层附近开挖一定深度灌水孔形成柱状辐射水流入渗，形成类似椭球体湿润体[1]。由于涌流灌溉入渗通常在地下进行，无污染，水分蒸发损失小，可以有针对性的湿润作物根部区域土壤，作物能够有效吸收并利用水分。涌流灌溉方式有效储水空间大，毛管直接深入灌水孔，一旦堵塞，容易发现并更换，从而有效解决了滴灌、渗灌在山地经济林灌溉时易堵塞和老化的问题，同时它对地形的适应性也很强。因此，涌流灌溉作为一种新型的节水方式，有一定的发展前景。研究不同灌水情况下土壤入渗湿润体特性有重要的意义。通过室内外试验可以获得土壤不同分布界面上水分运动状况，但由于试验的局限性无法对各种状况进行全面分析，且试验过程存在较大误差。近年来，以非饱和水运动理论为基础的数值模拟方法被用来定量分析微灌条件下土壤水分分布状况。计算机模拟方法简便、灵活和快捷，在确定的初始、边界条件下，可以模拟不同因素组合下土壤湿润体的水分状况，从而为合理确定灌水设计参数提供方便的研究手段。

国内外学者通过大量试验手段和数值模拟方法研究了土壤水分及溶质运移[2-6]。Cote等[7]人运用HYDRUS-2D软件模拟了滴灌入渗条件下土壤水分运移；许迪等[8]人通过建立点源土壤水分运动方程与溶质对流-弥散方程，并采用HYDRUS-2D求解了模型；周青云等[9]用有限差分法求解了蓄水坑灌条件下土壤水分入渗数学模型；依据非饱和土壤水分运动理论，李光永等[10-12]人建立了点源滴灌条件下土壤水分运动动力模型，并对模型进行求解；李明思[13]在圆柱体湿润体模型的基础上，建立了不同流量入渗下滴灌入渗的数学模型，并通过室内试验研究了不同滴头流量下各因子的变化过程，模型研究为滴灌系统的设计提供了重要的参考依据；本文依据非饱和土壤水动力学理论，试图建立涌流灌灌溉条件下土壤入渗数学模型，并通过室内试验对模型正确性予以验证。

1 涌流灌溉条件下土壤水分运动的数学模型

1.1 模型条件假定

模型都是在理想条件下提出的，与实际有一定出入，须作如下假定：假定土壤为固相刚性介质，各项同性，土壤水分不可压缩，忽略容质势、温度势作用；忽略土壤水分运动滞后效应；土壤基质势和土壤含水率的关系在土壤水分特征曲线中呈单值函数，土壤的密度(体积)和机械组成决定函数的形状。

1.2 土壤水分运动方程

水分运动边界为柱状界面，属于三维轴对称入渗，依据达西定律及质量守恒定律，通过以下偏微分方程来描述土壤水分运动：

(1)

式中：h为土壤总水势[L]；kx(θ)，ky(θ)，kz(θ)分别为土壤在x、y、z三个方向上的土壤非饱和导水率[LT-1]；θ为体积含水率[L3L-3]；t为计算时间[T]；x、y、z分别为三维坐标[L]。

1.3 初始及边界

(1)初始条件[14]。

h0(x,y,z)=aθ-b0，t=0

(2)

式中：h0为土壤剖面初始负压水头，cm；θ0为土壤初始体积含水量，cm3/cm3；假定模型计算区土壤水分均匀分布，通过选定的土壤水分特征曲线将已知的土壤初始含水率转化为需要模拟时初始输入的土壤负压剖面。其中a，b为土壤水分特征曲线常数。

(2)边界条件。图1为数值模拟计算区示意图。 均质土壤水分运动可以视为轴对称的三维入渗问题，假定计算区域较大，则水分在入渗过程中没有到达边界a′a″b″b′，a″c″b″无水量交换，为零通量界面。

图1 计算区示意图

则a′a″b″b′边界：

h0(x,y,z)=aθ-b0，t>0，x=xa，y=yb

(3)

a″c″b″边界为下边界，处理时考虑地下水埋深较大，看作第一类边界：

h0(x,y,za″c″b″)=h0，t>0″

(4)

a′a″c″c′与c′b′b″c″为对称边界，看作不透水边界，水通量为零。

(5)

a′b′c′边界为上边界，竖直方向水流通量为降雨强度或蒸发强度q0：

(6)

边界o′o″d′d为柱状入渗界面，由试验观测可知，有明显的饱和积水区，为运动的边界，较为复杂，一是灌水初期，灌溉出流速率小于土壤入渗速率时，灌水器所供应的水分很快渗入到土壤中，孔洞中未形成积水，土壤入渗速率等于边界流量，在模型中根据所模拟计算结果中累积入渗量来推求边界入渗速率；二是随着入渗的进行，孔洞周围土壤含水率达到饱和含水率，孔洞中形成积水，积水深度随时间呈渐进增加趋势，一定时间后达到稳定，此时边界按动水头边界处理计算。ho′o″d′d=h′(h′通过试验实测获得，试验时孔外安装测压管来观测入渗水头数据)。

1.4 模型求解

1.4.1 HYDRUS-3D

HYDRUS-3D是一款可以用来模拟土壤水分及溶质运移的有限元分析软件。模型中的水流状态为非饱和三维达西水流，计算过程忽略温度和空气对水流的影响，控制方程采用Richards方程，方程中可以通过嵌入汇源项来考虑作物根系对水分的吸收。程序可以灵活处理多种复杂边界，有固定流量边界、动水头和定水头边界、大气边界及排水沟等。计算区域可以由不规则或各项异性非均质土壤构成。水流计算区域可以通过三棱柱网格划分，采用伽辽金线状有限元法求解控制方程。时间离散对饱和及非饱和均采用隐式差分法。离散化后的非线性控制方程组可采用迭代法将其线性化。

1.4.2 模型求解

利用HYDRUS-3D软件求解以上所建模型[15]。模型计算区域为长方体，灌水孔置于中间，尺寸根据试验决定，采用棱柱体(1.25×1.25×1.25×2.5)对计算区进行划分，如图2所示。由于孔洞周围水势梯度变化较大，为了提高计算的精确度，该处适当对网格体加密。模拟时间采用变时间步长计算，初始时间步长定为1 min，最大和最小时间步长定为10和3 min。

图2 模拟计算流区的网格划分

2 试验验证

2.1 土壤水动力学特性

试验土壤土质为砂壤土，土壤容重为1.32 g/cm3，土壤稳定入渗率为1.5 mm/min，土壤平均含水率为10.4%。试验土壤取自河南农科院试验基地，物理性状为：粒径>0.02 mm占55.29%，0.02～0.002 mm占33.87%，粒径<0.002 mm占10.84%。

根据经验采用压力膜仪来测定土壤水分特征曲线，定水头法测定土壤饱和导水率。采用Van Genuchten模型对土壤水分特征曲线θ(h)和土壤导水率K(θ)进行拟合。

(8)

其中：

(9)

式中：以上方程中有5个独立参数θs,θr,α,n,K，其中：θs为土壤饱和含水率， cm3/cm3；θr为土壤残余含水率，cm3/cm3；m、n、α为模型的拟合参数，其中：m=1-1/n，α的取值与土壤物理性质有关，cm-1；Ks为土壤渗透系数，cm/min；Se为有效含水量，cm3/cm3；L为孔隙连通性参数，对大多数土壤一般取为0.5。供试土壤水分特性参数为：θr=0.028 3 cm3/cm3，θs=0.467 2 cm3/cm3，α=0.035 7 cm-1，h=1.48，Ks=0.125 cm/min。

2.2 测试方法

为了验证模型的正确性以及所选参数的合理性，采用室内土箱入渗进行物理模型试验。试验装置由土箱和供水系统两部分组成，土箱由厚度10 mm的有机玻璃组成，规格为80 cm×80 cm×100 cm(长·宽·高)。试验前将土壤风干过筛，按试验要求密度分层装土。供水采用马氏瓶，土箱中部孔洞周围用套管固定，管壁开孔，用纱网包裹，防止土壤进入孔洞堵塞出水孔，毛管伸入孔洞内进行灌水，灌水量通过开关进行调节。试验结束后在以孔洞为中心一侧对称断面不同位置取土测定含水率。

2.3 结果验证

模拟计算了地表以下不同位置处(观测点布置如图3)涌流灌溉下土壤中心断面含水率分布情况。图4为模拟最终值与实测结果在灌水量14 L结束后距灌水孔不同位置处土壤含水量对比情况。由于模拟计算条件均为理想状态，在初定的参数求得计算结果后与实测数据进行比对，利用实测数据反推模型中与土壤物理性质有关的参数，发现影响模拟结果的主要参数是m、α，通过修正参数计算结果，直到模拟值能较好地与实测值吻合为止。为了保证试验数据的准确性，同一组试验各做三次，取平均值作为最后的试验结果。图中对灌水量为14 L灌水420 min情况的模拟与实测结果对比分析表明，不管是水平方向或者垂直方向，灌水孔周围水势梯度变化大的点处土壤水含水率计算值与实测值比较接近，表明用数值模拟方法可以较好地表示湿润体内水分分布状况。

图3 观测点位置图

图4 土壤含水率模拟与实测值对比

表1为土壤入渗420 min后距离灌水孔不同位置观测点处含水率分布状况实测与模拟对比情况，观测点布置如图3所示。通过实测值与模拟值对比验证，误差定义为：(实测值-模拟值)/实测值。结果表明，相对误差基本保持在±10%以内，两者一致性较好。

3 模型的应用

通过模型对以上土壤参数条件下地下涌流灌溉土壤水分状况进行初步模拟分析。

3.1 土壤湿润体内水分分布状况

作为一种新型的节水灌溉方式，掌握涌灌土壤湿润体内含水量分布规律，对于科学合理施行作物水分管理，合理布置监测仪器有重要的意义。试验实测获得土壤水分在不同土质、流量、初始含水率下的动态分布状况，耗时费力，且试验过程繁琐。通过数值模拟可以较好地解决以上问题，且模拟能够为延伸试验提供可靠的结果。图5为灌水量为14 L，灌水420 min情况的入渗湿润体内不同监测点处土壤含水率及负压水头随时间的动态变化情况，从中可以得知，靠近灌水孔处土壤含水量最大，增加速率最快，当土壤含水量接近土壤饱和含水量时，增加速率明显减缓，甚至不再增加。计算结果能连续反映土壤水分动态分布规律。

模拟土体具有相同的土壤初始含水率。从模拟得到的结果可以得出，灌水孔周围的土壤水压力几乎接近正压，含水率越接近饱和含水率。距离灌水孔水不同位置处对土壤水分的影响主要是湿润锋到达的时间不同，含水率增长速率随着远离灌水孔越小，且最终含水率也越小。

图5 涌流灌溉下入渗剖面土壤模拟水分分布状况

3.2 土壤湿润体分布

在砂壤土条件下分别对灌水流量为2.0、3.6、5.4 L/h 3种情况进行模拟计算，分析灌水流量对土壤水分分布状况的影响，其中，模拟的土体土壤初始含水率及灌水量均一致。图6给出了3种不同流量下灌水结束后土壤水分等值线图及土壤含水量分布，由图6可知，灌水量相等，随着流量的增加，湿润圈形状呈水平方面束窄，垂直方向伸长的变化趋势，这是由于土壤水分在重力势及基质势的共同作用下运动，流量增大后，土壤湿润体内部含水率迅速增大，灌水孔中积水深度也变大，所以湿润体上半部分扩大较快；土壤含水率较大的范围增加，在土壤湿润体内部水分分布上，相同位置处土壤含水率随流量的增大而增大。综上可知，灌水流量较大时，土壤中靠近灌水孔附近局部土壤含水率在较短时间内达到较高值，从而水势梯度变大，促进土壤水分扩散。但是流量过大会破坏周边土壤结构，在来不及扩散的情况下孔洞内产生积水，向上入渗湿润体扩散至地表，造成水分蒸发损失。所以选择流量时，要充分考虑所选流量既不会破坏周围工作条件亦不会发生深层渗漏和扩散到地表造成水分蒸发损失。

图6 不同流量湿润体含水率等势线图(单位：cm3/cm3)

4 结论与讨论

(1)本文依据非饱和土壤水动力学理论，建立了涌流灌溉条件下土壤水分运动数学模型，通过采用HYDRUS-3D对所建模型进行模拟，模拟结果采用室内试验实测数据进行了验证，通过实测数据对模型参数进行修正，发现参数α、m对模型计算影响较大。最后经试验数据修正以后的参数代入模型计算，结果与实测数据吻合较好。

(2)应用所建模型分析了地下涌流灌溉土壤水分运动的一些规律，从而可以为涌流灌溉系统的合理设计和运行提供理论指导。

(3)由于涌流灌溉系统流量一般较大，根系吸收水分对土壤水分运动有一定影响，在今后的研究中有待进一步增加根系吸收模型。

[1] 吴普特，朱德兰.涌泉根灌技术研究与应用[J].灌排机械工程学报，2010，28(4)：355-357.

[2] Thorburn P J,Cook F J,Bristow K L.Soil-dependent wetting from trickle emitters:implications for system design and management[J].Irrigation Science,2003,22(3/4):121-127.

[3] Provenzano G.Using HYDRUS-2D simulation model to evaluate wetted soil volume in subsurface drip irrigation systems[J].Journal of Irrigation and Drainage Engineering,2007,133(4):342-349.

[4] Zhou Q Y,Kang S Z,Zhang L,et al.Comparison of APRI and Hydrus-2D models to simulate soil water dynamics in a vineyard under alternate partial root zone drip irrigation[J].Plant Soil,2007,291(1/2):211-223.

[5] 赵伟霞，蔡焕杰，陈新明，等.无压灌溉土壤湿润体含水率分布规律与模拟模型研究[J].农业工程学报，2007，23(3)：7-12.

[6] 康银红，马孝义，李 娟，等.黄土高原重力式地下滴灌水分运动模型与分区参数研究[J].农业机械学报，2008，39(3)：90-95.

[7] Cite CM.Bristow K L,Charlesworth P B.et al.Analysis of soil wetting and solute transport in subsurface trickle irrigation[J].Irrigation Science.2003，22(3/4)：143-156.

[8] 许 迪，程先军. 地下滴灌土壤水运动和溶质运移数学模型的应用[J]. 农业工程学报，2002，18(1)：27-31.

[9] 李光永，郑耀泉，曾德超，等. 地埋点源非饱和土壤水运动的数值模拟[J]. 水利学报，1996，(11)：47-56.

[10] 周青云，孙西欢，康绍忠. 蓄水坑灌条件下土壤水分运动的数值模拟[J]. 水利学报，2006，37(3)：342-353.

[11] 李光永，曾德超，段中锁，等. 地埋点源滴灌土壤水运动规律的研究[J]. 农业工程学报，1996，12(3)：66-71.

[12] 李光永，曾德超，郑耀泉. 地表点源滴灌土壤水运动的动力学模型与数值模拟[J]. 水利学报，1998，(11)：1-25.

[13] 李明思，康绍忠，孙海燕. 点源滴灌滴头流量与湿润体关系研究[J]. 农业工程学报，2006，22(4)： 32-35.

[14] Warrick A W. Soil water dynamics[M]. Oxford: Oxford University Press, 2003.

[15] Simunek J，Sejna M, van Genuchten M T. The HYDRUS-2D Software Package for Simulating the Two-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media[M]. California: U.S. salinity laboratory, Agricultural Research Service, U.S.Department of Agriculture, Riverside, 1999.