华北山前平原区厚包气带次降水入渗补给定量估算*

刘美英 ,高 雅 ,张俊芝 ,沈彦俊 ,吴 林 ,李明良** ,闵雷雷**

(1.中国科学院遗传与发育生物学研究所农业资源研究中心/中国科学院农业水资源重点实验室/河北省节水农业重点实验室石家庄 050022;2.中国科学院大学 北京 100049;3.河北省水文勘测研究中心 石家庄 050051)

降水入渗补给是水文循环的重要环节之一,是地下水资源的重要组成部分,华北平原降水入渗补给占浅层地下水总补给量的70%左右[1],因此深入研究降水入渗补给过程、正确评估降水入渗补给量是地下水资源评价、可持续利用和科学管理的基础[2-3]。另一方面,在我国区域水资源评价工作中,年降水入渗补给量是由年降水量乘以年降水入渗补给系数得出的。但在实际条件下,年内的补给并非是随着降水同步均匀发生的,而往往只是由年内几场降水引起;由于不同年份降水年内分配具有较大差异,即使年降水量相同,年降水入渗补给量也会存在一定差异[4],因而年降水入渗补给系数无法反映降水年内分配不均导致的补给量差异。开展次降水入渗补给过程研究,深入分析次降水入渗补给机理,定量计算次降水入渗补给量,才能在机理上推动区域入渗补给量评估的科学性和准确性[5]。

次降水入渗补给是指一次大气降水补充到地下水的水量,该水量与本次降雨量的比值称为次降水入渗补给系数[6-7]。次降水入渗补给受到诸多因素的影响,例如: 降水条件(降水总量、降水强度、降水历时等)、地形地貌、地表植被类型、包气带岩性、降水前期土壤含水量和地下水埋深等。目前,学者针对次降水入渗补给开展了大量研究。部分学者研究了地下水埋深、次降水特性和包气带岩性对次降水入渗补给量的影响[6,8-11]及次降水入渗补给系数的空间变异性[7]。也有部分学者研究了次降水入渗补给机理,并重点分析了次降水入渗补给过程随地下水位埋深增加而发生的坦化和滞后现象[4-5,12-20]。以上研究表明: 降水入渗补给具有滞后性,且滞后时间随地下水埋深增加而增大,补给过程线随地下水埋深增加而坦化严重,埋深大于一定程度时,多个降水场次形成的补给过程逐渐合并成为一个总的补给过程。当前,土壤水均衡法、蒸渗仪法、地下水位动态法、环境组分(同位素)示踪法、数值模拟法等手段广泛应用于年降水入渗补给量的估算[21-25]。但在地下水位埋深较大的地区,由于包气带较厚,利用蒸渗仪等观测手段直接测定实际补给量不太现实,地下水位动态法也受到补给事件以外因素的干扰而无法使用,因此通常采用观测或者估算的潜在补给量(深层渗漏量)来代替实际补给量。总之,在地下水深埋区,究竟多少深度土层的渗漏量能够代表潜在补给量、如何明确次降水入渗补给事件和确定相应的潜在补给量以及次降水入渗补给量与次降水量间存在怎样的关系是亟需厘清的问题,但目前相关的研究还比较少。

本文以华北山前平原典型植被类型区[冬小麦(Triticum aestivum)-夏玉米(Zea mays)一年两熟轮作]为例,基于典型站点的厚包气带观测数据,深入研究了地下水深埋区次降水入渗补给机理,定量计算了次降水入渗补给量,并建立次降水入渗补给量与次降水量间的定量关系,研究结果可为区域地下水水资源管理提供理论依据。

1 材料与方法

1.1 研究区概况

华北平 原位于113.172°～119.845°E,34.981°～40.777°N。山前平原位于华北平原西部和北部,沿太行山、燕山山麓呈带状分布,宽度30～60 km,由形成于不同时期、规模不等的河流冲洪积扇叠加和连接而成,主要为第四系冲积的砂砾石层、亚砂土及亚黏土的松散堆积物,包括大清河水系山前平原和子牙河水系山前平原等在内的多个亚区[26]。

华北山前平原气候属于暖温带半湿润大陆性季风气候,全年7 月温度最高,l 月温度最低;雨热同期,降水量主要集中在7、8 月份[27]。区域内,冬小麦-夏玉米一年两熟轮作是主要种植制度[28],冬小麦一般10 月上旬播种,翌年6 月上旬收获;夏玉米一般在冬小麦收获后播种,9 月下旬收获。由于降水量无法满足一年两熟制的作物需水,作物高产依赖地下水灌溉,冬小麦生育期内一般灌溉3～5 次,夏玉米生育期一般灌溉1～2 次,每次灌溉量约60～100 mm,高强度开采地下水导致地下水位持续下降,包气带不断增厚。

1.2 试验方案与数据观测

本研究依托华北山前平原两个试验站开展的次降水入渗补给试验,分别为位于子牙河流域的中国科学院栾城农业生态系统试验站(以下简称“栾城站”)和位于大清河流域的河北省水文局冉庄水资源实验站(以下简称“冉庄站”) (图1a)。两个站点年平均降水量分别为506.7 mm和466.1 mm (1992—2016 年),研究时段内两个站点的包气带厚度均大于8 m;两个站点的沉积物类型分别为亚黏土和亚砂土,在山前平原具有较好的代表性;种植制度均为冬小麦-夏玉米一年两熟轮作,代表了山前平原的主要种植制度。

图1 栾城站和冉庄站位置 (a)及栾城站土壤水分要素观测竖井 (b)和冉庄站地中蒸渗仪(c)示意图Fig.1 Locations of Luancheng Station and Ranzhuang Station (a),schematics of soil water item observation instrumentation in Luancheng Station (b) and lysimeter in Ranzhuang Station (c)

栾城站地理位置为37°53′N、114°41′E,位于河北省栾城区,海拔50.1 m,2016 年地下水位埋深为45 m。该站点的田间试验水肥管理均参照当地平均水平,具有区域代表性;监测数据包括0～8 m 深土壤含水量、土壤水势。土壤含水量监测深度间隔为0.2 m;土壤水势0～1 m 监测深度间隔为0.2 m,1～8 m监测深度间隔为1 m (图1b)。土壤含水量和水势监测的时间间隔均不超过10 d。此外,该站内还对降水量(1985—2016 年)、灌溉量(1985—2016 年)和蒸散发(2007—2016 年)等数据进行观测。

冉庄站地理位置为38°40′N、115°23′E,位于河北省清苑县冉庄镇,距离保定市40 km,2016 年地下水位埋深为20 m。冉庄站系列监测数据依托站内大型地中蒸渗仪观测,该蒸渗仪于1987 年建成,面积为10 m2,地下水埋深为8 m,蒸渗仪中设有中子仪测管和负压计(图1c),用于观测土壤含水量和土壤水势,蒸渗仪还设有补给量观测装置。蒸渗仪内田间水肥管理也参照当地平均水平。土壤含水量监测深度为0～8 m,深度间隔为0.2 m,时间间隔为10 d 左右;土壤水势监测深度为0～5 m,间隔为0.2 m,时间间隔为10 d左右;补给量监 测深度 为0.3～8 m,间隔为0.7～2 m。土壤含水量、土壤水势和补给量的监测时间段为1988—2016 年,其中1988—1998 年数据资料较好,用于本研究。

1.3 研究方法

1.3.1 长时间序列降水入渗补给量模拟

本研究利用Hydrus-1D 模型模拟降水入渗条件下土壤水分运动和降水入渗补给量,Hydrus-1D 模型是由美国盐土实验室于1991 年开发的用于模拟土壤水流和溶质运移的数值模型软件[29]。该模型由改进的一维Richards 方程模拟水分运动,采用Feddes 模型模拟根系吸水[29]。本研究采用均方根误差(RMSE)和平均绝对误差(MAE)对模型的拟合效果进行评价。模型的详细模拟方程、初始和边界条件设置及模拟效果评价方法见附录(见官网www.ecoagri.ac.cn 此文的“资源附件”)。

Min 等[30]已对栾城站的模型参数进行了校准和验证,本研究采用已校验模型参数进行了长时间序列模拟。在冉庄站,利用土壤含水量、土壤水势和补给量等实测数据对模型参数进行校准和验证,结果表明: 基于Hydrus-1D 模拟的土壤含水量、土壤水势、土壤储水量及8 m 深水分通量与实测值均吻合较好,模型率定详细结果见附录(见官网www.ecoagri.ac.cn 此文的“资源附件”)。

1.3.2 数据处理与分析方法

本研究基于Microsoft Excel 2020 和Origin 2021处理数据并绘制图表,采用SPSS 21.0 软件进行数据统计分析。

2 结果与分析

2.1 潜在补给界面发生深度的确定

土壤中水势变化是水分下渗和水分耗损协同作用的结果,土壤水势梯度是决定水分运动方向的驱动力[31-32]。定义水势梯度为垂直剖面上两深度之间水势(φ)之差(φ1-φ2)与两深度(h)之间距离(h1-h2,h1位于点h2下方)的比值。当水势梯度为正值时,代表自h2深度以下的土壤水分运动方向向下;当某一层水势梯度在研究时段内均大于0 时,则表明该层以下土壤的水分仅向下运动而不再向上运动,该深度也可称为零通量面出现的最大深度。

利用栾城站和冉庄站土壤水势观测资料刻画土壤水势梯度变化特征。在栾城站,当土层深度大于2 m后,土壤水势梯度均大于0,即地表以下2 m 为该站点零通量面存在的最大深度,该深度以下土壤水分不再向上运动仅向下运动(图2a,b,c);在冉庄站,当土层深度大于3 m 后,土壤水势梯度才基本上大于0,即地表以下3 m 为该站点零通量面存在的最大深度(图2d,e,f)。值得注意的是,冉庄站土壤类型为亚砂土,栾城站土壤类型为亚黏土,而通常土壤质地越细,蒸发极限埋深(即零通量面存在的最大深度)越大,但冉庄站的零通量面存在的最大深度却大于栾城站,这可能是由地表水分输入条件差异引起的。栾城站1992—2016 年平均降水量和灌溉量分别为506.7 mm·a-1和358.7 mm·a-1,冉庄站1992—2016 年平均降水量和灌溉量分别为466.1 mm·a-1和355.5 mm·a-1,栾城站总水分输入比冉庄站多43.7 mm·a-1,进而导致栾城站蒸散发影响到的土层深度相对较浅。

图2 栾城站(a,b,c)和冉庄站(d,e,f)土壤水势梯度随深度变化(正值代表水分向下运动)Fig.2 Variations of the total soil water potential gradients along depth in Luancheng Station (a,b,c) and Ranzhuang Station (d,e,f)(positive values represent the downward movement of soil water)

利用校准和验证完成的Hydrus-1D 模型,模拟了栾城站和冉庄站长时间序列(1992—2016 年)地表以下不同深度水分通量(图3)。分析两个站点水分通量随深度变化特征可发现: 零通量面存在的最大深度以下土层年平均水分通量随深度波动极小。在栾城站,零通量面存在的最大深度(2 m)以下水分通量趋于不变(图3a)[30];在冉庄站,零通量面存在的最大深度(3 m)以下水分通量随深度变化也很小(图3b)。

图3 栾城站(a)和冉庄站(b)土壤层年平均水分通量随深度变化Fig.3 Variations of average annual recharge fluxes along depths in Luancheng Station (a) and Ranzhuang Station (b)

尽管零通量面存在的最大深度以下土层年平均水分通量随深度波动极小,但随着深度增加,深层补给的滞后和坦化效应更明显。例如,在栾城站模拟的2000—2002 年2 m、5 m 和8 m 水分通量结果显示2 m 深度处水分通量变化较为剧烈,而5 m 和8 m 深度处水分通量波动状况依次减缓,且5 m 和8 m 深度处水分通量峰值相较于2 m 均存在不同程度的滞后(图4a);在冉庄站模拟的1994—1996 年3 m、6 m和8 m 水分通量结果表明,3 m 处水分通量变化较为剧烈,而6 m 和8 m 水分通量波动状况依次减缓,且6 m 和8 m 水分通量峰值相较于3 m 均存在不同程度的滞后(图4b)。因此,考虑到深度越大包气带对补给过程的调蓄和缓冲作用越强,补给过程越平缓,更不利于次降水入渗补给量的确定。因此,本研究将2 m 和3 m 分别定为栾城站和冉庄站地下水潜在补给的发生深度,并利用这两个深度的水分通量开展次降水入渗补给量分析。

图4 基于Hydrus-1D 模型模拟的栾城站(a)和冉庄站(b)不同深度地下水补给过程Fig.4 Recharge process at different depths in Luancheng Station (a) and Ranzhuang Station (b) simulated by Hydrus-1D model

2.2 次降水入渗补给量的确定

本文考虑作物类型、包气带特性和降水条件等因素,提出了适宜的次降水入渗补给事件识别的标准: 连续7 d 无降水为次降水的划分间隔,日平均降水量大于2.00 mm (详细划分标准见附录文件,官网www.ecoagri.ac.cn 此文的“资源附件”)。根据该标准将研究时段内(1992—2016 年)栾城站和冉庄站分别划分出67 场和69 场次降水入渗补给事件,且确定了相应的降水入渗补给量。栾城站67 场次降水入渗补给量为0.1～421.7 mm,平均值为22.2 mm,标准差为57.7 mm;冉庄站69 场次降水入渗补给量为0.4～261.9 mm,平均值为32.0 mm,标准差为65.1 mm。

在确定了次降水入渗补给量的基础上,分析次入渗补给量与次降水量以及前期土壤含水量影响下次总有效水量(次降水量与降水前期0～200 cm 土层储水量的和)间的关系,并进一步分析次降水入渗补给系数与上述两者水分输入量之间的关系。

栾城站和冉庄站次降水入渗补给量与次降水量均存在较好的相关性,相关系数(R2)分别为0.894 和0.801 (图5a,c);当考虑土壤前期含水量影响时,次降水入渗补给量与总有效水量相关关系更佳,R2分别为0.962 和0.853 (图5b,d)。这是因为土壤前期含水量对降水入渗补给量具有重要影响,降水发生时,如果降水前期土壤含水量很高,甚至在饱和水平,则土壤剖面仅需少量水分补充甚至无需补充就可形成补给,而当降水前期土壤含水量很低时,降水首先需补充土壤剖面水分亏缺,然后再形成补给,如果降水量很低,甚至不足以补充水分亏缺,则不会形成补给。

图5 栾城站(a,b)和冉庄站(c,d)次降水量和次总有效水量(次降水量+前期200 cm 土壤储水量)与次降水入渗补给量的关系Fig.5 The relationship between sub-precipitation amount and total effective sub-water input (sub-precipitation amount and previous soil water storage in 200 cm soil profile),and sub-precipitation infiltration recharge amount in Luancheng station (a,b)and Ranzhuang station (c,d)

栾城站和冉庄站两个站点次降水入渗补给系数与次降水量相关关系较差,R2分别仅为0.073 和0.192(图6a,c);如果考虑土壤前期含水量的影响,总有效水量与次降水入渗补给系数相关关系虽比不考虑前期土壤含水量有所改善,但仍然较差(栾城站和冉庄站R2分别为0.128 和0.260,图6b,d)。这是因为降水入渗补给系数是次降水入渗补给量与次降水量的比值,其分母并未考虑降水前期土壤含水量的影响。本研究中降水入渗补给量为地表以下2 m 或3 m 水分通量,如果降水前期土壤含水量较低,水分亏缺比较严重,那么即使降水量较大,入渗补给量也可能很低甚至没有,则入渗补给系数将会很小甚至为零;而如果降水前期土壤含水量较高,甚至高于田持水平,那么即使降水量相对较小,也会有较大的补给量,则入渗补给系数将较大甚至达到1,这是次降水量和次降水入渗补给系数相关关系差的主要原因。

图6 栾城站(a,b)和冉庄站(c,d)次降水量和次总有效水量(次降水量+前期200 cm 土壤储水量)与次降水入渗补给系数的关系Fig.6 The relationship between sub-precipitation amount and total effective sub-water input (sub-precipitation amount and previous soil water storage in 200 cm soil profile),and sub-precipitation infiltration recharge coefficient in Luancheng station(a,b) and Ranzhuang station (c,d)

综上所述,利用统一的降水划分标准对次降水场次进行划分,结合合适的次降水入渗补给量确定方法,再进行次降水事件的合理筛选和次降水入渗补给量的校核,可进行次降水入渗补给机理研究。研究发现: 研究时段内次降水入渗补给量具有较大变异性,主要受到水分输入条件的影响,次降水入渗补给量与次降水量和次总有效水量(次降水量+前期200 cm 土壤储水量)存在显 著相关关系,R2可达0.801～0.962。但次降水入渗补给系数与次降水量和次总有效水量(次降水量+前期200 cm 土壤储水量)的相关关系较差,两个站点的R2为0.073～0.260。这是因为次降水入渗补给系数没有考虑降水前期土壤含水量对补给的影响,因此在实际应用中,利用次降水入渗补给系数可能无法准确计算次降水入渗补给量。

3 讨论

次降水入渗补给研究是年或汛期降水入渗补给研究的基础,开展次降水入渗补给机理和定量研究对于提高年降水入渗补给量的估算精度具有重要意义。本研究利用土壤水分、土壤水势和补给量等数据对Hydrus-1D 模型进行校准和验证,然后模拟了长时间序列降水入渗补给通量,进一步构建了适宜的次降水划分标准,深入分析了次降水入渗补给机理,研究了次降水入渗补给量及其与次降水量之间的关系。

次降水事件的合理划分是开展次降水入渗补给研究的前提,研究发现次降水划分受流域特性和地下水位埋深[4-5,33]等因素的影响,其划分标准可能具有不可移植性,划分需视具体研究区域和研究目标而定。已有研究表明,降水入渗补给具有明显的滞后性,且如果两次降水事件相隔时间较短时,入渗补给过程线有可能重叠甚至合而为一,只形成一次补给事件。因此,对次降水进行划分时,需要考虑次降水的时间间隔和当次降水的持续时间[33],一方面要保证次降水划分间隔足够长,即要考虑补给的滞后效应和持续时间[13,34]的影响,避免相邻几次降水事件的补给过程相互叠加无法区分;另一方面划分间隔也不能太长,否则将出现每年只有1～2 场降水事件,且大多数降水事件中包括灌溉事件,导致无法开展次降水入渗补给机理的研究。本研究次降水事件划分从补给过程线入手,在总结前人相关研究成果和区域特性引起的相关限制条件下,通过试错法确定了适用于本地区的次降水划分标准,划分标准对相关研究具有指导意义。

本研究中的降水入渗补给量并非是实际地下水补给量,是地表以下2 m 或3 m 深度水分渗漏量,即潜在降水入渗补给量,但是该补给量可近似等于实际地下水补给量,在地下水深埋条件下,实际地下水补给滞后时间更长且更坦化[13-14,18,35],本研究确定的潜在补给界面发生深度对相关研究具有借鉴意义。

本研究栾城站和冉庄站分别划分出67 场和69场次降水事件,且这些降水场次主要发生在汛期,即年内补给主要是由汛期几次降水事件引起的,这与以往研究结果相似,例如李娜等[21]利用HYDRUS-1D 模型模拟了北京大兴区厚包气带水分通量(1951—2005 年,55 年),结果表明年内入渗补给主要发生在6—10 月份。

次降雨入渗补给过程极其复杂,本研究中提出的次降水划分标准和次降水入渗补给量确定方法是针对冬小麦-夏玉米灌溉农田确定的,次降水划分标准受到土壤质地、地下水埋深、降水、植被类型等综合影响,该研究中的方法和结果可为相关研究提供参考,但需要视流域特性等进行调整。

仍需要注意的是,尽管本文在研究降水入渗补给时剔除了灌溉影响场次,但灌溉对补给的形成起到很大的促进作用。在栾城站,Pei 等[36]研究冬小麦、夏玉米雨养条件下补给量时发现,由于没有春季灌溉对土壤剖面水分亏缺的补充,玉米季汛期补给量很低,雨养条件下年平均补给量(2011—2015 年)仅为1.3 mm·a-1,远小于灌溉条件下年平均补给量(116 mm·a-1),因此,灌溉为夏季的降水入渗奠定了更好的土壤水分基础。当前对于灌溉和降水对补给的各自贡献仍不清楚,仍需要深入开展降水和灌溉的补给机理研究。

4 结论

本文基于土壤含水量、土壤水势和蒸渗仪补给量等观测数据,结合Hydrus-1D 模型模拟了长时间尺度日降水入渗补给通量。结合区域特征并总结前人相关理论基础建立了次降水和次入渗补给的划分标准,次降水间隔定为7 d,该划分标准具有一定的物理意义。基于上述方法在栾城站和冉庄站(1992—2016 年,25 年)分别划分出67 场和69 场次降水入渗补给事件,次降水入渗补给量具有较大变异性。次降水入渗补给量与次降水量呈现较好的线性相关性,当考虑前期土壤含水量影响时,相关关系更佳。但是次降水入渗补给系数与次降水量相关关系较差,这是因为次降水入渗补给系数为次降水入渗补给量和次降水量的比值,无法反映降水前期土壤含水量对补给的影响,因此,利用降水入渗补给系数可能无法比较正确地估算次降水入渗补给量。