考虑可动-不动水体假设的河套灌区适宜秋灌定额UBMOD模拟

2022-04-19伍靖伟杨金忠

农业工程学报 2022年2期

毛威，朱焱，伍靖伟，杨金忠

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

0 引言

西北农业灌区由于干旱的气候条件与不合理的灌排措施而面临土壤盐渍化的威胁。为了控制土壤盐分，灌区每年在作物收获后采用较大的灌溉定额进行灌溉压盐，以便淋洗累积在作物根系层的土壤盐分。以内蒙古河套灌区为例，现状情况下河套灌区年均引黄水量约为47亿m，其中约1/3的年引黄水量被用作秋灌淋盐。随着国家水资源“三条红线”管理制度的实施，河套灌区未来的水资源供需矛盾将更加突出。秋灌作为现状情况下水资源的重要消耗项，必将大幅削减。因此，探究以河套灌区为代表的西北干旱-半干旱灌区作物根系层的土壤盐分积聚规律及淋洗特征，对灌区节水与灌溉农业可持续发展具有重要的现实意义。

目前对农业灌区土壤水盐运动规律已有大量研究，主要方法有田间试验研究、数值模型计算等。现阶段土壤盐分模拟与计算中广泛采用的是对流弥散方程。该方程物理概念清晰，但该方程是在均质情况下推导得到，无法解释溶质穿透曲线在非均质情况下的早期穿透和拖尾等物理非平衡现象。在实际的淋洗试验中，大量水分淋洗之后的土壤仍然保持着一定量的土壤盐分，而对流弥散方程难以描述此类现象。从多孔介质空间结构分析，土壤孔隙可以分为有效孔隙与死端孔隙，有效孔隙连通性较好，水盐运移通畅，而死端孔隙中的水体及盐分难以运动。因此在研究土壤盐分淋洗的问题时，将土壤视为均质介质的对流弥散方程容易高估盐分淋洗效率。不同学者提出了众多分区类模型，以描述溶质运移的物理非平衡现象。有学者在总结了大量模型的实际应用后指出，分区类模型对溶质运移实际问题的应用效果好于对流弥散方程。两区模型将土壤水分划分为可动水体与不可动水体，以可动水体描述土壤有效孔隙中的水分及溶质迁移转化运动，以不动水体描述土壤死端孔隙中的溶质转化过程，因此可以有效地描述非均质介质中的溶质运移转化问题。相关试验及研究表明，自然状态下土壤水体中0～80%的水分属于不动水体，且不动水体对溶质运移的影响显著。然而在现阶段水盐运移实际问题的研究中，较少有研究考虑不动水体的存在及其对水盐运移过程的影响。

本研究主要针对内蒙古河套灌区土壤水盐运移的实际问题。在已开展研究中，笔者基于水盐运移的均衡原理与数值求解方法，考虑可动-不动水体假设，开发了土壤水盐运移模型UBMOD。在对河套灌区土壤水盐运移实际问题的分析中发现，考虑可动-不动水体假设时的模拟效果比不考虑可动-不动水体假设的对流弥散方程更好。本研究基于已开发的数值模型，分析考虑可动-不动水体假设时河套灌区永联试验区根系层盐分累积与淋洗规律，并探究研究区现状及未来可能的根系层盐分累积与淋洗规律及适宜秋灌定额。

1 材料与方法

1.1 研究区概况与数据说明

永联试验区位于河套灌区义长灌域内，属于内蒙古自治区巴彦淖尔市五原县永利公社辖区（107°59′09″～108°01′51″E，41°00′50″～41°08′06″N），如图1a所示。试验区南北狭长约为12 km，东西宽约为3 km，控制面积29.75 km。试验区地势南高北低，海拔从西南面的1 028.9 m逐渐降至东北面的1 025.4 m。试验区南部为皂火干渠，区域内灌溉引水量均由皂火干渠取水引入。试验区其余边界则为排水沟。根据区域内土地利用情况，可以划分为3种土地利用类型，分别为耕地、村庄和荒地，其空间位置分布如图1a所示。试验区地质构造以冲湖积层和冲积层为主，是由湖沼河流相向黄河泛流相过渡的沉积物，主要岩性为黄色亚黏土与粉细砂互层，土质变化复杂。

图1 永联试验区示意图与监测点垂向土质分布 Fig.1 Location of the Yonglian test area and the vertical distributions of soil texture at different monitoring points

永联试验区属于北温带干旱地区，年均降雨量173 mm，主要集中在7—9月，由20 cm蒸发皿实测得到的年均水面蒸发量为2 088 mm，约为降雨量的12倍，因而该地区的灌溉农业非常依赖于地表引水。试验区冬季寒冷干旱，因而只在每年的5—9月种植一季作物，在10月进行秋灌，其余月份则休耕。区内主要种植作物有小麦、向日葵、玉米等。

试验区在2008年进行了区域土壤水盐运移的监测试验，共有13个土壤水分与盐分监测点，其分布如图1a所示。其中第1、2、3、5、6和13号监测点位于耕地，4号和8号监测点位于村庄，第7、9、10、11和12号监测点位于荒地。各监测点均进行了详细的土壤粒径分析，依据国际土壤质地分类标注的土壤质地垂向分布情况见图1b所示，各监测点土质展现出强烈的非均质性。各土质所采用的土壤水力参数如表1所示，不同土质的水力参数的取值是基于美国水土保持局土壤参数数据库获得，田间持水率（θ）基于实测数据获得。在该年的5—11月，每月对监测点采样3次进行土壤水分的观测。对于盐分而言，一共有3次监测，分别是在4月26日，7月26日和10月6日，分别对应作物种植前、作物生育期和作物收获后的土壤盐分状况。土壤水盐监测点在垂直方向上近地表处取样较密集，取样深度分别为0～20、>20～40、>40～60、>60～80、>80～120和>120～160 cm，所有监测结果均为2次重复的平均值。2019年，监测点13于9月和11月开展了秋浇期间土壤水盐情况监测。永联试验区的气象数据由附近的五原气象站获得，灌溉数据由灌区管理机构提供。

表1 各土质土壤水力参数 Table 1 Soil hydraulic parameters with different soil textures

1.2 UBMOD水盐均衡模型

UBMOD模型基于传统水盐均衡模型，结合统计方法与数值方法构造得到。相比于求解Richards方程的数值模型，UBMOD模型具有更好的计算效率与模型稳定性，对土壤水分运动过程的描述仅需4个稳定易获取的土壤物理参数，且对时空计算步长没有严格的要求，适合于模拟农田水盐运移的实际问题。UBMOD代码网址为https://github.com/Weiwei-Mao/UBMOD_S。

在水分运动模块中，UBMOD模型综合考虑了影响土壤水分运动的重力势、源汇项等外力和基质势的影响。将研究区域进行垂向离散，在一个计算时间步中，可以将水分运动划分为4步，即入渗水的分配、重力势作用下的土壤水分运动、蒸散发等源汇作用下的土壤水分运动和基质势作用下的土壤水分运动。当存在水分入渗时，入渗水从上到下依次按照饱和含水率对各分层进行填充，对特定土层，其填充水量为

式中为单位面积的该层土壤所分配的水量，m；是该层土壤的厚度，m；是该土壤的含水率，m/m；是该土壤的饱和含水率，m/m；是单位面积的入渗水的总量，m；是在入渗过程中已经分配过的入渗水量，m。

非饱和带土壤水分的运动由以下控制方程描述：

基于可动-不可动水体假设，构建溶质运移模块。在可动水体部分，考虑溶质运移的对流运动、化学反应等源汇项、一阶质量交换运动和扩散运动，在不动水体部分仅考虑溶质的化学反应等源汇项和一阶质量交换运动。

可动水体部分存在水分的运移，溶质随之发生对流运动，其控制方程为

式中是土壤中可动水体含量，m/m；为可动水体的浓度，kg/m；为可动水体的平均孔隙流速，m/s，其与土体中水分通量的关系为

可由水分模型求得，m/s。

溶质运移方程中包含多种源汇项与物理化学反应的影响，可动水体与不动水体的控制方程分别为

式中是土壤中不可动水体含量，m/m，=-；为不动水体含水量占总含水量的比例，=/；为不可动水体的浓度，kg/m；和分别为可动水体部分和不可动水体部分土壤的干密度，kg/m；和分别为可动水体部分和不可动水体部分吸附在土壤颗粒上的固体溶质浓度；和分别为可动水体部分和不可动水体部分液体的一阶反应速率系数，s；和分别为可动水体部分和不可动水体部分固体的一阶反应速率系数，s；和分别为可动水体部分和不可动水体部分液体的零阶反应速率常数，kg/(m·s)；和分别为可动水体部分和不可动水体部分固体的零阶反应速率常数，s；为结晶态盐分溶解速率，s；为该盐溶解平衡时的浓度（或溶解度），kg/m。零阶和一阶速率常数表示一系列的反应与转化过程，包括生物降解、挥发、析出与放射性衰减等。

水平向不动水体与可动水体之间溶质的交换作用，不可动区域的溶质仅以扩散形式与可动区域发生质量交换，且假设质量交换速率与两区的溶质浓度成正比，其控制方程如下：

式中为质量交换系数，s。

之后，计算可动水体的垂向扩散运动，其控制方程为，

式中为可动水体的弥散系数张量，在垂向一维情况下仅考虑其垂向分量，m/s。

1.3 模型率定与验证

对每种土地利用类型，各选择两个监测点的数据进行参数率定，而采用剩余监测点的数据进行验证。需要率定的参数有土壤盐分运移的纵向弥散度（）、不动水体含量（）和一阶质量交换系数（）。

采用相对误差（Mean Relative Error，MRE）和均方根误差（Root Mean Square Error，RMSE）判断模型拟合结果，其计算公式如下：

式中Y为参照值或实测值，y为模型模拟值，为样本数量。采用盐分存储指数（Soil Storage Index，SSI）量化根系层盐分的变化。盐分存储指数是盐分储存变化量与初始存储量的比值，其计算公式为

式中和分别为根系层在研究初期与末期所含盐分量，kg。SSI为正值表示土壤盐分增加，根系层处于积盐状态；SSI为负值表示土壤盐分减小，根系层处于脱盐状态。

1.4 情景设置说明

以0～50 cm土层作为根系层，根据文献资料，河套灌区根系层土壤盐分小于2 g/kg时为非盐碱化土壤，基本满足绝大部分作物生长发育的需求；当根系层土壤盐分在2～3 g/kg之间时，为轻度盐渍化土壤，会影响作物生长发育；当根系层土壤盐分大于3 g/kg时，对大部分作物生长发育均有严重的影响。因此，本研究将根系层盐分的淋洗标准设置是：秋灌过后根系层土壤盐分控制在2 g/kg以下。设置不同水文年型、不同生育期灌溉定额、不同灌溉水矿化度，探讨其对根系层秋灌淋盐的影响。

不同水文年情景：2008年降雨量较大，为典型的丰水年。采用Pearson-Ⅲ型曲线对五原气象站58a的实测降雨资料进行排频分析，并取25%、50%和75%保证率的年型作为丰水年、平水年和枯水年的代表，其各自对应的降雨量分别为209、162和125 mm。研究区在作物生育期的灌溉量采用2000—2013年的均值，因而取生育期灌溉量238 mm。秋灌定额设为60、80、100、120、140 mm，结合3种水文年型，共15种情景。

不同生育期灌溉定额情景：选择水文年型为平水年，生育期灌溉定额设为178、208、238、268、298 mm，并秋灌定额设为60、80、100、120、140 mm，结合5种生育期灌溉定额设置，共25种情景。

对于灌溉水矿化度，以上情景中均采用引黄水进行灌溉，灌溉水矿化度较小（0.65 g/L）。河套灌区地下水矿化度多年平均值为4.08 g/L，其中43%～49%区域的地下水矿化度小于3.0 g/L。矿化度小于2 g/L的地下水水源满足公认的农业灌溉用水水质标准，而只要满足作物对水分的需求并控制盐分危害，利用微咸水（2～3 g/L）进行农业灌溉也能获得高产。因此在未来引黄水量减小的背景下，有可能混合矿化度较高的地下水与引黄水进行灌溉。设置不同的灌溉水矿化度浓度，分别为1.0 、1.5、2.0、2.5和3.0 g/L。选择水文年型为平水年，生育期灌溉定额为238 mm，秋灌期灌溉定额同样设置5种不同的情况，共25种情景进行模拟计算。

2 结果与分析

2.1 模型率定和验证结果

各观测点在时间上平均的率定与验证的水分及盐分MRE和RMSE见表2所示，盐分剖面在2008年10月6日的对比如图2所示。

表2 监测点率定及验证结果 Table 2 Calibration and validation results at monitoring points

图2 10月6日土壤盐分模拟值与实测值 Fig.2 Observed and simulated soil salt content on October 6

由图2可知，模型模拟结果与实测结果较为接近。在率定情况下，耕地、村庄和荒地的剖面平均含水率MRE均小于20%，RMSE值小于0.06 m/m，剖面平均全盐量的MRE小于30%，RMSE小于0.23 g/kg。基于剖面的模拟结果，率定得不动水体含量为0.2，纵向弥散度为0.4 m，一阶质量交换系数为0.001 d。在验证情况下，耕地、村庄和荒地的剖面平均含水率MRE均小于30%，RMSE值小于0.07 m/m，剖面平均全盐量的MRE小于25%，RMSE小于0.17 g/kg。可见，模型可以模拟试验区土壤水盐运移的实际情况。需要说明的是，图2显示的盐分剖面在垂向上存在锯齿状的突变，主要由垂向复杂的土质分层所导致。此外，土壤水盐的部分模拟结果与实测值有一定差距，如8号点的表层结果，12号点的深层结果等。产生该结果的原因有，土壤具有强烈的非均质性与空间变异性，使得模型模拟难以刻画所有水盐运动的细节。此外，位于村庄、荒地的监测点（4号和8号监测点位于村庄，第7、9、10、11和12号监测点位于荒地）处于道路、房屋附近，受人为活动的干扰较大。

此外，采用2019年第13号监测点秋灌期间的实测土壤水盐结果进一步验证模型的模拟效果。以9月15日采样结果作为初始值，11月18日土壤含水率与全盐量的模拟值与实测值对比情况如图3所示。可以发现，模型模拟的土壤含水率与实测结果较为接近，剖面含水率MRE为9.4%，RMSE值为0.10 m/m，剖面全盐量的MRE为11.8%，RMSE为0.13 g/kg。综上，模型能可靠地模拟试验区土壤水盐运移的实际情况。

图3 2019年13号监测点秋灌结束后土壤水盐模拟值与实测值 Fig.3 Observed and simulated soil water content and salt content in the No.13 measuring point after autumn irrigation in 2009

2.2 耕地土壤盐分运移规律分析

基于实测值计算2008年各监测点平均的根系层土壤盐分和含水率，其动态变化和剖面分布如图4。由图4a可知，年内耕地根系层的盐分变化情况大致分为2个阶段：作物生育期和秋灌期。在作物生育期，当不存在灌溉或降雨时，可能由于作物蒸腾与土壤蒸发的作用，根系层盐分缓慢增加。当存在灌溉或降雨时，在入渗水分的淋洗作用下，根系层盐分下降。8月16日出现了极端降雨，单日降雨量高达64 mm，导致了根系层盐分含量的显著下降。根系层总体含盐量下降了0.18 g/kg，SSI为-6.3%。到作物生育期末（9月30日），SSI为5.3%。从土壤盐分剖面分布来看，从初始时刻到秋灌前的作物生育阶段，距地表50和100 cm处的土壤含盐量显著增加。总体而言，根系层盐分在作物生育期稍有增加。10月和11月主要为秋灌淋盐过程。秋灌总的灌溉量为140 mm，然而由于灌溉量比较集中，土壤盐分淋洗较为充分，根系层总体含盐量下降了约0.58 g/kg，SSI为-21.1%。秋灌前后，剖面表层60 cm的土壤盐分显著下降，土壤盐分主要被淋洗至100 cm及以下的深层土壤中，秋灌期洗盐效果明显。

图4 根系层土壤盐分及土壤含水率变化 Fig.4 Variation of soil salt content and water content in root zone

8月16日极端降雨天降雨量为64 mm，雨水的盐分浓度为0，秋灌总灌溉量为140 mm，盐分浓度为0.65 g/L。然而对于根系层而言，8月16日极端降雨情况下，根系层SSI仅为-6.3%，而秋灌期SSI为-21.1%。极端降雨与秋灌前后土壤剖面含水率对比如图4c所示。极端降雨在短时间内有大量水分入渗进入土壤，使得表层土壤的含水率急剧变大。但是降雨量的大部分均留存在0～70 cm的表层土壤中填补水分亏缺，深层入渗水量较小。与之相对应，秋灌水量充足且持续时间较长，使得剖面整体含水率都有显著的增大，可以充分填补非饱和带的水分亏缺，且尚有剩余水量持续下渗，因而可以将大量根系层盐分淋洗至深层土壤甚至地下水中，取得显著的盐分淋洗效果。因而，对于本研究区而言，秋灌是作物根系层盐分淋洗最直接有效的方式。

2.3 模型模拟适宜秋灌定额

不论枯水年、平水年或是丰水年，在作物生长季，由于较大的土壤蒸发和作物蒸腾量，使得根系层土壤盐分在生育期具有逐渐累积的趋势，因此一定需要秋灌淋洗土壤盐分。不同水文年型对根系层淋盐过程影响明显（图5）。不同水文年型不同秋灌定额下的SSI见表3。

图5 不同水文年型根系层土壤盐分淋洗情况 Fig.5 Soil salt leaching results in root zone of different hydrogeological years

由表3可知，当秋灌定额为60～140 mm时，在枯水年不同情景的SSI值为-10.3%～5.7%，在平水年不同情景的SSI值为-12.7%～3.4%，在丰水年不同情景的SSI值为-17.0%～4.0%。可见，秋灌定额相同时，降雨越多，秋灌期土壤盐分淋洗越充分。究其原因，生育期降雨量并未引入更多的盐分，而较大的降雨量会使得土壤盐分累积在根系层较深的位置，秋灌更容易淋洗。总体而言，当水文年型相同时，随着秋灌定额的增大，根系层在秋灌作用下的脱盐程度变大。但在不同的水文年型当秋灌定额相对较小时，随着秋灌定额的增加，也可能会出现SSI上升，根系层土壤盐分淋洗效果变差的情况。如表3中，当秋灌定额由60 mm增加至80 mm时，SSI值从2.6%变为5.7%，即根系层盐分淋洗的效果反而变差。而进一步增大到100 mm及更大的秋灌定额后，根系层盐分淋洗效果会随着淋洗水量的增大而越来越好。可能因为秋灌之前根系层存在较大的水分亏缺，当秋灌定额较小时，秋灌引起的深层渗漏水量较小，而该水量带离根系层的土壤盐分量小于秋灌定额本身所含有的盐分量，因而造成了根系层土壤盐分不降反升。

表3 不同情景下秋灌前后根系层盐分存储指数 Table 3 Soil Storage Index (SSI) values of the root zone before and after the autumn irrigation for different scenarios %

生育期不同灌溉定额与秋灌定额情景下根系层土壤盐分变化和SSI值如图6所示和表3所示。当灌溉定额为298 mm时，秋灌定额从140 mm减小至60 mm时，不同情景的SSI值为-16.8%～3.7%；当灌溉定额为268 mm时，秋灌定额从140 mm减小至60 mm时，不同情景的SSI值为-14.8%～4.3%；当灌溉定额为238 mm时，秋灌定额从140 mm减小至60 mm时，不同情景的SSI值为-12.7%～3.4%；当灌溉定额为208 mm时，秋灌定额从140 mm减小至60 mm时，不同情景的SSI值为-10.9%～5.1%；当灌溉定额为178 mm时，秋灌定额从140 mm减小至60 mm时，不同情景的SSI值为-9.9%～5.8%。可以发现，生育期灌溉定额的影响较为复杂。较大的生育期灌溉水量引入了更多的盐分，因而在生育期末，根系层所含有的盐分量也较大。如图6所示，因趋势一致，238 mm时图省略，当灌溉定额为298 mm时，生育期结束时根系层全盐量为2.28 g/kg，而当灌溉定额为178 mm时，生育期结束时根系层全盐量为2.11 g/kg。因此对于生育期灌溉定额较大的情景，其秋灌所需淋洗的土壤盐分更多。在该算例中，同样出现了当秋灌定额增大，根系层盐分淋洗效果变差的情况，如图 6a所示。此外，对比不同情景可以发现，当秋灌定额大于100 mm时，不同灌溉定额情景下的作物根系层淋盐结果非常接近。可能主要是因为虽然生育期灌溉定额较大时，根系层所含土壤盐分更多，但是生育期较大的灌溉定额也将整个根系层盐分冲洗至根系层下部，因而更利于淋洗。

图6 不同灌溉定额根系层土壤盐分淋洗情况 Fig.6 Soil salt leaching results in root zone of different irrigation quota

不同灌溉水矿化度与秋灌定额条件下根系层土壤盐分变化情况如图7所示。当灌溉水矿化度为1.0 g/L时，秋灌定额从140 mm减小至60 mm时，不同情景的SSI值为-15.5%～6.8%；当灌溉水矿化度为1.5 g/L时，秋灌定额从140 mm减小至60 mm时，不同情景的SSI值为-14.2%～8.0%；当灌溉水矿化度为2.0 g/L时，秋灌定额从140 mm减小至60 mm时，不同情景的SSI值为-13.0%～9.0%；当灌溉水矿化度为2.5 g/L时，秋灌定额从140 mm减小至60 mm时，不同情景的SSI值为-12.0%～10.0%；当灌溉水矿化度为3.0 g/L时，秋灌定额从140 mm减小至60 mm时，不同情景的SSI值为-11.1%～10.8%。相同秋灌定额情况下，灌溉水矿化度与SSI呈现明显的正相关关系，灌溉水矿化度每升高1 g/L时，SSI平均增加0.02。随着灌溉水矿化度的增大，通过灌溉引入的盐分越来越多，秋灌淋洗效果变差。当灌溉水矿化度为1.5 g/L，秋灌定额为120 mm时，秋灌结束后根系层全盐勉强可以达到2.03 g/L。当灌溉水矿化度大于等于2.0 g/L时，秋灌定额为140 mm亦不足以将土壤根系层盐分淋洗至2 g/L之下。

图7 不同灌溉水矿化度根系层土壤盐分淋洗情况 Fig.7 Soil salt leaching results in root zone of different salt concentration of irrigation water

本文研究表明，对于作物生育期降雨较大的丰水年，采用100 mm的秋灌定额即可将根系层盐分淋洗至2 g/kg左右，对于平水年，则需要120 mm的秋灌定额方可将根系层盐分淋洗至2 g/kg，而对于枯水年，则需要140 mm的秋灌定额。可见，不同水文年型对秋灌淋洗作用影响较大。为了保证秋灌对根系层盐分淋洗的效果，在考虑不同水文年型的情况下建议本试验区多年秋灌淋盐的平均水量定额采用120 mm。

在不同生育期灌溉定额情景下，秋灌定额为120 mm时可以将根系层盐分淋洗至2 g/kg以下，满足根系层盐分淋洗的要求。

在不同的灌溉水矿化度情景下，当秋灌定额为120 mm时，需要控制灌溉水矿化度小于1.5 g/L，才可达到根系层盐分淋洗要求。

综上，不同水文年型、灌溉水矿化度对秋灌定额的影响较大，而不同生育期灌溉定额对盐分淋洗所需秋灌定额的影响较小。总体而言，当秋灌定额为120 mm，灌溉水矿化度小于1.5 g/L时，可以满足根系层的淋盐要求。