基于分层水量平衡模型的蓄水坑灌果园土壤水分动态模拟

2023-08-28张甜甜郭向红雷涛孙西欢马娟娟郑利剑

节水灌溉 2023年8期

张甜甜，郭向红，2，雷涛，孙西欢，马娟娟，郑利剑

（1.太原理工大学水利科学与工程学院，太原 030024；2.中国水利水电科学研究院流域水循环模拟与调控国家重点实验室，北京 100038）

0 引言

土壤水分与作物生长息息相关，农田土壤水分的动态变化直接影响作物的生长发育，农田土壤水分含量长期处于过高或过低的水平都会对作物的生长产生不利的影响，因此对田间土壤水分含量进行动态监测与预测对保障农作物的产量具有重要的意义[1]。

对于土壤水分动态的获取方法主要有仪器监测和模型模拟。仪器监测采用专有土壤水分测量仪器长期对田间水分进行测量，比较费时费力，监测成本较高。而近些年随着计算技术的不断发展，模型模拟逐渐成为土壤水分动态研究的重要手段。常见的土壤水分动态模型有三类，即机理模型、神经网络模型、水量平衡模型。机理模型[2-4]是根据土壤水分运动方程并结合实际边界条件，建立不同条件下的土壤水分运动模型，可以模拟土壤水分的空间分布和动态变化，但由于模型为非线性偏微分方程，数值模型的收敛性较差，需要精确水分运动参数，计算较为复杂，比较费机时。神经网络模型[5-7]是目前较为流行的基于数据驱动的预测方法，可以实现对土壤水分的实时动态预测，但该类模型属于黑箱模型，不能解释模型背后的机理，且需要大量的田间数据支撑，才能保证预测精度。水量平衡模型[8]是根据根区土壤水量平衡原理，建立的反映根区水量动态变化的模拟模型，原理简单，所需数据较少，虽然不能精准模拟土壤水分空间分布，但模拟的根区土壤水分动态变化可以满足农田水分管理的需求，已经广泛应用于农田水分动态变化的研究。

目前土壤水量平衡模型结构已由单层水量平衡模型发展到多层水量平衡模型[9-11]，主要涉及的作物有冬小麦、夏玉米、牧草等[12-14]，模型模拟的灌水方法有地面灌溉包括漫灌、沟灌、喷灌等全面灌水方法，但对于如何模拟局部灌溉对土壤水分的影响尚未研究。蓄水坑灌是一种典型的适用于果林的局部灌溉方式，水分通过坑壁侧向渗入根区土壤，入渗界面并非普通灌溉方式的垂向入渗，而是侧向入渗[15]，因此现有的分层水量平衡模型无法直接应用于蓄水坑灌土壤水分动态模拟，需要对现有的模型进行改进。基于此，本文将考虑蓄水坑灌的水分分布特征，对现有分层水量平衡模型进行改进，建立蓄水坑灌条件下的分层水量平衡模型，以期为蓄水坑灌果园土壤水分动态预测提供新途径。

1 材料和方法

1.1 试验区概况

试验于山西省晋中市太谷县农业科学院果树研究所（112°32′E，37°23′N）进行，该地平均海拔约781.9 m，无霜期175 d，平均气温一般在5～10 ℃之间。全年降雨量分布不均匀，全年降雨量最多的时期为6-8月，降雨量最少的时期为12-2月，多年平均降雨量约为460 mm，气候类型属于典型的温带性大陆气候。

依据试验区土壤质地对土体进行分层，将地表以下200 cm的土体分为3层，第一层厚度为70 cm，第二层厚度为50 cm，第三层厚度为80 cm。土壤的主要物理参数如表1所示。

表1 土壤物理参数表Tab.1 Table of physical parameters of soils

1.2 试验设计

试验所用果树品种为红富士，南北方向种植，行距4 m，株距2 m。试验用树选择树龄与树势基本一致、无病虫害的苹果树，并在树下布置蓄水坑，蓄水坑的直径为30 cm，深度为40 cm。本试验共设2个处理，分别为处理一（灌水上下限为田间持水量的50%～80%）和处理二（灌水上下限为田间持水量的60%～90%），每个处理3个重复。试验于2015年4-10月进行，如表2所示为试验期间灌水安排的具体情况。试验期间每隔7～10 d对各层土壤含水率进行监测，并在降雨、灌溉前后进行加测。

表2 灌水时间及灌水定额Tab.2 Irrigation time and irrigation quota

1.3 测试项目

1.3.1 土壤水分的测定

各测点土壤体积含水率采用TRIME-PICO IPH土壤水分测量系统进行监测，监测范围为果树的树干周围地表以下0～200 cm深度。

土壤水分测点分过坑、不过坑、处于过坑与不过坑的中间位置3种情况。过坑测点布置于树干与蓄水坑中心连线上，距离树干分别为30 cm（g30）、50 cm（g50）、100 cm（g100）、120 cm（g120）；不过坑测点布置于相邻两蓄水坑分界线与树干连线上，距离树干分别为30 cm（b30）、60 cm（b60）、90 cm（b90）、120 cm（b120）；中间位置测点距离树干分别为90 cm（z90）、120 cm（z120）。在垂直方向上每20 cm取一个水分测点，深度取到地表之下200 cm处，水平方向的水分测点布置如图1所示。

图1 含水率监测点位置示意图（单位：cm）Fig.1 Schematic diagram of water content monitoring points

1.3.2 叶面积指数的测定

在试验期间每隔一周采用LAI-2200冠层分析仪对叶面积指数进行测量。

1.3.3 气象数据的测定

本次试验采用精度较高的小型自动气象监测站对试验地的参考作物需水量ET0以及降雨量进行了实时监测，图2为试验期间ET0及降雨量变化图。

图2 试验期间ET0及降雨量变化图Fig.2 ET0 and rainfall variation during the experiment

2 模型建立

2.1 蓄水坑灌果园土壤分层水量平衡模型的建立

土壤水量平衡模型根据一定时段内水分的输入和输出来确定土壤水分的动态变化，反映了作物根系层水量变化与水分收、支之间的关系。根据蓄水坑灌的布置以及土壤质地的分层情况，将0～200 cm土层分为上、中、下3层，如图3所示为蓄水坑灌果园土壤分层水量平衡模型概化图。第一层土壤的水分收入有降雨量P去除植被截留量C之后的净雨量(PC)、灌溉量I和第二层土壤分配至第一层的水量q1,2，水分支出有地表蒸发量Es、蓄水坑的坑壁蒸发量Ep、下渗至第二层的水量Q1、根系吸水量S1；第二层土壤的水分收入有第一层的下渗水量Q1和第三层土壤分配至第二层的水量q2,3，水分支出有下渗至第三层的水量Q2、分配至第一层的水量q1,2和根系吸水量S2；第三层土壤的水分收入为第二层的下渗水量Q2，水分支出有第三层的下渗水量Q3、分配至第二层的水量q2,3和根系吸水量S3。不考虑地下水补给，结合土壤水量平衡模型可得蓄水坑灌果园的各层土壤水量平衡方程如下：

图3 蓄水坑果园土壤水量平衡模型概化图Fig.3 Soil water balance model generalization of water storage pit orchard

由各层土壤的储水量变化可得出每一层的土壤含水率计算公式：

2.2 下渗水量

蓄水坑灌是一种典型的局部灌溉技术，在布设蓄水坑并进行灌水之后，水分将通过坑的侧壁渗入土壤，土壤水分呈现出非均匀特征[16]，且表层土壤可能没有湿润或者仅很小部分湿润。因此，蓄水坑灌条件下降雨和灌溉之后的水分运移情形是不同的，降雨之后水分从地表开始渗入土壤，逐渐湿润下层土壤；灌溉时水分侧向入渗，并且水分并未完全湿润上层土壤就向下运动。因此在本模型中对于第一层土壤的下渗水量而言，降雨和灌溉需要分开考虑。

2.2.1 降雨之后的下渗水量

雨水经过植被拦截之后由地表向下层土壤入渗，根据满溢渗透的假定[17]，当土层中的水分达到田间持水量之后再向下渗透，直到降雨量和下渗水分渗透完为止。降雨之后超出饱和含水量的部分全部下渗，由于实际情形中即使没有达到田间持水量，也会发生土壤水分的下渗，所以本文认为未超出饱和含水量的部分当其达到临界含水量[18]时就发生下渗。

在实际计算下渗水量时，分为两种情况，一种是该层土壤的水分增加量超过该层的容水量即该层所能保持的水分（时段初土壤储水量与饱和含水率对应的土壤储水量的差值），此时各层的下渗水量公式计算如下[18,19]：

式中：Pi为第i天的降雨量，cm；为第i天第j-1层土壤下渗水量，cm；为第i天第j层土壤容水量，cm；n为土壤剖面的排水常数[18]，取值范围一般为0.30～0.35；θsj为第j层土壤的饱和含水率，cm3/cm3；θj为入渗之前第j层土壤的含水率，cm3/cm3；θfj为第j层土壤的田间持水率，cm3/cm3；LAI为作物叶面积指数；LAImax为某一时段内作物叶面积指数的最大值；θkj为第j层土壤的临界含水率（当土壤水分低于某一临界点时，导致作物的生长状况发生显著变化时所对应的土壤含水率[20]；θkj与作物叶面积指数及各层土壤的田间持水率呈线性关系[18]），cm3/cm3。

另一种是该层土壤的水分增加量并未超过该层的容水量，若此时该层的土壤含水率小于该层土壤的临界含水率时则下渗水量为0，当该层的土壤含水率大于临界含水率时按下式计算下渗水量：

2.2.2 灌溉之后的下渗水量

灌水之后，蓄水坑中的水并未完全湿润表层就向下运动，不满足满溢渗透假定，为了使模型反映此种现象，假定第一层土壤所能保持的水分即容水量减少，即第一层入渗水量超过容水量就向下层土壤入渗。因此，本文引入局部灌溉系数ξ对灌溉之后第一层土壤的容水量进行了修正，修正之后第一层土壤的容水量为ξM1(ξ＜1)。在计算灌溉之后土壤的下渗水量时只用修正第一层土壤0～70 cm土层的容水量，如公式（7）所示，剩余两层土壤的下渗水量仍按照降雨时的下渗水量公式进行计算。

2.3 土壤重分配水量

由于各层的土壤含水率差异较大，需考虑在相邻土壤层间扩散的水量，引入下述土壤水重分配公式[19]对重分配的水量进行计算：

式中：qj,j+1为第(j+1)层的水分再分配量（以水分向上扩散为正），cm；θwj为第j层土壤的凋萎含水率，cm3/cm3；ηj为土壤水分扩散系数，与含水量有关，按下式计算：

植物不能从土壤中吸水时的土壤含水率称为凋萎含水率[21]，本文蓄水坑灌果园土壤凋萎含水率对应的负压水头为-150 m[22]。

2.4 蒸发量

蓄水坑灌法是将灌溉水注入树冠下的若干个蓄水坑内，比起传统灌溉多了坑壁蒸发，蓄水坑灌的棵间蒸发包括地表蒸发和坑壁蒸发。在蓄水坑灌条件下，蒸发强度主要与土壤含水率因素有关，本文采用谷琼琼[23]拟合得到的蓄水坑灌果园地表蒸发强度和坑壁蒸发强度计算公式：

式中：es和ep分别为蓄水坑灌地表蒸发强度和坑壁蒸发强度，mm/d，通过已知初始含水率计算得到初始的地表蒸发强度和坑壁蒸发强度。

2.5 植被截留量

当发生降雨时，部分雨水经过植物的叶片时会被截留，这部分雨水不会渗入土体，超过截留能力的雨水才会到达土壤表面。植被的最大截留降水量指的是单位面积上植被所能截留的最大降水量，即聚积在叶面所在水平面上的水层深度。植被截留量[24]的大小主要与植被的叶面积指数LAI和降雨强度有关。叶面积指数为单位面积上的植物叶面积总和，即叶面积总和与占地面积的比值。本文采取马生军等人的苹果树冠层降雨截留模型[25]：

式中：C为截留量，mm；P为计算时段内的降雨量，mm。

2.6 根系吸水量

本文采用的根系吸水模型为[26]：

式中：S(x,y,z,t)为根系实际吸水量，L；Vmax(x,y,z,t)为根系最大吸水速率，L/min；γ(h)为水分胁迫系数。

γ(h)可表示为：

式中：h0、h1、h2、h3为θs、80%θf、60%θf、θw对应的基质势水头。Vmax(x,y,z,t)采用Vrugt[27]的表述。

式中：Tpot为潜在蒸腾强度，cm/min；Kc为作物系数，采用FAO[28]发布的作物蒸散量、作物需水量计算指南中的数据；ET0为参考作物需水量，cm/min，采用气象监测站采集的数据。

依据上述根系吸水模型和气象站收集到的数据，采用MATLAB进行编程，计算得出模拟区域各层各个节点处的根系吸水速率，根据求得的节点处的平均根系吸水速率进而求得各层的根系吸水速率以及根系吸水量。

2.7 模型精度评价指标

精度评价的指标为均方根误差RMSE、平均相对误差MAPE和平均绝对误差MAE，计算公式如下：

3 结果分析

3.1 模型率定与验证

将2015年4月13日测得的土壤含水率作为计算的初始含水率，采用上述的模型对蓄水坑灌条件下4-10月份的田间土壤含水率进行模拟。采用处理一的数据率定局部灌溉系数ξ，并采用处理二的数据对模型进行验证。系数ξ的率定过程为，采用分层水量平衡模型，分别计算ξ为0.6、0.7、0.8、0.9、1.0时的土壤含水率，通过精度评价比较得出在本模型下参数ξ的最优值，表3为不同ξ取值下处理一模拟计算精度评价表。

表3 不同ξ取值下处理一模拟计算精度评价表Tab.3 Evaluation of the accuracy of treatment 1 simulation calculation under different values of ξ

从表3的计算精度评价表来看，随着ξ的取值从0.6变化到1.0，模型的计算误差先减小后增大，当ξ=0.8时的模拟效果最好，所以ξ的率定值为0.8。同时由表3可知，模型各个土层深度的计算误差随深度增加基本呈先减小后增大的趋势，70～120 cm土层模拟精度最高。这可能是由于0～70 cm土层直接受到降雨、灌溉、地表蒸发等气象条件和蓄水坑边界等多因素影响，造成该层模型模拟的精度比70～120 cm低；120～200 cm土层在垂向上处于模型模拟的最深层，在模型中只考虑了70～120 cm 土层和120～200 cm土层之间的水分交换，忽略200 cm以下区域与120～200 cm土层之间的水分交换，造成该层模拟精度比70～120 cm土层低，但从表3的误差指标来看，忽略这一部分的水量虽然模拟精度降低但在可接受范围之内。进一步将ξ=0.8代入模型，进行处理二模拟计算，采用处理二的实测数据对模型进行验证，表4为处理二模型计算精度评价表。

表4 处理二模拟计算精度评价表Tab.4 Evaluation of the accuracy of treatment 2 simulation calculation

从表4来看，当设定局部灌溉系数ξ为0.8时，处理二各土层土壤含水率模拟值与实测值的3种误差都较小，模型模拟的精度较高，可以用于蓄水坑灌土壤水分动态模拟。

3.2 土壤水分动态变化对比分析

图4和图5是2015年4-10月蓄水坑灌果园的分层水量平衡模型模拟的土壤含水率随时间变化的动态。表5和表6是灌水和降雨之后各层土壤含水率模拟值与实测值的增加幅度情况。表7和表8是两种处理下各层土壤含水率的极差值以及平均含水率值。

图4 处理一各土层土壤含水率模拟值与实测值比较Fig.4 Comparison of simulated and measured soil water content values for each soil layer in Treatment 1

图5 处理二各土层土壤含水率模拟值与实测值比较Fig.5 Comparison of simulated and measured soil water content values for each soil layer in Treatment 2

表5 灌水之后各层土壤含水率模拟值与实测值的增加幅度Tab.5 Increase in simulated and measured soil water content of each layer after irrigation

表6 降雨之后各层土壤含水率模拟值与实测值的增加幅度Tab.6 Increase in simulated and measured soil water content of each layer after rainfall

表7 各层土壤含水率极差值Tab.7 Value of extreme difference of soil water content in each layer

表8 各层土壤平均含水率 cm3/cm3Tab.8 Average water content of soil in each layer

从图4和图5土壤含水率的动态变化来看，两种处理条件下3层土壤含水率模拟值与实测值的变化趋势较为一致，结合表3和表4的计算精度评价指标表明实测值与模拟值吻合较好。两个水分处理的水分动态模拟结果表明，由于蒸发、水分下渗、果树的根系吸水等耗水因素，随着时间的推移各层土壤含水率呈下降趋势，当有灌溉和降雨时呈上升趋势，总的来说各层土壤含水率在一定范围内不断波动。在灌水日期附近，两种处理下各层土壤含水率均出现了较大幅度的增加，符合实际情况。处理一在模拟期间于5月19日和7月12日进行了两次灌水，处理二在模拟期间于5月4日、6月20日和8月26日进行了3次灌水，灌水之后各土层土壤含水率的模拟值与实测值分别较上个测量时间的增加幅度如表5所示。综合表5的数据来看，灌水之后70～120 cm和120～200 cm土层土壤含水率模拟值与实测值的增加幅度明显大于0～70 cm土层，且70～120 cm土层的增加幅度最大，这是由于蓄水坑灌是中深层立体的局部灌溉方式，灌溉水更集中分布在整个土体的中层，这也是果树根系分布相对密集的区域。说明采用蓄水坑灌的灌溉方式可以有效减少地表蒸发，给果树生长提供更充足的水分。在整个模拟期间降雨量较大的时段出现在9月8日-10日，累计降雨量为57.4 mm，这次降雨之后两种处理下各层土壤含水率模拟值与实测值也出现了较大幅度的增加，具体增加幅度如表6所示。从表6的数据来看，降雨之后0～70 cm土层土壤含水率模拟值与实测值的增加幅度明显大于70～120 cm和120～200 cm土层，这是由于降雨增加的水量是直接从地表开始入渗的，因此上层土壤含水率的增加幅度大于中下层土壤。在蓄水坑灌条件下灌溉和降雨之后0～70 cm土层的下渗机理不同，因此灌溉和降雨发生之后从浅到深各土层土壤含水率的增加幅度呈现了不同的趋势，具体表现为灌溉之后从浅到深各层土壤含水率增加幅度呈先增大后减小的趋势，降雨之后从浅到深各层土壤含水率的增加幅度呈逐渐减小的趋势。从表5及表6的数据来看，灌水和降雨之后各土层土壤含水率实测值的增加幅度通常低于模拟值，这是由于实际水分的变化具有一定的滞后性，没有模型模拟对灌水和降雨的响应迅速。

从3个土壤层水分的动态模拟情况和表7的数据来看，最接近地表的0～70 cm土层由于受到外界气象因素的影响较大，比如降雨、灌溉以及地表蒸发和坑壁蒸发，相比较其他土层，0～70 cm土层的含水率变化趋势波动幅度较大，70～120 cm土层以及120～200 cm土层含水率波动趋势相较上层土壤更为平缓。从表8的数据来看，模拟期间两个处理的平均含水率均表现为0～70 cm土层最低，这是因为表层直接受到太阳辐射，地表蒸发和坑壁蒸发较为强烈。处理二的灌水上限高于处理一，且处理二的灌水次数多于处理一，处理二各土层的土壤平均含水率均高于处理一。

4 讨论

土壤分层水量平衡模型已广泛应用于农田各层土壤水分动态的模拟，陈皓锐[19]、张雪飞[29]等人根据水量平衡原理建立了土壤分层水量平衡模型，模拟了各土层土壤含水率的动态变化。但现有的分层水量平衡模型无法直接应用于蓄水坑灌土壤水分动态模拟，本文基于土壤的水量平衡原理，在现有模型的基础上引入局部灌溉系数ξ，建立了蓄水坑灌条件下0～200 cm土层的分层水量平衡模型，通过精度评价，表明模型的精度较高，证实本文的分层水量平衡模型可适用于蓄水坑灌果园土壤水分动态模拟。局部灌溉系数ξ的合理取值十分关键，ξ的大小会影响到模型的下渗水量，在模拟计算前需要对此参数进行率定。

由于0～70 cm土壤与外界的水分交换较为频繁，本模型模拟的该层土壤含水率变化对外界的响应更敏感，波动幅度更明显，这与刘雪静[30]等人的研究一致。由于测量仪器在浅层更容易产生误差以及灌溉和降雨会引起含水率的突变，导致该层土壤的模拟误差较大。表层土壤直接受到日光照射和温度的影响，水分流失较快，模型模拟易高估0～70 cm土层土壤含水率，这与刘昭[10]等人的研究结果一致。下层土壤最接近地下水，在实际中即使超过了田间持水率也有可能并不会立即发生水分下渗，因此在用本模型模拟时会出现低估120～200 cm土层土壤含水率的现象，这与郭其乐[9]等人的研究结果一致。当长时间无降雨或者灌溉，无明显的水分收入之后，模型的模拟能力会下降[9]，在降雨以及灌溉之后，模型的水分收支情况更新，会使各层土壤水分不断达到新的平衡，使模型的模拟效果更好。