基于FlexPDE的膜孔灌土壤水分运动数值模拟

赵文刚，范严伟，刘晓群，石 林，邢旭光，马孝义，宋 雯

(1.湖南省水利水电科学研究院，湖南 长沙 410007； 2. 兰州理工大学能源与动力工程学院，甘肃 兰州 730050；3.西北农林科技大学水利与建筑工程学院 旱区农业水土工程教育部重点实验室，陕西 杨凌 712100)

土壤水分运动是陆地水循环的重要环节，也是土壤圈物质循环的重要组成部分，还是土壤中能量和溶质传输的主要载体[1]。掌握土壤水分运动规律对改进灌溉技术、合理利用水资源以及提高作物产量具有重要意义[2]。

目前，土壤水分运动规律研究主要依赖室内试验、数值模拟。室内试验因人工填土、夯实导致土壤各向同性受到影响；同时，对于多因素多水平试验方案设置难以实现完全试验，且完全试验时无法保证因素绝对固定。数值模拟试验依托Richards方程利用有限元法求解，已形成的软件主要有HYDRUS和SWMS，其被广泛应用于灌水方式[3-5]、水质状况[6]、土壤分层情况[7-8]等对土壤水分运动的影响研究中，但缺乏实测资料支撑。基于二者优缺点，通常采用室内试验、数值模拟相结合的方法。

SWMS软件模拟土壤水分运动精度较高，但其计算前期需手动剖分网格、计算单元节点控制长度等工作，工作量大且易出错；此外，其基于Fortran 语言开发导致用户界面不友好。HYDRUS软件虽已实现网格自动划分，但网格划分依赖经验且网格密度单一，对于复杂计算区域，无法局部加密处理，影响计算精度。此外，为简化计算，SWMS和HYDRUS等软件模拟土壤水分运动规律时采用饱和导水率代替非饱和导水率，而研究发现非饱和导水率受孔隙状况、温度、含水率影响[9-10]。

FlexPDE是一款专门求解偏微分方程的软件，具有网格自动划分、网格数量和密度随误差自动调整等优势，能通过编辑问题描述脚本文件将偏微分方程描述的系统转化为有限元模型求解，用户界面友好。Bader[11]利用FlexPDE对粘土薄膜的水分和溶质运移特性进行了模拟，且效果理想；Albrieu[12]采用FlexPDE模拟了定体积土壤中根系生长的吸水模型；Rühaak[13]通过FlexPDE模拟冰川效应下的三维水力机械耦合地下水流动；唐延贵[14]通过压力板仪和非饱和土三轴仪开展试验研究，结合FlexPDE分析了非饱和土-水特征和变形特征；此外，FlexPDE在温度场、电场和热传导[15-17]等方面也有所应用。

相比常规地面灌溉而言，膜孔沟灌在重力影响下的垂向运移与水势梯度作用下的侧向运移共同作用对土壤进行局部湿润的灌溉，使得灌水更加均匀、高效，也更节水。综合上述，本文采用土壤扩散率与土壤比水容量的函数关系推求土壤非饱和导水率；建立土壤储水量模型，分析土壤储水量与累积入渗量的关系，验证FlexPDE软件在模拟膜孔沟灌土壤水分运动规律可靠性的同时也为合理确定灌溉技术要素提供参考依据，以便更好地指导农业生产。

1 材料与方法

1.1 试验设计

供试土壤取自杨凌三道塬、渭河滩0～30 cm土层的粉粘壤土和砂壤土(颗粒组成见表1)。初始含水率取60%田间持水率[18]，膜孔间距一般为12～30 cm、膜孔直径3～8 cm、灌溉水深2～10 cm[19-20]。据此，研究共设置4个处理，具体布设方案列于表2，两种土壤的水力参数由RETC软件拟合获取(表3)。其中处理T1、T2采用A土箱进行试验，处理T3、T4采用B土箱进行试验。

表1 供试土壤颗粒组成

注：土壤质地分类依据国际制进行分类，下同。

Note: The soil texture was classified by international system, the same below.

表2 土箱试验方案

表3 供试土壤水分特征参数拟合

注：ρ、θr、θs、α、n分别代表土壤容重、残余体积含水率、饱和体积含水率、土壤进气值的倒数、参数，下同。

Note:ρ、θr、θs、α、andnrepresent the bulk density, residual volumetric water rate, saturated volumetric water rate, the reverse of air entry value, and parameters, the same below.

1.2 试验方法

膜孔灌溉装置由马氏瓶、有机玻璃土箱和膜孔装置组成(图1)：马氏瓶用于维持恒定水头；土箱规格(长×宽×高)为10 cm × 10 cm × 60 cm(A)和15 cm × 15 cm × 60 cm(B)，土箱A水平向从距边壁2.5 cm间隔2.5 cm、垂向从距上边壁5 cm间隔5 cm设置取土孔，土箱B水平向从距边壁3 cm间隔3 cm、垂向从距上边壁5 cm间隔5 cm设置取土孔，在土箱相邻两面对应位置均设置取土孔(可视为重复)；对于膜孔装置，将1/4膜孔(即土箱一角)作为研究对象。土样经风干、过2 mm筛，按设定容重每5 cm厚度装填，土壤装填完成后，除膜孔处均采用薄膜覆盖以防止水分散失；试验过程中监测马氏瓶水位获得入渗量；实验结束时，从取土孔取土，采用干燥法测定不同深度土壤含水率及土壤剖面水分分布特征。

土壤扩散率采用水平土柱法测定，有机玻璃土柱材质，直径(内径)5 cm、长36 cm。土样经风干、过2 mm筛，按设定容重将含水率为风干含水率土壤均匀装入水平土柱中，试验过程中忽略重力作用，马氏瓶供水，记录马氏瓶水位变化及湿润锋每过1 cm所用时间；试验结束时，从湿润锋处开始取土，采用干燥法测定土壤含水率。

a.膜孔装置；b.取土孔；c.通气孔；d.土面；e.马氏瓶 a. film-hole settings; b. soil collecting hole; c. air vents; d. soil surface; e. Markova bottle图1 膜孔灌实验装置Fig.1 The experimental equipment of film hole

2 模型建立

2.1 土壤水分运动模型

膜孔灌溉条件的Richards方程如下[20]：

(1)

式中，h为土壤基质势(cm)；k(h)为非饱和土壤导水率(cm·min-1)；t为时间(min)；x为水平距离(cm)；z为垂直距离(cm)。

土壤基质势和扩散率与含水率的关系可分别采用van-Genuchten方程[1]和常用的经验公式[21]来拟合，即式(2)、(3)；土壤非饱和导水率与土壤基质势的关系可采用式(4)拟合。

(2)

式中，θr为残余体积含水率(cm3·cm-3)；θs为饱和体积含水率(cm3·cm-3)；n、α、m为经验常数，其中m=1-1/n。

D(h)=D0×exp(-β×(θ0-θ(h)))

(3)

k(h)=C(h)×D(h)

(4)

式中，D0、β、θ为经验常数，θ0为某一特征含水率，本文中取θ0=0；C(h)为土壤比水容量(cm-1)；D(h) 为土壤扩散率(cm2·min-1)。

本研究中，利用FlexPDE软件对土壤水分运动特性进行模拟，且模型可通过式(1)、(4)获得二维条件下的Richards方程：

(5)

2.2 模型初始条件和边界条件

将膜孔直径所在平面作为研究对象(图2)，由于关于z轴对称，故可取OABC围成的区域作为研究对象，s为膜孔半径，L为土箱的长度，H为土箱的高度，则初始条件和边界条件可表述为(T为入渗时间)：

初始条件：

h(x,z,0)=h0(x,z) 0≤x≤L,0≤z≤H

上边界条件：

下边界条件：

左边界条件：

右边界条件：

图2 膜孔灌土壤水分运动简化模型示意图Fig.2 The diagram of simplified model for water flow of film-hole irrigation

2.3 土壤储水量模型

取OABC面作为研究对象(图2)，t和t+Δt时刻的土壤储水量按式(6)、(7)计算，W(t)为t时刻土壤剖面整体含水量，则Δt时间间隔内土壤储水量ΔW可按式(8)计算：

W(t)=∬θ(x,z,t)dxdy

(6)

W(t+Δt)=∬θ(x,z,t+Δt)dxdy

(7)

ΔW=∬θ(x,z,t+Δt)dxdy-∬θ(x,z,t)dxdy

(8)

2.4 FlexPDE软件简介

FlexPDE软件求解偏微分方程问题的可读脚本文件主要由下面几部分组成[22]：TITLE——输入描述标签项；SELECT——控制FlexPDE参数；VARIABLES——偏微分方程变量；DEFINITIONS——定义偏微分方程中出现的参数；EQUATIONS——偏微分方程；BOUNDARIES——边界条件；MONITORS and PLOTS——数值解的图形；END——脚本结束。

3 结果与分析

3.1 土壤水分动态模拟模型有效性验证

土壤剖面含水率分布特征可以反映膜孔灌溉的入渗规律[2]，为了验证膜孔灌溉条件下土壤水分运动模型的有效性，对土壤剖面含水率的实测值与模拟值进行对比分析。

图3为通过FlexPDE软件模拟膜孔沟灌灌水结束时刻的土壤剖面水分分布状况。从图中可以看出：靠近膜孔位置，土壤水分以膜孔为中心呈椭圆形分布；离膜孔位置越远，膜孔对土壤水分分布特征影响越小，呈水平形态分布。

4种处理按照设计的灌水定额所需的灌水时间分别为120、642、209 min和568 min，灌水结束后，通过对比土壤剖面含水率实测与模拟值可以发现，除T3中土壤水分模拟值小于实测值之外，其余处理中土壤水分模拟值均高于实测值，处理T1～T4土壤含水率的平均相对误差分别为9.600%、5.930%、8.340%和14.040%。进一步分析发现，粉粘壤土FlexPDE的模拟效果较砂壤土理想，即对于粉粘壤土，土壤含水率相对误差小于5%、10%和20%的比例分别为40%、68%和84%；砂壤土，其相对误差小于5%、10%和20%的比例分别为12%、40%和72%。而范严伟[23]基于SWMS对此研究得出，各处理平均相对误差分别为10.230%、5.830%、7.660%和8.540%；其中对于粉粘壤土，土壤含水率相对误差小于5%、10%和20%的比例分别为48%、72%和96%，与本研究结果对比发现，FlexPDE的模拟误差较SWMS略高，但二者差异很小，亦满足精度要求，同时鉴于FlexPDE的优点(前已述及)，因此该软件更适合对粉粘壤土的水分入渗进行模拟。

土壤含水率模拟与实测值之间存在一定差异，分析其原因认为：(1)土壤水分特征曲线测定依赖有限点的土壤基质势与含水率拟合，且用土壤脱湿过程来代替土壤吸湿过程，存在一定滞后性；(2)非饱和土壤水分运动参数的测定导致误差；(3)缺少实时监测设备，灌水结束时土壤含水率测定有所延迟；(4)假设土壤各向同性，但实际土壤各向异性；(5)人工装土与设计容重存在差异。

3.2 土壤储水量与累积入渗量的关系

由图4可以看出，各处理累积入渗量、储水量与入渗时间均呈幂函数关系，累积入渗量与土壤储水量的关系可以用线性模型描述，且所有拟合模型均具有较高精度(R2>0.900)。其中，累积入渗量随时间延长而增大，但入渗速率随时间延长而减小；储水量随累积入渗量增大以恒定速率增大。基于储水量与累积入渗量的线性关系，可通过储水量的变化反映累积入渗量的变化情况。

综合图3、4可以看出，基于实测土壤非饱和扩散率和土壤水分特征曲线，FlexPDE对粉粘壤土水分运动模拟效果优于砂壤土，但均满足精度要求。这可能因为砂壤土大孔隙较多，饱和含水率相对较低，相同入渗量使得扩散率变化量增大，土壤湿润锋运移速度加快，含水率分布变化加快，进而导致模拟比较困难。

4 模型应用

膜孔灌，就单个膜孔而言，属于局部不交汇灌溉；但在实际田间管理中，随灌溉水量、历时增加，灌溉水发生交汇入渗，从而导致土壤水分再分布。水分运动及再分布特性受入渗影响显著，而膜孔灌溉水分入渗特性受膜孔直径和间距等因素影响，因此，为提高田间灌溉水利用率，有必要深入研究膜孔直径、间距和灌溉水深等因素对土壤水分运动特性的影响，从而实现田间水的有效管理。

图3 不同处理土壤剖面水分分布Fig.3 The moisture distribution in soil profile under different treatments

图4 储水量和累积入渗量与入渗时间的关系Fig.4 The relationship among soil water storage, accumulative infiltration, and infiltration time

4.1 模拟方案

试验土壤为杨凌一级阶地的粉粘壤土(表4)，利用RETC拟合的水力参数性质见表5。膜孔交汇灌溉水分模拟，重点研究膜孔间距、入渗水深和膜孔直径3个因素对入渗的影响，初始含水率取60%田间持水率，入渗时间设置为300 min，具体模拟试验方案见表6。

4.2 不同处理交汇时间、水分分布和储水量分析

4.2.1 膜孔灌单向交汇时间 由表7可知，对于膜孔间距、直径相同但水深不同的处理，交汇时刻基本相同，这表明入渗水深对膜孔灌单向交汇时间的影响不显著，这可能因为入膜水深虽增大，但只增加了初渗时刻土壤表面膜孔处的水势梯度，起到瞬时加速湿润土壤的作用，一旦膜孔处土壤水分达到饱和时，土壤水势梯度就变为水分饱和土壤与土壤初始土水势之间的差异，入渗水深失去作用。此外，膜孔直径、水深相同条件下，膜孔间距越小，交汇所需时间缩短，这主要因为间距减小，缩短了水分水平运动距离；膜孔间距、入渗水深相同条件下，膜孔直径越大，达到交汇时间越短，这主要由于膜孔孔径增大一方面增大了初始时刻的湿润面，另一方面缩短了水分水平运动距离。

表4 模拟方案供试土壤颗粒组成

表5 土壤水力特性参数拟合

4.2.2 土壤储水量变化 从图5可以看出，储水量随时间延长而增大；相同膜孔直径、膜间距条件下，入膜水深对土壤储水量影响微弱；相同膜孔直径、入膜水深条件下，膜孔间距对于土壤储水量的影响也较微弱。分析其原因认为，储水量变化主要与土壤含水率变化相关，含水率变化与进入土体水量相关，膜孔直径相同的情况下，相同时间周期内水深或膜间距变化使得进入土体水量差异不大，进而导致土壤储水量变化不大。相同膜间距、入膜水深条件下，膜孔直径对储水量影响较显著，且随膜孔直径增大而增大，这主要因为随膜孔直径增大，水分进入土壤通道增大，使得入渗水量增加，进而导致储水量增大。通过对土壤储水量和入渗时间拟合发现，各处理均呈现出良好的幂函数关系(表8)。

4.2.3 土壤水分分布状况 图6显示，在土壤入渗过程中当湿润锋未发生交汇时为自由入渗，水分分布呈椭圆形；当湿润锋发生交汇时，随时间推移土壤水分运动逐渐趋向水平向下推进。对比M1和M2、M3和M4可以看出，膜孔间距较大时，湿润锋发生交汇时间延后，土壤水分分布不均匀；而膜孔间距较小时，土壤水分分布均匀程度高。通过对其它处理对比分析也均反映出与上述类似的特性。由此可知，膜孔间距对土壤水分分布影响显著：膜孔间距越小，土壤水分分布不均匀程度越小；而其他条件相同时，入渗水深越大，相同入渗时间，垂向湿润锋会有所增大，但增大幅度不明显；膜孔间距相同条件下，灌溉水深、膜孔直径对水分入渗的影响比较微弱，这与范严伟[3]、费良军[6]等的研究结果一致。

表6 膜孔交汇灌溉模拟方案

表7 各处理单向交汇时刻

图5 储水量与入渗时间的关系Fig.5 The relationship between soil water storage and infiltration time

处理Treatment拟合公式Fitting functionR2处理Treatment拟合公式Fitting functionR2M1W=0.381t0.7170.999M10W=0.468t0.6940.999M2W=0.339t0.7310.992M11W=0.467t0.6950.999M3W=0.354t0.7300.999M12W=0.457t0.6980.999M4W=0.400t0.7030.999M13W=0.456t0.7180.993M5W= 0.349t0.7310.998M14W=0.453t0.7180.993M6W=0.397t0.7030.999M15W=0.543t0.6860.999M7W=0.461t0.7000.999M16W=0.460t0.7140.992M8W=0.459t0.6980.999M17W=0.541t0.6850.999M9W=0.468t0.6960.999M18W=0.531t0.6880.998

注：鉴于篇幅所限，只列举M1～6，8。 Note: Due to space limitation, only treatments 1 to 6 and 8 are listed.图6 膜孔单向交汇水分分布Fig.6 The moisture distribution of film-hole single-line convergence infiltration

5 结 论

基于对60 cm土层深度的土壤水分分布监测，同时测定了土壤扩散率，并结合FlexPDE软件对膜孔灌土壤水分以及土壤储水量进行模拟，主要得到以下结论：

1)通过Richards方程对土壤水分运动进行模拟可知FlexPDE软件适合于粉粘壤土的土壤水分入渗模拟，满足精度要求；

2)基于建立的土壤储水量模型，储水量与累积入渗量呈现较好的线性关系，且二者均随入渗时间的延长而增加，但增长速率逐渐减小；

3)基于粉粘壤土建立的土壤水分、储水量模型，对膜孔单向交汇入渗进行预测，结果表明，膜孔间距对于交汇初期土壤的水分分布影响较大，膜孔直径对于土壤交汇时间以及储水量影响显著。

本文利用FlexPDE对于膜孔灌溉土壤水分状况进行初步模拟，对于实际膜孔灌溉技术要素的确定具有指导意义。在本研究中，只是对于土壤初始含水率较高的情况进行了模拟，对于土壤初始含水率较低的情况需要进一步研究。