基于HYDRUS的深施作业液肥水分入渗数值模拟

张泽龙，何建炜，张彬，张颖兴，吴子彬，刘海军，曹鹏，薛秀云

（华南农业大学电子工程学院、人工智能学院，广州510642）

0 引言

液肥凭借其施肥方式方便、节肥增产效果显著等特点被广泛应用。液肥中所包含的水分和养分（盐分）是作物生长的最基础条件（李俊峰，2017）。土壤中水分跟养分的实际分布与作物根系最大化地吸收营养具有密切的联系（Wang，2018）。因此充分了解不同条件下液肥的水分迁移规律，对提高作物肥料利用率的积极性，进一步为制定合理的灌溉施肥制度具有重大意义。

数值模拟计算是研究液肥入渗特性的一个重要方法，在数值模拟中HYDRUS模型应用最为广泛（J.Si-munek等人，2012）、Simunek等人（2010）结合室内试验和田间试验，验证了HYDRUS-2D模拟地下滴灌土壤水分分布具有较高可靠性。Bufon等人（2012）利用HYDRUS-2D进行滴灌条件下的棉花含水量模拟，结果表明模拟值与测量值数据拟合程度较好。HYDRUS模型可对不同条件下水盐运移规律进行模拟计算，为区域田间水肥管理提供理论依据和参考价值。

当前我国大部分地区仍采用施肥效率低下、化肥损耗严重的施肥机械或非机械化施肥方式。为此，提出了一种能够定量精准均匀施肥的液肥变量深施控制系统。为给液肥深施系统精准深施效果提供科学依据，进一步提高系统深施效率，对液肥深施系统深施作业条件下的注肥压力、注肥深度、土壤质地3个因素进行HYDRUS建模，探究各因素对水分入渗的影响。

1 试验方法

试验采用单因素试验法进行研究，不考虑各个因素间的交互作用。试验中包括了三个主要影响因素（A系统注肥压力、B系统注肥深度、C土壤质地类型），其中分别对系统压力和注肥深度设置了3个水平，土壤质地设置了2个水平。为更准确地观测土壤湿润体特征值变化规律，系统压力控制水平为：0.03MPa、0.04MPa、0.08MPa。由于成年的南方柑橘树的根系主要集中分布在地表下15～55cm处，故系统注肥深度水平控制为：25cm、35cm、45cm（薛秀云等人，2020）；通过查阅中国科学引文数据库相关文献，了解华南地区主要土壤质地类型，壤质地控制水平为砂壤土、壤土。在数值模拟计算试验中，对以上3因素组合，注肥时间均设定为10min，并在注肥开始后记录垂直剖面下不同时刻9个方向土壤湿润锋运移规律、湿润体形态大小，如图1所示。对模拟试验得到的仿真数据利用Origin-Pro 2018软件进行数据处理及图表绘制，以此来分析不同因素条件下湿润体时空变化规律。

图1 垂直观测面下湿润锋记录示意图

2 数值模拟模型

选择HYDRUS（2D/3D）对液肥深施系统深施作业下土壤水分入渗规律进行动态模拟。HYDRUS模型是由美国盐土改良实验室（US Salinity Laboratory）研发的一款用于模拟饱和-非饱和多孔介质中水分、能量及溶质运移的数值模型。该模型适用于不同初始条件与各种恒定或可变边界条件，能较好地模拟土壤中水分、热以及溶质的时空变化、空间分布及运移规律（迟卉等人，2014）。

2.1 土壤水分运动方程

本研究所探讨的基于液肥变量深施系统的注肥深施入渗为恒定流量条件下的三维柱状面源水分入渗。假定土壤为均质，固相骨架不变形的刚性多孔介质，忽略土壤水分运动中的滞后效应，同时不考虑温度及空气阻力对土壤水分运动的影响，土壤水分运动可用Richards方程来描述：

式中：θ为土壤体积含水率，cm3/cm3；ψ为土壤负压水头，cm；K（θ）为非饱和土壤导水率，cm/min；t为时间，min；x、y为水平面坐标，cm；z为垂直面坐标，cm。

2.2 模拟区域

将试验土壤区域的几何模型规格设置为80cm×80cm×100cm（长×宽×高），注肥枪由直径2cm的圆柱枪体和圆锥形枪头组成，系统注肥处设置为模型计算区域中央，模型入渗面是三维柱状出流边界。采用垂直剖面法观测竖直平面内湿润体形态与特征值，efgh垂直剖面设定为YZ平面，垂直观测面设为距注肥枪口0.1cm处。模型示意图如图2。

图2 模型示意图

2.3 初始条件

本研究模型计算的初始条件为初始条件为初始含水率，因此求解土壤水分运动的初始条件为：

式中：H0（x,y,z）为土壤初始负压，假设计算区域内土壤初始负压处处相等，设置为-100000cm（对应不同土壤质地初始含水率为：砂壤土—0.06cm3/cm3；壤土—0.09cm3/cm3）土壤初始含水量由土壤水分特性曲线转换得到。

2.4 边界条件

如图3所示，本研究模型的上边界A’B’C’D’为土壤表层，为大气边界条件，垂直通量为零；前后左右四面边界AA’B’B、CC’D’D、AA’D’D、BBC’C为自由排水边界。底面边界ABCD，其垂直通量为零。而对于注肥枪末端O及枪管圆柱侧面OE为出流边界，是一个比较复杂的运动边界。试验设置的系统注肥压力较大，注肥枪出口流量远大于土壤入渗速率。因此要将系统注肥压力转化为模型能够设置的边界条件，即为可变压力水头边界处理。对于注肥枪入渗面的可变压力水头按因素组合处理去不同数值还实现。通过系统压力与水柱高度转换公式将系统压力转化为水柱高度：已知一个大气压A0=0.1MPa等于h0=10m高水柱，故利用公式（3）将系统压力转化为压力水头。计算得出压力水头见表1。

表1 三种系统压力对应压力水头

式中：A0为标准大气压，0.1MPa；h0为一个大气压下水柱高度，1000cm；A为系统压力，MPa；h为压力水头，cm。

2.5 数值求解方法

HYDRUS-3D流体建模采用Galerkin有限单元分析法求解三维土壤水分运移问题，时间采用隐式差分法，迭代参数和残差参数均按默认值设置。模型模拟土壤水流区域选用Triangular Prism（四面体）网格进行划分。模型边界离散化采用软件自动生成的边界密度节点。

3 结果与分析

3.1 系统压力对水分入渗的影响

图3为注肥深度为35cm，土壤质地为砂壤土条件下，不同系统压力下湿润锋动态运移过程。由图3可知，系统压力对土壤湿润体形状大小的影响较为明显，湿润锋垂直向上扩散距离、湿润锋水平扩散距离、湿润锋垂直向下扩散距离均随注肥时间的增加而增加，但各方向湿润锋运移速率在注肥后期明显逐渐减小。在整个注肥过程中，由于液肥深施系统施肥为大压力施肥，供水强度大于土壤入渗能力，又因在竖直方向上水分受重力势的影响，导致湿润锋垂直向下运移速率略大于水平运移速率。而在竖直向上方向土壤孔隙相对较大，水分在运移过程中受到阻力相对较小，且在土壤水分入渗过程中土壤基质势起到重要作用，所以垂直向上湿润锋运移速率略大于垂直向下湿润锋。至注肥结束，三个方向的湿润锋运移距离差异很小，所以深施作业下湿润体形状近似“球体”，随着系统压力的增大，湿润体大小也随着增大。由于液肥深施系统系统压力决定了系统输出流量，当系统压力越大时，系统输出流量越大，湿润体形状因而随着增大。

图3 不同系统压力条件下湿润锋运移动态变化图

对在垂直向上、水平扩散距离和垂直向下入渗深度三个方向的湿润锋运移进行拟合，如图4所示。从中可以发现不同系统压力条件下垂直剖面下的垂直向上湿润锋运移距离X、水平湿润锋运移距离R和垂直向下湿润锋运移距离H与注肥时间t均存在显著的幂函数关系，R2均大于0.98，表明拟合结果与模拟值吻合度较高。其中，对于同一湿润锋运移方向上幂函数指数差异性较小，其平均值分别为（a）0.2487、（b）0.33、（c）0.3213。而其幂函数系数均随系统压力的增大而增大，这表明各方向湿润锋运移距离受系统注肥压力影响，均呈现正相关。对于同一系统压力条件下不同方向湿润锋幂函数系数均呈现相同规律，具体表现为：垂直向上湿润锋的幂函数系数>水平湿润锋的幂函数系数>垂直向下湿润锋的幂函数系数。

图4 不同系统压条件力下湿润锋随时间运移距离

3.2 注肥深度对水分入渗的影响

图5为注肥深度为35cm，土壤质地为砂壤土条件下时，不同注肥深度湿润锋动态运移过程。从图5可知，注肥深度对湿润体形态和大小影响较小，不同处理下深施作下土壤湿润体垂直剖面均近似于圆形且形态大小差异很小，湿润锋水平方向最大运移距离略高于注肥枪注肥枪口所在深度。在整个注肥过程中，各处理垂直向上、水平方向、垂直向下湿润锋运移距离接近，至注肥结束，各方向湿润锋运移距离差距始终不大，所以导致了深施模式下土壤湿润体形状最终近似为圆形。虽然不同注肥深度处理对土壤湿润体形态大小影响较小，但对土壤湿润体的空间分布影响很大，随着注肥深度的增加，土壤湿润体分布所处深度越大。因此，针对不同作物种类制定合理的灌溉制度时，要根据作物种类根部生长状况来确定合适的注肥深度，以达到在有效防止地表蒸发和深层渗漏前提下满足作物水肥需求。

图5 不同注肥深度条件下湿润锋运移动态变化图

土壤湿润体的垂直方向（向上和向下）湿润距离和最大水平湿润半径是衡量液肥入渗效果的关键指标。对注肥深度分别为25cm、35cm和45cm三种梯度条件下各方向湿润锋随时间的变化过程进行拟合。从表2中可知，在相同注肥时间内，土壤湿润锋垂直向上、水平和垂直向下运移距离随注肥深度的增加不发生明显变化，说明注肥深度对湿润锋运移影响较小。

3.3 注肥深度对水分入渗的影响

图6为施肥深度为35cm，系统压力为0.03MPa条件下的不同土壤质地湿润锋运移动态变化图。由图6可知，对于不同土壤质地条件下的液肥深施作业而言，土壤湿润体在注肥结束时均近似为“圆形”，湿润体形状差异不大；而湿润锋运移距离却差异较大，各方向湿润锋运移距离在砂质壤土条件下均大于在壤土条件下，由此，导致了土壤湿润锋形状在大小上差异较大。这主要是因为土壤水分入渗主要受土壤孔隙度的影响，土壤孔隙度越多，气相比例越多，土壤内部连通性越好，越有利于土壤水分入渗的进行，刑杰等人（2013）研究表明砂质壤土黏粒含量约为10.43%，壤土黏粒含量约为18.63%，孔隙度随土壤黏粒含量的增加而减小。因此，土壤湿润体体积在砂质壤土条件下大于在壤土条件下。

图6 不同土壤质地条件下湿润锋运移动态变化图

表2 不同注肥深度各方向湿润锋拟合结果

对湿润锋随时间的变化过程进行处理与分析，发现垂直向上湿润锋运移距离X，水平方向湿润锋运移距离R和垂直向下湿润锋运移距离与注肥时间t均呈较好的幂函数关系，即：

式中：a、c、e分别为垂直向下、水平、垂直向上的入渗系数，b、d、f分别为下、水平、垂直向上的入渗指数，数据拟合结果见表3。

表3 各方向湿润锋运移距离与注肥时间拟合结果

从表3中可以看出，无论垂直向下，水平方向还是垂直向上湿润锋的入渗系数均受土壤质地类型影响较大，随土壤黏粒含量的增加而下降。而入渗系数的物理意义表征的是土壤入渗开始后第一单位时间内湿润锋入渗的能力，这也说明了在不同土壤质地在第一单位时间段末湿润锋运移距离差异较大。

4 结语

利用HYDRSUS-3D模型建立了系统深施作业下土壤水分运动模型，对湿润体特征值进行观测，分析了系统压力、注肥深度、土壤质地对土壤水分入渗的影响，得出以下结论：

（1）系统压力对湿润体大小，均随系统压力的增大而增大；注肥深度增大对湿润体形态大小无显著影响，对湿润体空间分布位置有着显著影响，随着注肥深度的下移，土壤湿润体空间分布也随之下移；不同土壤质地条件对湿润体形态影响较小，湿润体形态均近似为圆形，而对对湿润体大小移影响较大，土壤质地越黏重，湿润体越小。

（2）深施作业下土壤湿润体在形状上差异较小，湿润锋运移距离随系统压力的增大而增大，随土壤黏重程度的增大而减小，而与注肥深度无显著关系。不同方向土壤湿润锋运移距离与注肥时间都呈显著幂函数关系，且在相同注肥时间内均符合垂直向上>垂直向下>水平方向的规律。

通过以上研究，可为液肥深施系统的深施作业提供理论依据，但未考虑实际情况下土壤水力特性的差异性，因而仍需通过室内土箱试验进一步研究系统压力、注肥深度、土壤质地对土壤水分入渗的影响。