基于HYDRUS-2D和SHAW模型联合模拟确定膜下滴灌农田秋浇定额

肖 姚，王 珂，赵 强，毛 俊，伍靖伟

（武汉大学水资源与水电工程科学国家重点实验室，武汉430072）

0 引言

在我国西北干旱-半干旱季节性冻融地区，农业生产实践中常采取秋浇、冬灌等非生育期淋洗灌溉制度来保墒压盐，随着近年来膜下滴灌节水灌溉技术在西北地区的广泛推行，如何对该地区膜下滴灌农田进行非生育期淋洗成为亟待解决的新问题。不同于传统地面灌溉，膜下滴灌灌水仅以滴头为圆心局部湿润根区范围内土壤，膜下区域和膜间区域土壤水盐分布极不均衡，具有其特殊的水盐运移规律和空间分布状态［1-4］；西北地区广泛存在的季节性冻融现象，又进一步加剧了该地区膜下滴灌农田非生育期水盐运动的复杂性［5-8］。因此，针对传统地面灌溉农田的非生育期淋洗制度难以直接应用至季节性冻融膜下滴灌农田。

尽管学者们针对非生育期淋洗制度已展开诸多研究［9-12］，但研究对象多为传统灌溉农田，并未考虑到膜下滴灌农田特殊的水盐分布情况，难以指导膜下滴灌区的非生育期淋洗。目前仅有部分学者针对膜下滴灌区非生育期淋洗进行了研究。孙贯芳等［14］通过田间试验研究了河套灌区玉米膜下滴灌农田土壤水热盐效应及秋浇洗盐灌溉效果，结果表明不同滴灌制度下土壤水盐剖面分布极不均匀，非生育期洗盐灌溉效果显著。毛威［15］等利用SaltMod 模型模拟河套灌区膜下滴灌土壤盐分演化规律，分析了根系层土壤盐分累积与灌溉水量、水质、地下水位等因素的关系，并建议每两年引黄河水秋浇一次。虎胆·吐马尔白等［16］则利用HYDRUS-3D 模型来模拟不同秋浇定额下秋浇期10月1日-11月1日膜下滴灌农田的水盐运移，通过分析秋浇前后的土壤脱盐率和底层土壤水分通量，给出了新疆生产建设兵团石河子121团的推荐秋浇定额。现有研究从土壤水盐运移规律、淋洗灌溉效果、推荐淋洗定额等方面展开了部分讨论，但现有研究时段局限在非生育期淋洗期，未考虑冻融作用对膜下滴灌农田水盐重分布的影响，不能有效确定作物播前土壤水盐状态；且有关淋洗定额的研究均为一年一淋洗条件下的定额优选，缺乏膜下滴灌农田年际间水盐运移规律的有关研究，对于膜下滴灌农田多年一淋洗的淋洗制度研究尚未进一步展开。

针对上述问题，本文综合考虑膜下滴灌农田水盐空间分布的特殊性和季节性冻融影响，采用多维水盐运移模拟模型HYDRUS-2D 和冻融模拟模型SHAW 进行联合模拟，同时在进行单年际水盐运移模拟基础上，进一步模拟了膜下滴灌农田多年际水盐运移情况，确定了季节性冻融区膜下滴灌土壤单年际和多年际的推荐秋浇定额，以期为季节性冻融区膜下滴灌农田确定合理的非生育期淋洗制度提供参考。

1 材料与方法

1.1 试验区概况

田间试验于2017-2018年在河套灌区隆胜井渠结合试验区膜下滴灌玉米大田内进行，试验区位于河套灌区中部永济灌域，地理坐标为107°28′～107°32′E，40°51′～40°55′N，海拔1 037 m。试验区位于干旱-半干旱气候带，属典型的温带大陆性气候、季节性冻土地区，土壤冻结期始于11月中下旬，在翌年4 中下旬至5月上旬冻层融通。试验区土壤质地主要为壤土、粉壤土，表层土壤容重1.45～1.55 g/cm3，黄河水矿化度0.4 g/L。

1.2 试验设计与数据观测

本文基于野外秋浇试验进行监测，试验田内滴灌带铺设间距为100 cm，地膜宽60 cm，膜间距离40 cm。试验设置1 350、1 800、2 250 m3/hm2三组秋浇定额，于秋浇前（2017年10月22日）、秋浇后（2017年11月10日）及春播前（2018年4月15日）各进行一次土壤取样，取样时在膜下（A 组）、膜间（B 组）位置用荷兰钻取土，取样深度为100 cm，由表层开始分层取样，取样分层为0～5、5～10、10～20、20～30、30～40、40～60、60～80、80～100 cm。膜下、膜间取样点位置及取样剖面形式见图1。取样后测定土壤颗粒级配、容重、含水率、含盐量等基本参数。

图1 取样剖面示意图Fig.1 Schematic diagram of drip irrigation under mulch

同时，对土壤温度、地下水埋深和基本气象数据进行监测。土壤温度由试验区内地下温度传感器自动采集。地下水埋深数据由试验区地下水位监测井测定，每5 天人工测量一次并记录。气象资料采用站内自动气象站观测，每1 h 自动记录一组数据，基本气象资料包括空气相对湿度、风速、风向、最高气温、最低气温、降雨量、蒸发及日照辐射等。

2 模型耦合方法及其率定与验证

2.1 HYDRUS-2D和SHAW模型耦合方法

本文结合HYDRUS-2D 模型和SHAW 模型进行数值模拟。HYDRUS-2D模型是用于模拟饱和-非饱和多孔介质中水分、溶质及能量运移的有限元计算模型，软件可以描述复杂的边界形状，设置多样化的边界条件，目前在对滴灌点源入渗条件下土壤水、热及溶质运移的数值模拟方面应用广泛。SHAW 模型对系统各层结构之间水、热、溶质运动的物理过程有明确的数学公式表述，常用于模拟在土壤冻结和融化过程中的水量、热量和溶质通量变化。

考虑到膜下滴灌条件下土壤水盐空间分布的差异性以及季节性冻土区冻融作用的特殊性，本文将全年划分为生育期、秋浇期、冻融期3 个阶段，联合使用HYDRUS-2D 模型和SHAW模型模拟农田水盐运移过程，其中HYDRUS-2D 模型模拟生育期和秋浇期，SHAW 模型模拟冻融期。为将二维的HYDRUS-2D 模型和一维的SHAW 模型耦合使用，根据最不利原则，采用HYDRUS-2D模型秋浇期末输出的膜间积盐最大剖面的一维数据作为SHAW模型模拟冻融期的初始值。

2.2 模型定解条件

HYDRUS-2D模型中取两滴灌带之间的一半土体作为模拟区域，模拟区域为一个水平宽50 cm，纵向深100 cm 的矩形研究区。模型初始条件为隆胜试验区实测值。模型区分生育期和非生育期来设定上边界，生育期模型滴头处设为变流量边界，膜内设为零通量边界，膜间设为大气边界；非生育期无地膜覆盖，设为大气边界。模型下边界选择为自由排水边界，因为试验区地下水位长期低于2.5 m，土壤水和地下水之间联系较微弱。由于模型概化的模拟区域为平面对称，故侧边界选择为零通量边界。

SHAW 模型中设置各层节点深度分别为0、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1.0 m。模型初始及边界条件由相应气象文件、土壤温度文件、节点含水率文件和模型参数文件输入。气象参数、土壤温度参数为隆胜试验区实测值，初始土壤含水率为HYDRUS-2D模型输出值。模型上边界条件由隆胜试验区实测气象观测数据确定，下边界条件由土壤底层体积含水量资料和土壤底层温度资料确定。

2.3 模型率定和验证

利用2017-2018年土壤含水率和土壤含盐量的野外实测数据进行参数率定和验证，用A组数据率定模型参数，并用B组数据对初步率定的结果进行验证。模拟值和实测值的拟合程度采用均方根误差RMSE进行评价，计算公式如下：

式中：Si为第i个样本的模拟值；Oi为第i个样本的实测值；n为观测样本数目。

HYDRUS-2D 模型中，不同土层土壤含水率的RMSE为0.021 5～0.031 3 cm3/cm3，秋浇后大田各取样点不同位置剖面含水率模拟值与实测值基本一致，土壤含水率模拟值和实测值吻合较好。不同土层土壤含盐量的RMSE为0.268～0.339 g/kg。含盐量的模拟效果比含水率的模拟效果略差，总体吻合较好。秋浇后土壤含水率和含盐量实测值与模拟值对比图见图2。

图2 秋浇后土壤含水率和含盐量实测值与模拟值对比图Fig.2 Comparison of measured and simulated values of soil moisture and salinity after autumn irrigation

SHAW 模型中，不同土层土壤含水率RMSE为0.021 2～0.031 1 cm3/cm3，春播前土壤含水率模拟值与实测值吻合较好。不同土层土壤含盐量RMSE为0.629～3.77 g/kg，春播前土壤含盐量模拟值与实测值总体吻合较好。春播前土壤含水率和含盐量模拟值与实测值对比见图3。

图3 春播前土壤含水率和含盐量实测值与模拟值对比图Fig.3 Comparison of measured and simulated values of soil moisture and salinity before spring planting

结果表明，HYDRUS-2D 和SHAW 模型分别对秋浇期和冻融期的土壤水盐模拟结果良好，基本能够正确的反应秋浇及冻融作用下的土壤水盐运移规律。故结合HYDRUS-2D和SHAW模型来制定膜下滴灌条件下的秋浇制度可行。

3 膜下滴灌区秋浇定额确定

3.1 单年际秋浇定额

本文采用HYDRUS-2D 和SHAW 模型模拟膜下滴灌农田引黄灌溉时秋浇期（10月1日-11月10日）和冻融期（11月11日-翌年4月15日）的土壤水盐运移过程，输出春播前的土壤含水率和土壤含盐量，与相应土壤水盐标准对比后进一步确定单年际秋浇定额。秋浇应保证春播前作物根系活动层内土壤的水盐条件满足作物生长的需求，因此以满足春小麦苗期正常生长发育的土壤水盐条件为指标来确定单年际秋浇定额。春播前土壤0～40 cm 体积含水率应为25%～28%，土壤表层0～10 cm土壤含盐量应小于2 g/kg。

为模拟不同秋浇定额对土壤水盐运移的影响，模拟过程中针对不同盐渍化程度的土壤分别设定0、450、900、1 350、1 800、2 250、2 700、3 150 及3 600 m3/hm29个等级灌溉量的灌溉处理，根据秋浇定额对春播前土壤含水率和含盐量的影响关系来确定适宜的秋浇定额。土壤盐渍化程度的划分采用童文杰等［17］结合河套灌区土壤盐分普查和典型作物耐盐性评价后给出的分类方法，初始含盐量为试验站在相应盐渍化区间内的实测值的平均值。图4～6 分别为轻、中、重度盐渍化土壤条件下，春播前0～40 cm 土壤含水率及0～10 cm 土壤含盐量随秋浇定额的变化关系图。从图中可以看出，秋浇具有良好的增水保墒、盐分淋洗效果，且秋浇定额越大，春播前土壤0～40 cm 含水率越高，土壤0～10 cm 含盐量越低，但随着秋浇定额的增加，相应的边际效益逐渐减小。

图4 轻盐土春播前0～40 cm土壤含水率和0～10 cm土壤含盐量随秋浇定额变化图Fig.4 Changes of 0～40 cm soil water content and 0～10 cm soil salt content with autumn irrigation quota before spring sowing in light saline soil

模拟结果表明：轻度盐渍化土壤对应于适宜含水率和含盐量范围的秋浇定额为1 100～2 300 m3/hm2，出于节水考虑，推荐秋浇定额为1 100 m3/hm2。中度盐渍化土壤对应于适宜含水率和含盐量范围的秋浇定额为2 100～2 300 m3/hm2，出于节水考虑，推荐秋浇定额为2 100 m3/hm2。重度盐渍化土壤对应于适宜体积含水率范围的秋浇定额为1 100～2 300 m3/hm2，而对应于春小麦耐盐极限的适宜秋浇定额应大于2 500 m3/hm2。在重度盐渍化土壤条件下，当使用较大灌水定额以满足盐分淋洗目的时，会使土壤含水率过大而造成潮塌，当使用较低灌水定额以满足适宜含水率要求时，却并不能将土壤盐分淋洗至春小麦耐盐极限以下，秋浇很难同时满足春播前土壤0～40 cm 含水率要求及0～10 cm 含盐量要求。在实际生产中，对于重度盐渍化土壤一般不种植小麦，而以种植葵花等较为耐盐的作物为主。在重度盐渍化土壤条件下，为淋洗盐分、改良土质，建议秋浇定额为2 500 m3/hm2。

图5 中盐土春播前0～40 cm土壤含水率和0～10 cm土壤含盐量随秋浇定额变化图Fig.5 Changes of 0～40 cm soil water content and 0～10 cm soil salt content with autumn irrigation quota before spring planting in medium saline soil

3.2 多年际秋浇定额

为得到更为节水的秋浇淋洗方案，本文进一步研究了黄灌时多年际秋浇淋洗制度。通过运用HYDRUS-2D 和SHAW 模型进行连续多年的土壤水盐运移模拟，得到多年不秋浇情况下的土壤盐分累积规律，并据此进一步确定秋浇频次和不同秋浇频次下对应的秋浇定额。

3.2.1 确定秋浇频次

首先根据作物苗期含盐量需求来确定最长秋浇年限，即苗期0～20 cm 土层土壤含盐量需低于2 g/kg。进行连续多年不秋浇数值模拟后，输出播种前的土壤含盐量数据，当土壤含盐量大于作物生长土壤盐分控制上限时，则必须进行秋浇。模拟结果显示，一年不秋浇情况下轻、中、重度盐渍化土壤翌年的表层土壤含盐量即超过了作物苗期耐盐极限，故研究多年际秋浇定额时，只考虑非盐土情况，初始含盐量为试验站实测值。

图6 重盐土春播前0～40 cm土壤含水率和0～10 cm土壤含盐量随秋浇定额变化图Fig.6 Changes of 0～40 cm soil water content and 0～10 cm soil salt content with autumn irrigation quota before spring planting in heavy saline soil

对非盐土进行连续多年不秋浇数值模拟，模拟结果见图7。模拟结果表明，当不进行秋浇淋洗时，土壤连年积盐，当连续三年不秋浇时，翌年生育期初土壤含盐量即超过了2 g/kg，无法满足作物生长要求，故最长秋浇年限为3年，多年际秋浇频次可设为3年1灌和2年1灌。

图7 非盐土多年不秋浇翌年生育期初0～20 cm土壤含盐量Fig.7 The salt content of the 0～20 cm soil at the beginning of the growth period of the non-saline soil when irrigated in autumn for many years

3.2.2 确定秋浇定额

针对3年1 灌和2年1 灌两种秋浇频次，分别设置多组秋浇定额进行数值模拟，以确定不同秋浇频次下的秋浇定额。为使制定的多年一灌的秋浇制度能应对长期膜下滴灌要求，即能长期连续实施三年一灌或两年一灌的秋浇制度，以秋浇之后翌年4月15日0～20 cm 土层的土壤平均含盐量小于或接近于模拟初始时设定的0～20 cm 土层的土壤平均含盐量为标准（即小于或接近0.771 g/kg），选定秋浇定额。

为确定三年一灌条件下的秋浇定额，设定1 200、1 400、1 600、1 800、2 000、2 200、2 400、2 600、2 800、3 000 m3/hm2十组秋浇定额进行数值模拟，第一年和第二年不秋浇，第三年进行秋浇，秋浇后得到的第四年4月15日0～20 cm 土壤平均含盐量见图8。由图8 可知，随着秋浇定额的增加，0～20 cm 土壤含盐量逐渐降低，且含盐量曲线下降趋势逐渐放缓，这说明秋浇的淋盐效率随秋浇定额的增加而减少。当秋浇定额为2 000 m3/hm2时，0～20 cm 土壤平均含盐量开始小于0.771 g/kg，出于节水考虑，选择2 000 m3/hm2为三年一灌时的秋浇定额。

图8 三年一灌0～20 cm土壤含盐量随秋浇定额变化图Fig.8 Map of 0～20 cm soil salt content change with autumn irrigation quota under irrigation every three years

为确定两年一灌条件下的秋浇定额，设定800、1 000、1 200、1 400、1 600、1 800、2 000、2 200、2 400、2 600 m3/hm2十组秋浇定额进行数值模拟，第一年不秋浇，第二年进行秋浇，秋浇后得到的第三年4月15日0～20 cm 土壤平均含盐量见图9。由图9可知，随着秋浇定额的增加，0～20 cm土壤含盐量逐渐降低。当秋浇定额增加到2 200 m3/hm2时，含盐量曲线急剧变缓，这说明秋浇定额大于2 200 m3/hm2后秋浇的淋盐效率会显著减少。当秋浇定额为1 350 m3/hm2时，0～20 cm 土壤平均含盐量开始小于0.771 g/kg，出于节水考虑，选择1 350 m3/hm2为两年一灌时的秋浇定额。

图9 两年一灌0～20 cm土壤含盐量随秋浇定额变化图Fig.9 Map of 0～20 cm soil salt content change with autumn irrigation quota under irrigation every two years

4 结论

本文结合HYDRUS-2D 模型和SHAW 模型，构建了膜下滴灌农田土壤水盐运移数值模拟模型，进行了不同秋浇定额下的土壤水盐动态模拟，定量分析得到了膜下滴灌农田单年际和多年际秋浇定额。主要结论如下：

（1）秋浇定额越大，秋浇的增水保墒效果越好，盐分淋洗效果越明显，但秋浇增水保墒和淋盐的边际效率降低。

（2）单年际秋浇情况下，模拟了膜下滴灌农田黄灌条件下轻度、中度、重度盐渍化土壤不同秋浇定额下的土壤水盐运移，结果表明：轻度盐渍化土壤，推荐秋浇定额为1 100 m3/hm2；中度盐渍化土壤，推荐秋浇定额为2 100 m3/hm2；重度盐渍化土壤，秋浇难以同时满足春播条件下土壤含水量和含盐量要求，推荐种植耐盐作物，在以淋洗盐分、改良土质为目标时推荐秋浇定额为2 500 m3/hm2。

（3）多年际秋浇情况下，模拟了膜下滴灌农田黄灌条件下连续多年的土壤水盐运移过程，结果表明：非盐土至少每3年需进行一次秋浇，3年1 灌时推荐秋浇定额为2 000 m3/hm2，2年1灌时推荐秋浇定额为1 350 m3/hm2；轻度、中度和重度盐碱土必须每年进行秋浇，否则作物无法正常生长。□