地下水浅埋下层状土壤波涌畦灌间歇入渗模型研究

陈 琳 费良军 傅渝亮 钟 韵

(1.西安理工大学水利水电学院, 西安 710048; 2.省部共建西北旱区生态水利工程国家重点实验室, 西安 710048)

0 引言

波涌灌溉条件下的土壤水分入渗属于间歇入渗，国内外学者在其节水机理、入渗特性、水流运动特性和数值模拟，以及波涌灌适应性、灌溉实施方案及灌溉效益分析方面，对均质土下波涌灌溉间歇入渗理论和技术进行了大量研究，这些研究主要在土壤质地相对简单的均质土进行[1-8]。在实际农业生产中，由于气象、水文、地质和生物过程的作用，土壤并非简单的均一物质，大都呈现为交错分布的层状结构。如在冲积平原地区，由于河流的砂粘交错沉积，形成复杂多变的土壤剖面结构，导致田间土壤大都呈现为层状结构，其土壤孔隙结构的不均匀性使得水力特性在土层面处发生了突变，因而与均质土连续入渗相比，对波涌灌条件下的层状土间歇入渗土壤水分运动特性的变化研究显得更加复杂[9]。根据土壤水能量原理，水在层状结构的土体内入渗时，无论夹层土壤的质地粗或细，土体内的夹层均会对下渗水流起到阻水的作用[10-15]。目前，对波涌灌条件下的夹砂层层状土壤间歇入渗水分运移特性研究较少，仅王丁[16]对不同含砂层埋深下波涌灌技术要素对层状土壤间歇入渗水分分布特性进行了初步研究和分析，且未对影响层状土入渗时水分运动参数的确定进行分析，对土壤水进入夹砂层时稳渗率的确定也仅局限于连续入渗过程。因此，本文针对层状土壤条件开展波涌灌间歇入渗水分运动特性及影响参数的研究，对进一步发展波涌灌灌溉理论与技术、提高水分利用率、减少水分流失具有重要的理论价值和实际意义。

1 材料与方法

1.1 试验装置与方法

试验在西安理工大学西北旱区生态水利工程国家重点实验室进行。采用傅渝亮等[5]开发的间歇入渗设备进行灌水试验，并观测地下水埋深下间歇供水入渗量及剖面水分动态。试验装置由土柱、地下水控制系统和入渗供水系统组成，如图1所示。

试验土柱由高210 cm、内径15 cm的有机玻璃管制成。土柱底部60 cm深度内装粒径为3～7 mm的砂砾石模拟地下饱和含水层，上部150 cm为非饱和层状土壤，两部分间垫以弹性较差的薄海棉或过滤纸，以免上层土壤进入砂砾层。沿土柱垂向设4排取土孔，孔径1.5 cm，用橡皮塞封堵，孔中心间距为5 cm，相邻两排取土孔以垂直孔距2.5 cm交叉布置。土柱上部设加水斗，供水开始时先向加水斗供水，以免水流冲击破坏土柱表土致密层，并且在供水停水时便于尽快排净土柱上表面水层。土柱下部20 cm范围内，每5 cm设排水口一个，共4个，试验前排水口用胶皮管和夹子封堵，避免漏水。

图1 试验装置示意图Fig.1 Schematic diagram of experimental apparatus1.马氏瓶1 2.加水口 3.夹砂层 4.土层 5.粗砂层 6.排水汇流管 7.地下水供水口 8.马氏瓶2 9.取土孔 10.进水口

1.2 供试土样

试验土壤采用粉壤土，夹砂层采用砂壤土，土样经过风干、碾压、过筛、配水得到试验土样，初始土壤饱和导水率采用南-55型渗透仪测定，土壤pH值采用定性pH试纸测定，观测过程中，入渗量通过马氏瓶进行控制和计量(单位：cm)；土壤采用小型土钻采取土样，用干燥法测定土壤含水率。土壤颗粒分析使用英国马尔文公司的Mastersizer-2000型激光粒度仪测定，颗粒组成和初始理化参数如表1所示。

1.3 试验方案

室内试验模拟地下水位为150 cm条件下，肥液质量浓度为200 mg/L，周期数i为3、循环率r为1/2，停水时间为40 min，含砂层埋深为20 cm，夹砂层厚度D为5 cm，并以均质土处理作对比，进行土壤入渗特性室内试验。每个处理3次重复，其观测指标取平均值。由马氏瓶2供给地下水，待上升毛管水分布稳定之后，进行间歇入渗试验。

表1 试验土壤颗粒组成和初始理化参数Tab.1 Particle gradation and initial physicochemical parameters

2 理论分析

2.1 湿润锋面水吸力与湿润深度的关系

当地表下存在埋深较浅的地下水位时，则地下水位以上土壤剖面的土壤含水率呈非线性分布。设地下水位为H，则

H=φ(θ)+z

(1)

式中z——位置水头，cm

φ——位置水头z处相应的基质吸力，cm

θ——含水率，cm3/cm3

为了确定土壤剖面初始土壤含水率分布，确定体积土壤含水率与基质吸力之间的关系，采用BROOKS等[6]所提出的土壤水分特征曲线拟合公式

(2)

其中

式中K(θ)——非饱和导水率，cm/min

Θ——相对土壤含水率

θr——土壤滞留含水率，cm3/cm3

φb——土壤进气吸力，cm

n——形状系数，反映土体孔径的分布特征

l——弯曲度，一般取1[18]

结合式(1)和式(2)，浅埋地下水条件下的均质土土壤初始含水率分布计算式可写成

(3)

则累积入渗量I可以表示为

I=(θs-θi)zf

(4)

式中zf——湿润锋运移距离，cm

由质量守恒定律可知，式(4)可写为

(5)

将式(3)代入式(5)，则有

(6)

变形得

(7)

其微分形式为

(8)

由Green-Ampt入渗模型可知

(9)

式中Sf——湿润锋面处水吸力，cm

h——压力水头，cm

q——地表水分通量(入渗率)，cm/min

结合式(7)和式(9)，并代入式(8)变形可得

(10)

(11)

当Sf确定时，便可以利用五阶Runge-Kutta法[19]求解微分方程式(11)中时间t以及相应湿润锋运移距离zf。

根据BOUWER[4]研究结果可将概化的湿润锋面水吸力Sf表示为非饱和土壤导水率的函数，即

(12)

式(12)并非理论推导得出，MEIN等[20]则进一步提出用土壤水吸力的加权平均值来表示，即

(13)

其中

Kr=K(θ)/Ks

式中Kr——相对渗透系数

MOORE等[21]将上述湿润锋面处概化水吸力Sf修正为

(14)

通过比较分析，可采用式(14)计算Sf，则联立式(2)和式(3)并代入式(14)可得

(15)

若φb、n、H一定，则式(15)为Sf与湿润锋运移距离zf为z时的函数关系式。其中式(15)的参数φb、n可根据不同土壤的性质进行确定。显然，当土壤剖面的初始土壤含水率θi沿深度方向为一常量时(均质土)，则式(15)中的H-z为一定值，即初始含水率θi对应的φi为一常数，此时对应的Sf也为常量。

2.2 层状土连续入渗模型

对于连续入渗过程中的层状土的入渗变化特性，入渗过程可分为非线性自由入渗与线性稳渗两个阶段[22]，第1阶段：在入渗锋面未到达层状土上层之前，实际上是一个均质土的积水入渗问题，其入渗量I随时间的变化均为一个非线性变化关系，可根据式(7)、(10)、(15)分别对连续入渗的zf、I与t进行求解，随着入渗湿润锋的下移，当土壤湿润锋剖面达到层状土上层时，由于阻水作用，导致湿润锋暂时停止运移，但湿润锋位置处的土壤含水率逐渐增加，相应基质水吸力逐渐降低，直到达到砂土的进水吸力时，湿润锋才开始穿过含砂层上表面并进入到含砂层，入渗过程迅速进入了第2阶段——稳渗阶段，稳渗率q*可改写为

(16)

式中q*——稳渗率，cm/min

K′(θs)——上层壤土近似土壤饱和导水率，cm/min

对于参数K′(θs)的确定，范严伟等[13]通过含砂层土壤Green-Ampt入渗模型的改进与验证结果发现，引入导水度系数Cw来量化上层壤土达到接近饱和时的饱和导水率K′(θs)与饱和时Ks之间的关系，可用K′(θs)=CwKs表示，并通过分析得到适宜的Cw平均值为0.95。

对于再次穿过含砂层下层并进入到壤土中的湿润锋运移距离变化定量计算，可对式(11)进行修正，因湿润锋运移速率实际受到上层含砂层稳渗率和壤土质地共同影响，其中式(11)的Ks值不应是其饱和导水率，而由含砂层确定了的稳渗率q*控制，则水分由含砂层进入壤土阶段的湿润锋模型修正式为

(17)

2.3 层状土间歇入渗模型

(18)

2.4 参数确定

2.4.1水分运动方程

假定土壤为均质、各向同性的多孔介质，忽略温度与土壤中的气相对土壤水分的影响，不考虑根系吸水与源汇项。考虑到研究饱和-非饱和流动问题，在直角坐标系对垂直积水入渗的水分运动数学模型以负压水头h作为基本方程因变量，其中，模拟过程中的土壤水力特征参数选择BROOK-COREY(B-C)[6]模型，其表达形式见式(2)。控制方程Richards方程为

(19)

式中C(h)——比水容重，cm-1

hf——土壤吸力，cm

2.4.2初始条件与边界条件

取地表为零基准面，坐标轴(z轴)向上为正。模拟剖面水流模型可概化为：大气蒸散量、降水入渗量取零。孙西欢等[23]、王志涛等[24]研究表明，压力水头在0～10 cm之间对累计入渗量的影响较小；另一方面，为达到短时间内及时排空土柱表面积水，避免因表面积水排出时间过长影响间歇阶段脱失过程的进行。故确定上边界定压力水头供水hg为2.5 cm；间歇停止时，上边界定通量边界q为0。下边界为负压水头，始终为0。模拟区域为：Z0≤z≤0，其中地下水埋深Z0为-150 cm。模拟供水时间ton为40 min，循环率r为1/2，周期数i为3。

第1周期供水阶段：

初始条件为

h(z,t)=h(z,tG1)=hG1(z) (Z0≤z≤0，t=tG1=0)

(20)

上边界条件为

h(z,t)=h(tG1)=hG(z=0，0

(21)

下边界条件为

h(z,t)=hb(Z0,tG1)=hb(z=Z0)

(22)

第1周期间歇阶段：

初始条件为

h(z,t)=h(z,tJ1)=hJ1(z) (Z0≤z≤0，tJ1=ton)

(23)

上边界条件为

(24)

下边界条件为

h(z,t)=hb(Z0,tJ1)=hb(z=Z0)

(25)

第i周期(i=2，3)供水阶段：

初始条件为

h(z,t)=h(z,tG1)=hGi(z) (Z0≤z≤0,tGi=0)

(26)

上边界条件为

(27)

下边界条件为

h(z,t)=hb(Z0,tGi)=hb(z=Z0)

(28)

第i周期(i=2，3)间歇阶段:

初始条件为

h(z,t)=h(z,tJi)=hJi(z) (Z0≤z≤0,tJi=ton)

(29)

上边界条件为

(30)

下边界条件为

h(z,t)=hb(z0,tJi)=hb(z=Z0)

(31)

式中hG1(z)——第1周期剖面初始土壤吸力，cm

hJ1(z)——第1周期供水结束时剖面土壤吸力，cm

hGi(z)——第i-1周期间歇阶段结束时剖面土壤吸力，cm

hJi(z)——第i周期供水阶段结束时剖面土壤吸力，cm

hb——地下水埋深处土壤吸力，cm

hb(t)——随时间变化的地下水埋深处土壤吸力，cm

tJ1、tG1——第1周期间歇阶段停水持续时间和第1周期供水阶段供水持续时间，min

tc——周期时间，min

2.4.3数值求解方法

Hydrus-1D模型[25-26]主要用于变量饱和多孔介质的水流和溶质运移。本文以各供水周期的土壤含水率剖面实测数据作为输入项，结合Hydrus-1D模型自带的Inverse反推模块，采用Levenberg-Marquardt迭代(LM算法)反推最优土壤水分运动参数，当压力水头允许偏差小于0.1 cm，且土壤含水率容许偏差小于0.001 cm3/cm3时，迭代结束，便得到各供水周期内的土壤水分运动参数最优解。模拟时段按间歇入渗周期进行设定。采用变时间步长划分方法，设定初始时间为0 min，进一步调整模拟时间步长，本文设定时间最小步长为0.001 min，最大步长为0.005 min。

3 结果与分析

3.1 模型验证

3.1.1稳渗率的确定

表2 稳渗率实测值与简化模型计算值比对Tab.2 Comparison of steady infiltration rate between measured and calculated values

3.1.2土壤水分运动参数的确定

利用Hydrus-1D对基于Brook-Corey模型下的均质土和层状土各供水周期对应的土壤水分运动参数θr、θs、α(α=1/φb)、n值进行求解，最终得到各供水周期间歇入渗土壤水分运动参数，见表3。由方差分析结果可以看出，层状土间歇入渗过程中的粉壤土与均质土间歇入渗过程中的土壤水分运动参数基本一致，其滞留土壤含水率θr和土壤形状系数n值随周期数变化的差异性均不明显(P>0.05)，但饱和土壤含水率θs、进气吸力及饱和导水率Ks均随周期数改变呈极显著性差异(P<0.01)，相比均质土，砂壤土(含砂层)仅对进气吸力随周期数的变化呈极显著性差异(P<0.01)，而对其他参数的影响差异性均不显著(P>0.05)。

表3 入渗阶段Brooks-Corey模型参数模拟值Tab.3 Parameter simulation value of Brooks-Corey model in infiltration stage

注：同一列数据后不同小写字母表示不同周期数的差异显著；** 表示p<0.01水平上差异极显著。

利用式(7)、(16)、(18)分别将均质土及层状土的各周期内累计入渗量计算值和实测值对比，并将基于电容充电理论模拟计算的入渗量同时与实测值

表4 供水周期稳渗率实测值与模型计算值对比Tab.4 Comparison of steady infiltration rate between measured and calculated values in each cycle

进行对比，结果如图2、3所示。图2为均质土第1、第2、第3供水周期基于Brook-Covey结合Green-Ampt模型(简写为BC-GA模型)计算的入渗量和基于电容充电理论模型计算的入渗量和实测值对比图，图3为层状土第1、第2、第3供水周期基于BC-GA模型计算的入渗量和基于电容充电模型计算的入渗量与实测值对比图。

由图2、3可看出，两个入渗模型对均质土和层状土计算的入渗量与实测值的变化趋势基本一致，BC-GA模型与电容充电模型计算的入渗量相比，整体更接近于实测值。进一步给出了两种方法的计算值与实测值对比结果。图4a表示基于BC-GA模型计算值与实测值之间的对比，图4b表示基于电容充电模型拟合值与实测值之间的对比。

图2 均质土各供水周期入渗量变化Fig.2 Changes of infiltration with time of homogenous soil under all water supply cycle

图3 层状土各供水周期入渗量变化Fig.3 Changes of infiltration with time of layered soil under all water supply cycle

图4 入渗量计算值、模拟值与实测值对比Fig.4 Comparison of calculated, simulated and measured infiltration volumes

从不同模型对比结果可以看出，根据BC-GA模型模拟的入渗量与实测值相比，其决定系数R2为0.993，相应平均绝对误差MAE为0.038 cm，均方根误差RMSE为0.065 cm，而根据电容充电模型拟合的入渗量计算值与实测值相比，其决定系数R2也达到0.973，相应平均绝对误差MAE为0.070 cm，均方根误差RMSE为0.122 cm，说明不论均质土还是层状图，两种模型均可以定量地反映波涌畦灌间歇入渗实际入渗量随时间变化规律，但与电容充电入渗模型相比，BC-GA模型更接近于实际情况，虽然基于BC-GA模型拟合度较高，但取决于Hydrus-1D反推迭代的好坏，特别是在模拟过程中，对土柱离散步长还有时间离散步长的选取是提高BC-GA模型估算精度的关键。

图5为基于BC-GA模型湿润锋运移距离随时间变化计算值。可以看出，与均质土相比，在湿润锋运移到含砂层上边界之前，属于均质土自由入渗阶段，因此湿润锋运移距离随时间变化规律应与均质土基本一致，当湿润锋运移距离达到含砂层上边界时，因阻渗作用提前达到稳渗状态，湿润锋运移距离陡然减小，周期数越大，与均质土湿润锋运移距离差距越大，当达到第2、第3供水周期时，层状土湿润锋面运移距离均未超过30 cm，而均质土相同供水周期的运移距离分别达到51 cm和73.5 cm，说明层状土间歇入渗可显著抑制湿润锋运移距离的发展。

图5 层状土与均质土各供水周期湿润锋计算值对比Fig.5 Comparison of wetting front calculated values between layered soil and homogeneous soil in each water supply cycle

为对比层状土不同供水周期对含砂层的入渗特征影响，将各周期湿润锋到达和穿过夹层界面的时间和累积入渗量列于表5。可以看出，到达含砂层上边界所经历的时间t1随周期数的增大逐渐减小，达到含砂层下边界所经历的时间t2同样随周期数的增大逐渐减小，结合图2、3和表5可以看出，湿润锋面运移距离由含砂层上边界到达含砂层下边界所经历的历时(Δt=t2-t1)随着周期数的增加显著减小，对应的累积入渗量变化量(ΔI=I2-I1)同时也减小。

3.2 模型应用

表5 各供水周期湿润锋到达土层界面的时间及对应累积入渗量Tab.5 Time and cumulative infiltration when wetting frontreached interface of layers under each water cycle

注：“已穿过”代表供水周期开始时的湿润锋运移位置大于含砂层上边界或下边界；“未穿过”代表供水周期结束时的湿润锋运移位置小于含砂层上边界或下边界。

为对比不同含砂层埋深对供水周期内间歇入渗的变化特性，将各供水周期不同含砂层埋深的湿润锋到达和穿过夹层界面的时间和累积入渗量列于表7。

表6 各入渗周期不同含砂层埋深对应的稳渗率Tab.6 Steady infiltration rate of each supply water cycle under different sand buried depths

表7 不同含砂层埋深处理供水周期湿润锋到达土层界面的时间及对应累积入渗量Tab.7 Time and cumulative infiltration when wetting front reached interface of layers under different sand buried depths in each cycle

图6 各供水周期湿润锋运移距离变化Fig.6 Variation of distance of wetting front with time under all water supply cycle

由表7可看出，相同供水周期内，无论是含砂层上边界还是下边界，到达对应边界位置的时间随着埋深的增加呈滞后的趋势，即埋深越深，到达上、下边界所需时间越久，第1供水周期，埋深15、20 cm在供水周期时间内均已穿过含砂层上边界，但25 cm并未到达含砂层上边界z1，因此含砂层水分未发生改变。但3个埋深水平在供水周期时间内均未穿过含砂层下边界z2；当达到第2供水周期时，到达含砂层上、下边界的历时时间要比第1供水周期的显著缩短，仅埋深25 cm入渗结束时未达含砂层下边界；当达到第3供水周期时，3个埋深水平均穿过含砂层上边界和下边界，且历时时间最短。

图7 各供水周期累积入渗量变化Fig.7 Variation of infiltration with time under all water supply cycle

根据BC-GA模型计算结果，图6为不同含砂层埋深下间歇入渗的湿润锋运移距离随时间变化关系对比图，图7为不同含砂层埋深下间歇入渗的累积入渗量随时间变化关系对比图。可以看出，相同含砂层埋深条件下，随着周期数的增加，层状土整体湿润锋运移距离随之增大，不同的是当供水周期结束时，不同含砂层埋深之间的湿润锋运移相对距离对比差异性变化随周期数的增加而逐渐减小，当达到第3供水周期结束时，不同含砂层埋深处理之间的湿润锋运移距离分别达到27.96 cm(埋深为25 cm)、27.03 cm(埋深为20 cm)、26.55 cm(埋深为15 cm)，含砂层埋深越浅，湿润锋运移距离越小。

相对于湿润锋，不论是均质土还是层状土，相同埋深条件下的累积入渗量整体随着周期数的增加而增大，而相同周期内，含砂层埋深越浅，其累积入渗量越小，当达到第3供水周期结束时，不同埋深下的累积入渗量基本接近平均值0.42 cm，差异不明显，说明在周期数和含砂层埋深两者共同影响下，当周期数不大时，埋深不同对累积入渗量差异性影响明显，当第3供水周期结束时，由于含砂层埋深导致的累积入渗量的差异性基本不大，含砂层埋深对累积入渗量差异性影响较小。

3 结论

(1)将层状土间歇入渗过程分为两阶段描述，即自由入渗和稳定入渗，建立了层状土间歇入渗各供水周期内入渗量与时间变化的函数关系，不论均质土还是层状土，拟合后的决定系数均高于0.97，均方根误差均在0.13 cm以内，符合实际情况。

(2)通过理论分析，结合Brook-Corey和Green-Ampt两模型共同建立了基于地下水浅埋影响下的层状土(均质土)间歇入渗模型，推导出各供水周期湿润锋面处水吸力与湿润锋运移深度的函数关

系式，并求得层状土各供水周期穿过含砂层时的稳渗率，通过与实测值对比验证，表明所建立模型和求解方法适用于描述波涌畦灌层状土间歇入渗过程。

(3)在周期数和含砂层埋深共同影响下，周期数增加，相同埋深下的累积入渗量越小，而湿润锋运移距离越大；随着含砂层埋深的增加，相同周期内穿过含砂层历时时间越久，对应的累积入渗量计算值越大，当达到第3供水周期结束时，含砂层埋深的不同对其湿润锋运移距离和累积入渗量变化差异性影响不如第1、第2供水周期结束时明显，因此周期数越大，含砂层埋深对其湿润锋运移距离及累积入渗量的影响越不明显。

(4)由于试验数据在室内理想环境下得到，且层状土间歇入渗水头为静水水头入渗，而实际大田土壤质地更为复杂，灌溉时灌溉水流始终处于流动推进状态，传统估算入渗量和湿润锋运移距离不易观测，且估算方法费时费力。可将观测点不同时刻土壤剖面含水率分布与Hydrus-1D反推模块相结合，间接求得层状土间歇入渗不同供水周期土壤入渗水分运动参数，并通过BC-GA模型估算出各供水周期的入渗量和湿润锋运移距离，为波涌灌技术的应用提供理论依据。