精准滴灌农田排水过程数值模拟

麦荣幸， 康跃虎， 裴 亮， 李晓彬

(1.中国科学院 地理科学与资源研究所 陆地水循环及地表过程重点实验室, 北京 100101； 2.中国科学院大学, 北京 100049)

1 研究背景

随着水资源的日益紧缺和现代科学技术的发展，精准灌溉是未来农业灌溉的发展趋势。滴灌作为现代灌溉技术，可根据作物需水需肥规律进行适时适量施肥灌溉，可以最大限度地减少蒸发、深层渗漏和地表径流等带来的水分和养分的损失，是实现精准灌溉的主要技术之一[1-2]。Wan等[3]、Kang等[4]和Wang等[5]在对滴灌农田水分养分循环过程、作物耗水规律、作物生长与耗水关系等进行了长期、系统、深入研究的基础上，提出了监测滴头正下方20 cm深度土壤基质势下限的滴灌农田水盐调控方法，为精准滴灌的实现提供了简单、可靠的方法。

根据土壤基质势精准滴灌农田条件下具有少量高频灌溉的特点，土壤水分得到高度控制，农田土壤整体处于水分非饱和状态，并且在时间和空间上均保持一个较为稳定的关系[6-8]。只在关键生育期灌溉的传统地面灌溉农田条件下，每次灌溉时土壤表层为饱和状态，长时间的灌溉间隔期内受降雨、作物耗水、蒸发等影响，土壤水分在时间和空间上表现出不稳定、变化剧烈的关系[9]。因此，与传统地面灌溉农田相比，精准滴灌农田的土壤水分运动与再分布过程、灌溉或降雨产生的深层渗漏量、灌溉产生的弃水、降雨产生的地表径流等均不相同，造成滴灌农田的排水过程和排水量也与传统地面灌溉农田有所不同[10]。随着滴灌在灌区农业中的大面积应用，以及灌区管理的现代化、信息化发展，亟需对精准滴灌农田的排水过程进行研究，为农田排水系统的规划设计与改造、灌区灌排信息化管理提供参考。

滴灌农田排水过程分为土壤水分入渗至浅层地下水和随地下水运动流入排水系统两个过程。目前国内外学者对滴灌条件土壤水分在非饱和带的运动、区域浅层地下水运动以及土壤水与地下水运动耦合方面做了大量的研究[11-13]，构建起了单点源和多点源的土壤水分在非饱和带的运动模型，即LINKFLOW、MODFLOW、FEFLOW等地下水运动模型以及SWMS-3D等耦合模型[14-17]。综观现有研究成果，对于精准滴灌农田条件下土壤水分随灌溉水、降水入渗至浅层地下水再运动到排水系统的整个排水过程的研究还未见到报道。本文在监测土壤基质势指导滴灌灌溉条件下，根据农田土壤水分分布特征，综合考虑浅层地下水运动的尺度问题，通过构建饱和带-非饱和带水分运动耦合模型，建立精准滴灌农田排水过程系统模型，对精准滴灌农田排水过程进行准确高效地数值模拟，对现代灌区排水系统规划及信息化管理具有重要的应用价值和意义。

2 精准滴灌农田排水过程的数学模型

随着灌溉和降雨，土壤水分垂向运动入渗至浅层地下水中，在水力梯度作用下随浅层地下水流入到排水沟中。理论上，农田中离排水沟越远的区域，地下水埋深越浅。这就造成了在排水过程中，非饱和带厚度在空间上的变化，不同位置的深层渗漏水量也不一样，排水过程也不相同。在非饱和带中，土壤水分除了存在垂向交换之外，还存在水平交换，水分在非饱和带的运动可采用三维Richards方程。在农田中，饱和带含水层厚度较大，浅层地下水向排水沟的运动主要以水平运动为主。水分在饱和带的运动可采用基于三维Richards方程沿饱和含水层垂向积分平均得到的水平二维方程[18]。在农田排水过程中，非饱和带通过深层渗漏水量补给饱和带。通过深层渗漏水量将非饱和带水分运动方程与饱和带水分运动方程耦合起来，构建精准滴灌农田排水过程系统模型。

2.1 模型的控制方程

2.1.1 非饱和带水分运动控制方程 在滴灌条件下，当滴灌带和滴头间距较小时，滴灌带间湿润锋会发生重叠，形成湿润带，土壤水分运动为三维运动。假设土壤为均质、各向同性的刚性多孔介质，不考虑温度、生物或化学反应对水分运动的影响，也不考虑土壤水分运动的滞后效应，土壤非饱和带水分运动可用三维Richards方程描述:

(1)

式中：K(h)为非饱和导水率，cm/min；θ为土壤体积含水率，cm3/cm3；h为土壤负压水头，cm；t为时间，min；x、y、z为空间坐标，x、y为水平方向坐标，z为垂直方向坐标，cm；规定z向下方向为正。

2.1.2 饱和带水分运动控制方程 在饱和带中，地下水流运动为三维运动，地下水流运动可用以 为因变量的Richards方程描述:

(2)

式中：Ks为土壤饱和导水率，cm/min；μ为比储水系数;I为源汇项。

在农田排水过程中，饱和带含水层的厚度较大，地下水流运动主要以水平方向运动为主，垂向运动较小。假设饱和带含水层隔水底板为水平，将公式(2)方程沿着饱和带含水层厚度垂向积分即可平均得到水平二维方程。

公式(2)左边的积分项为：

(3)

(4)

公式(2)右边的积分项为：

(5)

综合以上各式，饱和带地下水运动方程可采用以下水平二维方程形式：

(6)

2.1.3 饱和带和非饱和带的耦合项 在农田排水过程中，非饱和带通过深层渗漏水量补给饱和带，通过深层渗漏水量将非饱和带水分运动方程与饱和带水分运动方程进行耦合。公式(6)中的V|Z1便为非饱和带补给饱和带的深层渗漏水量，其值可用公式(1)从地表Z0沿着垂向积分至潜水水面Z1，并利用Darcy定理求得：

(7)

将公式(7)代入公式(6)中，可得农田排水过程系统模型：

(8)

2.2 定解条件

2.2.1 初始条件

h(x,y,z,t)=h0(x,y,z)

(9)

(Z0

h(x,y,t)=hs(x,y)

(10)

(Z1

2.2.2 边界条件

(11)

(z=Z0,t≥0)

h(z,t)=0 (z=Z1,t>0)

(12)

(13)

式中：E为作物蒸散发强度，cm/min；q=q0/A，q0为滴头流量，mL/min,A为土壤湿润面积，cm2；Ω为计算区域； ∂Ω为水平四周边界；n为(x,y)∈Ω的外法向量；g(x,y)为水平向排水沟流动的水流通量。

2.3 模型的求解

2.3.1 离散方程构建 利用Du Fort-Frankel有限差分格式对非饱和带水分运动方程和饱和带水分运动方程进行空间和时间离散，建立矩阵求解方程。采用矩形网格对空间进行划分，非饱和带含水率在垂向上急剧变化，需要较细的空间划分网格；饱和带空间尺度较大，采用较大的空间划分网格。

解 法 四 ：由f（x）=sin2x+a cos2x= （2x+φ）的周期为π，并于直线x=-π对称，可知f8，则，即a=-1.

2.3.2 程序编写和计算 计算程序编写采用的是MATLAB软件。首先，对非饱和带进行了网格划分和离散，形成非饱和带有限差分矩阵方程；再通过对饱和带进行网格划分和离散，形成饱和带有限差分矩阵方程，根据耦合关系形成饱和带-非饱和带统一矩阵方程，最后进行求解。具体计算流程见图1。

图1 程序设计流程图

首先输入节点、步长、定解条件和参数等信息，形成饱和带和非饱和带水流矩阵方程。假设浅层地下水初始水头为h1，确定非饱和带垂向上的节点数，利用非饱和带求解矩阵求得非饱和带土壤含水率分布，再利用饱和-非饱和统一矩阵方程求得当前时刻浅层地下水水头h2，若h2与h1收敛，说明h2为当前时刻水头值；若h2与h1不收敛，则将h2重新赋值给h1，重复上一步骤，不断迭代，直到h2与h1收敛为止，然后进行下一时刻的运算，直到运行结束。

3 模型验证

3.1 实测数据获取

在天津静海建立了精准滴灌农田排水过程浅层地下水位的田间观测系统。所选区域为45 m×95 m大小的四方形地块，地块的三周边界为排水沟渠，沟渠宽2 m；另一边界为水泥道路。所选地种植冬小麦，在地块上布置滴灌带，滴灌带间距0.6 m，滴头间距0.3 m，滴头流量1.38 L/h。所选地块及滴灌带布置见图2。

3.1.1 土壤基本特性 田间观测区域土壤类型为黏壤土，土壤的非饱和导水率K(h)与压力水头h为高度非线性函数，可采用Van-Genuchten[19]模型所给出的表达式。土壤的饱和导水率Ks采用环刀法进行测定，比储水系数 采用田间试坑法测定。具体参数见表1。

图2 所选地块及滴灌带布置示意图(单位：m)

土壤类型θr/(cm-3·cm-3)θs/(cm-3·cm-3)α/cm-1nKs/(cm-3·min-1)μ黏壤土0.0580.6360.00891.5770.00470.042

注：θr为土壤残余含水率；θs为土壤饱和含水率；α为土壤进气值的倒数；n为形状系数；Ks为土壤饱和导水率；μ为土壤比储水系数。

3.1.2 观测孔布置 在精准滴灌灌溉条件下，灌溉水分在整个地块上均匀分布，因此在相同土壤地质条件下，整个地块土壤的地下水分运动是轴对称的，可选地块的一半作为观测区域进行观测点布置。横向、纵向分别均匀布置10行(列)观测点，共布置100个观测点。

所选地块地下水位观测孔布置方式及编号如图3所示。

图3 所选地块地下水位观测孔布置示意图(单位：m)

观测孔编号：从下到上为1～10，从左到右为1～10，即左边第1列的最下方的点编号为NO 101，从下往上为NO 101～NO 110；右边第1列的最下方的点编号为NO 1001，从下往上为NO 1001～NO 1010；中间各点编号以此类推。

3.1.3 滴灌灌溉频率和灌溉水量 根据Kang等[6～7]提出的滴灌农田精准水盐调控方法，将滴头正下方20 cm深度处的土壤基质势阈值下限定为-20 kPa进行灌溉时，随灌溉水进入土壤中的盐分得到有效淋洗，土壤盐分平衡，作物获得高产。因此在整个地块内选择3个有代表性的点分别在距离滴头正下方20 cm深度处埋设负压计，当负压计显示的土壤基质势低于-20 kPa时，对地块进行灌水，每次灌水量为10 mm。

3.1.4 观测内容 每天早上7:00记录观测孔潜水水位、排水沟水位和冠层上20 cm蒸发皿水量蒸发量。在非饱和带水分运动过程中，需要考虑作物的蒸散发对水分运动的影响。目前测定作物蒸散发量的方法主要有水量平衡法、梯度法、蒸渗仪法、波文比法和涡度相关法等，这些测定方法所需参数较多，对参数要求精度高[20]。在田间实际观测过程中，由于实验条件有限，难以应用以上方法获取日尺度上的作物实际蒸散发量。Liu等[21]通过在华北地区对冬小麦蒸散发量的研究发现作物冠层顶部蒸发皿的蒸发量与作物实际蒸散发量有良好的经验关系，两者的回归系数趋近于1，根据冠层顶部蒸发皿的蒸发量能够较为准确地估算冬小麦的日蒸散发量。在已有的实验条件下，本文采用该方法对作物实际蒸散发量进行估算。

在拔节到灌浆期，冬小麦蒸散发量与冠层上20 cm蒸发皿水量蒸发量的关系式为[21]:

ETc=1.124ET20

(14)

式中：ETc为冬小麦蒸腾蒸发量,mm;ET20为冠层上20 cm蒸发皿蒸发量,mm。

图4为冠层上20 cm蒸发皿蒸发量在观测期间内随时间的变化过程。

图4 冠层上20 cm蒸发皿蒸发量随时间的变化

3.2 模拟结果与实测结果的对比与分析

将精准滴灌农田排水过程系统模型潜水水位模拟值与实测值进行比较，初步识别和验证模型。田间地下水位观测一共持续了30 d，系统模型模拟运行时间分为两个时间段：第3 d到第20 d和第21 d到30 d。因为在田间观测的第3 d和第21 d有强降雨，降雨量分别为62和48 mm。

3.2.1 灌溉周期内浅层地下水位的变化和模拟结果 任取田间观测1个灌水周期内的潜水水位与模拟水位值进行对比分析。图5为灌溉周期内浅层地下水位实测值变化过程，图6为灌溉周期内浅层地下水位模拟值变化过程，图7为灌溉周期内浅层地下水位模拟值与实测值之差的变化过程。其中田间观测的第7、8、9 d正好为1个灌溉周期，且没有降雨。从图5、6中可以发现，实验区地块靠公路端中部区域潜水水位最高，并沿着向排水沟方向逐渐降低，即离排水沟越近，潜水水位越低。潜水水位降低梯度沿着农田中间区域至排水沟方向逐渐减小，农田中间区域水位降低速度快，而排水沟渠附近区域水位降低速度慢。灌溉周期内的第1 d水位最高，并且随着时间的推移，水位逐渐降低且水位降低速率逐渐减小，即灌溉后随时间推移农田各点水位逐渐趋于稳定。结合图5和6，系统模型潜水水位模拟值与实测值在时间和空间上的变化趋势一致。

图7中，系统模型浅层地下水位模拟值与实测值之差很小，随时间推移，两者的误差由正值逐渐向负值转化，误差绝对值有增大趋势。中间区域水位误差大于两侧区域水位误差，受边界条件约束，离排水沟渠越近，水位模拟值越趋近于实测值。表2、3、4为灌溉周期内部分观测孔潜水水位模拟值与实测值的具体数值。表2、3、4中潜水水位模拟值与实测值的绝对误差在-13～21 mm范围内，潜水水位模拟值与实测值吻合效果较好。

3.2.2 降雨后浅层地下水位的变化和模拟结果 取NO 101、NO 501、NO 505这3个地下水位观察孔的潜水水位模拟值与实测值进行拟合分析。这3个观测孔分别位于地块公路端中间、中部中间和中部侧边，能一定程度上反映整个农田的模拟情况。图8为NO 101、NO 501、NO 505观测孔浅层地下水位模拟值与实测值的拟合结果。从图8中发现，第2次降雨后的潜水水位下降速度明显高于第1次降雨后的下降速度，这是因为第2次降雨后的作物蒸散发和沟渠水面蒸发速度高于第1次降雨后的蒸发速度。在模拟初始阶段，浅层地下水位模拟值与实测值拟合效果。

图5 浅层地下水位实测值(图中展示的是整个实验区域的一半，左边为整个研究区域的中轴线，下图类同)

图6 浅层地下水位模拟值

图7 浅层地下水位模拟值与实测值之差

很好，随时间推移模拟值与实测值的误差有增大趋势。部分原因是模型在处理作物蒸散发量时采用的是一个阶段的平均值，所以在第1个时间段先模拟值大于实测值，而后模拟值小于实测值；而第2个时间段为模拟值一直大于实测值。在模型运行第1个时间段，NO 101观测孔潜水水位模拟值与实测值的绝对误差在-27～43 mm之间；NO 501观测孔潜水水位模拟值与实测值的绝对误差在-22 ～45 mm之间；NO 505观测孔潜水水位模拟值与实测值的绝对误差在-27 ～47 mm之间。在模型运行第2个时间段，NO 101观测孔潜水水位模拟值与实测值的绝对误差在-34 ～-77 mm之间；NO 501观测孔潜水水位模拟值与实测值的绝对误差在-27 ～-76 mm之间；NO 505观测孔潜水水位模拟值与实测值的绝对误差在-27 ～-79 mm之间。总体上浅层地下水位模拟值与实测值拟合效果较好。虽然随着模型运行，潜水水位模拟值与实测值的误差有增大趋势，但是对于强降雨后短期内的农田排水过程，系统模型能较好地模拟精准滴灌农田的潜水水位变化过程。

表2 第7 d部分观察孔水位实测值与模拟值 m

表3 第8 d部分观察孔水位实测值与模拟值 m

表4 第9 d部分观察孔水位实测值与模拟值 m

注：表2～4中括号外的值为潜水水位模拟值，括号内的值为潜水水位实测值。

图8 NO 101、NO 501、NO 505观测孔地下水位模拟值与实测值拟合结果图

4 结 论

(1)本文根据精准滴灌农田土壤水分的分布特征，利用Richards方程建立了基于有限差分法的精准滴灌农田排水过程系统模型，模拟了精准滴灌农田排水过程。

(2)潜水水位降低梯度沿着农田靠公路端中间区域至排水沟方向逐渐减小，农田中间区域水位降低速度快，而排水沟渠附近区域水位降低慢。灌溉周期内的第1 d水位最高，并且随着时间的推移，水位逐渐降低且水位降低速率逐渐减小，即灌溉后随时间推移农田各点水位逐渐趋于稳定。

(3)在降雨后农田排水过程的模拟中，潜水水位模拟值与实测值在时间和空间上变化趋势一致，两者差异较小，吻合较好。虽然随着模拟时间变长，潜水水位模拟值与实测值的差异有增大趋势，但总体上，系统模型能较好地真实反映强降雨后短期内精准滴灌农田的排水过程，可为大面积应用滴灌技术的现代灌区排水系统规划设计和改造以及灌排信息化管理提供参考依据。