黄土高原地区因其特有的地理环境以及土壤质地特性,成为全球土壤侵蚀和水土流失最为严重的地区之一[1]。降雨引起的水土和养分流失加剧了土地退化和面源污染,因此需要更好地认识和量化分析径流的产生和土壤侵蚀过程。准确地预测降雨径流过程对进一步预测土壤侵蚀至关重要[2]。针对这一现状,国内外专家学者们开展了大量的野外试验和模型模拟研究,一些与径流过程和泥沙输移过程相关的数学模型已经被用来预测降雨条件下的水土流失过程。通常使用Saint-Venant方程(即连续性方程和动量方程)来描述地表径流过程[3]。但由于Saint-Venant方程是高度非线性的,因此很难获得解析解,这意味着只能使用数值方法对其进行求解。但当忽略Saint-Venant方程的加速度项时,可以使用扩散波方程来对其进行简化[4]。若同时忽略Saint-Venant方程的加速度和压力项时,Saint-Venant方程可表示为运动波方程。Luce等[5]通过使用Philip入渗方程修改了运动波方程来预测超渗产流的降雨过程,并将该模型应用于野外试验;Tsai等[6]使用特征值法与三次样条插值法相结合的方法来计算一维和二维运动波模型。

现有的土壤侵蚀模型一般可分为经验统计模型和物理模型。经验统计模型主要有USLE模型和RUSLE模型[7,8]等,因其考虑的影响因子全面且计算简单而被广泛采用,但经验统计模型并不能描述土壤侵蚀产沙的物理过程;现有物理模型主要是以泥沙运输机制为基础,结合水力学、土壤水动力学等来描述侵蚀产沙过程,有WEPP模型、EROSEM模型、汤立群模型和陈国祥模型等。WEPP模型[9]是一种基于侵蚀过程的模型,可以模拟气候、水文等变化与产沙之间的关系;李文杰等[10]将坡面水文模块和侵蚀产沙模块嵌套,构建了涵盖植被截留、入渗和径流等多因素的一维坡面侵蚀产沙模型;戚隆溪等[11]用特征线法和有限差分法求解了径流量和土壤侵蚀量;Haster等[12]通过实测资料和土壤侵蚀模型确定了侵蚀模型参数a1=0.155和b1=0.110;陶汪海等[13]采用Preissmann四点差分格式求解了运动波模型和产沙模型,从而得到坡面降雨条件下的产流过程,并进一步模拟了各种初始条件下的产流产沙过程;邹宏荣等[14]耦合新安江模型与土壤侵蚀概念模型,构建基于网格的淮河上游泥沙负荷模拟模型;Tao等[15]建立基于运动波方程和稳态产沙方程的半解析模型,求解了坡面产流产沙过程,并通过参数分析确定了侵蚀模型参数的取值。

本文通过耦合降雨条件下的Green-Ampt模型累积入渗量显函数、运动波方程和稳态的输沙方程,建立了基于入渗、产流和产沙等多水动力过程的坡面降雨产沙近似解析模型,通过设置不同雨强的模拟降雨试验对模型进行检验,并对模型参数进行敏感性分析,以期为黄土区坡面产沙过程预测和模拟提供一种可靠且便捷的方法。

1、理论与模型

1.1降雨条件下的Green-Ampt模型

通常情况下用Green-Ampt模型来描述降雨条件下的土壤入渗过程[15],其表达式为:

式中:i为入渗率(cm/min);ks为饱和导水率(cm/min);Sf为湿润峰处基质吸力(cm);I为累积入渗量(cm);θs为土壤饱和含水量(cm3/cm3);θ0为土壤初始含水量(cm3/cm3);t为入渗时间(min)。

由于公式(2)是关于入渗时间t的显式方程,而对于累积入渗量I则是隐函数方程,其求解需要借助数值分析方法求其近似数值解[16],为了简化其计算过程,本文使用刘继龙等[17]借助Philip入渗模型和Green-Ampt入渗模型[18]推导Green-Ampt模型的累积入渗量显函数的方法来求解累积入渗量I,表达式为:

式中:b是常数,主要受土壤水扩散率函数的形状控制,取0.55[19]。

当地表积水很浅或者入渗历时较长时,积水深度H可以略去,公式(3)可以简化为:

Mein等[20]在1973年提出了基于传统Green-Ampt模型的稳定降雨入渗过程的模拟方法;降雨初期,降雨强度小于土壤的入渗能力,土壤入渗率等于降雨强度,随着降雨过程的推移,地面开始产生地表径流,其过程可以表示为:

式中:p为降雨强度(cm/min);tp为降雨开始到开始产流所经历的时间(min);

降雨过程中地表产流之前,土壤入渗率等于雨强,即:

式中:t1为假设t=0时刻开始积水,累积入渗量I=Ip所用的时间(min)。

当雨强较大时,降雨至开始产流所经历的时间tp较小,假设降雨条件下产流开始后的入渗过程线即为将积水入渗条件下的入渗过程线向右平移t0个单位,公式(4)可以改写为:

式中:h为径流深(cm);q为单宽流量(cm2/min);

由于坡面水深与入渗率之间存在关联关系[25],径流的水量是降雨超渗所产生的,Yang等[26]假设坡面水深变化率与入渗率之间呈线性关系,即:

根据质量守恒原理,单位宽度小区在开始产流后任意时段内水量平衡关系为:

p(t-tp)l=Ql+Il+Hl(21)

式中:Ql为累积径流量(cm3/cm);Il为累积入渗量(cm3/cm);Hl为坡面累积积水量(cm3/cm);l为坡长。其中Ql、Il、Hl可以分别表示为:

式中:s为水流体积含沙率(g/cm3);Di为雨滴动能侵蚀率(g/(min·cm2));Dr为径流侵蚀率(g/(min·cm2));ρ为土壤密度(g/cm3);a为径流侵蚀系数(g/(min·cm3));b为雨滴侵蚀系数((g2·min)/cm7)。

公式(25)为非稳态产沙过程,为了方便计算,一般情况下假定产沙过程是恒定的[28],故忽略时间项,将其表示为:

2、材料与方法

2.1试验区概况

试验于2019年5月在中国科学院水利部水土保持研究所长武黄土高原农业生态试验站担水沟流域野外模拟降雨小区(35°12′N,107°10′E)进行,试验区海拔为1200m,气候属暖温带半湿润大陆性季风气候,年平均气温为9.1℃,年平均降水量580mm,典型土壤为粉砂质壤土,具体土壤物理特性见表1。母质为深厚的中壤质马兰黄土,流域塬面面积占35%,梁坡占35.6%,沟谷占29.4%,各约占1/3;地貌属典型的高原沟壑区。无霜期171d,地下水位50～80m,无灌溉条件,属典型的旱作雨养农业区。该区人口密度大,人均耕地少,水土流失严重,降水分布不均,干旱频繁,作物产量不高[29]。

表1试验区土壤物理化学特性导出到EXCEL

2.2试验方法

试验用地为3年撂荒地,为了便于控制和减小环境影响因素,本次模拟降雨试验的小区设置尺寸为1.0m×1.0m,坡度为15°;使用西安理工大学自主研制的针孔式人工模拟降雨装置(图1)进行降雨试验,针孔式降雨装置主要由主体支架(可调节高度),含有针孔的水槽(可根据雨强更换不同孔径的针孔),供水装置3部分组成,有效降雨面积为1.0m2,每次降雨试验开始前需反复率定雨强,通过控制水槽内水位高度来控制雨强大小。经过测试:该套人工降雨器的平均雨滴直径为2mm,降雨均匀度在80%以上,雨滴终速符合天然降雨特征。

图1针孔式降雨器示意

为了消除土壤前期含水率对试验结果的影响,在每次开始降雨24h前以25mm/h的雨强在试验小区进行预降雨,直至开始产流时停止降雨。当晴天蒸发作用较强时,在小区表面覆盖遮阳网,以防止表层土壤过分失水;若发生自然降雨事件,则使用防雨布对小区进行覆盖处理;为保证每次试验均保持相同的初始条件,开始正式降雨试验前测定小区内土壤表层0-20cm剖面的初始含水率,使用EC-5传感器测得体积含水率为15%±0.5%时开始试验。试验开始前在小区内覆盖防雨布,防止雨水进入小区内,反复率定雨强,通过调节水池中水位高度来调节雨强大小,直至其稳定为设计雨强时开始试验。根据试验区侵蚀性降雨实测资料及降雨分级标准[12],设计5种试验雨强分别为30,45,60,75,90mm/h,设计总降雨时间为60min,开始产流后按照前10min每1min1次和10min以后每2min1次的频率用量筒承接出口处径流,并采用过滤风干称重的方式测定径流泥沙含量。

2.3模型评价指标

用相关系数R2,均方根误差RMSE和纳什效率系数NSE对模型进行评价,其表达式具体为:

式中:yi为第i个试验测量值;为n组试验测量值的平均值;xi为第i个模型模拟值;n为试验数据的个数。

2.4数据处理

所有试验实测数据均为2次重复试验的平均值,使用Matlab2015b进行参数求解和模型模拟;使用SPSS24.0软件进行数据分析,使用Origin2017和MicrosoftExcel2016软件进行图表绘制和函数拟合。

用Dif(计算值与实测值差值平方的均值)来判断预测的准确性[30],表达式为:

式中:pi为第i个实测值;oi为第i个计算值;M为试验观测数据的个数。

采用相关系数R2,均方根误差RMSE对2次重复试验数据进行分析,由表2可知,2次重复数据的差异较小,进一步求2次平均值进行模型模拟分析。

表2实测2次单宽流量、产沙率的R2和RMSE导出到EXCEL

3、模型验证及参数敏感性分析

3.1径流过程模拟

通过不同雨强下的径流实测数据分别推求入渗率参数c和曼宁糙率系数n,采用相关系数(R2)和均方根误差(RMSE)对模拟结果进行评价。从表3可以看出,曼宁糙率系数n随雨强的变化均分布在0.00026～0.00029,而入渗率参数c随着雨强的增大呈减小趋势,且分布在0.0031～0.0060;R2均在0.89以上,均方根误差RMSE很好地刻画了实测值与计算值之间的偏差,随着雨强的增大,RMSE也随之增大,取值分布在0.406～1.052,NSE均大于0.397;而当雨强为60mm/h以上时,NSE则进一步增大到0.783以上,说明随雨强增大模型计算值和实测值的匹配性也随之变得更好;进一步揭示入渗率参数c随雨强的变化关系,采用指数函数对其进行拟合,从图2可以看出,指数函数可以很好地拟合参数c与雨强之间的关系,相关系数R2为0.9765,表达式为:c=9×10-4p-0.644R2=0.9765(35)

图2入渗率参数c随雨强的变化过程

将30～90mm/h的5个雨强梯度下的单宽流量实测值与模型计算值分别绘于图3。可以看出本文建立的模型可以较好地模拟设定雨强范围内的坡面产流过程,且随着雨强的增大,模拟趋势也随之愈加准确。整体来看,模型模拟的开始产流时间均较实际有所滞后,这可能是模型中Green-Ampt公式假设入渗过程中湿润锋后部的土壤是完全饱和状态,使得模型计算结果产生一定的滞后作用;同时在产流前期径流增大阶段模型计算值上升趋势较实际值缓慢。进一步分析可以看出,30,45mm/h下,产流后期的模拟值均大于实测值;而在60,75,90mm/h下,稳定产流阶段的实测值与计算值匹配度较好,这可能是雨强较小时,雨滴在未产流时段对地表的溅蚀增大了地面粗糙度,且在产流时段内坡面水深较浅,径流流经过程中极易形成洼地微地形截流,使得径流动力学特征的作用不显著,从而使得模拟结果与实测结果误差较大,这也说明本文建立的模型更适合于模拟>60mm/h下的产流过程。

表3不同雨强下的拟合入渗率参数c、R2、RMSE和NSE导出到EXCEL

对不同雨强下的径流实测过程进一步分析可以看出,随着雨强的增大,开始产流时间均开始前移,5种雨强下的实际开始产流时间分别为21.5,8.5,4.8,3.6,2.5min,90mm/h较30mm/h提前19min;且径流量在产流开始后的一段时间内呈上升趋势,雨强越大,上升速率越快,而后均保持为波动缓慢增加趋势,这也与模型模拟过程相一致;进一步分析可以看出,90mm/h下产流稳定阶段的单宽径流量较其他雨强下分别依次增大6.3,2.7,1.6,1.2倍,说明随着雨强的变大,其对径流的贡献作用也随之增强。

图3不同雨强下的单宽径流实测值与计算值

3.2泥沙输移过程模拟

在上述径流过程模拟的基础上,结合实测资料对不同雨强下径流产沙过程进行模拟,从表4可以看出,参数a和b随雨强的增大整体上均呈增大趋势,相关系数R2均大于0.76,且分别分布于0.13～0.15和0.10～0.12,RMSE和NSE的取值变化均与上述径流过程较为一致,其取值分别分布于2.07～3.37和0.635～0.819,说明本文建立的产沙模型能够很好地模拟不同雨强下的产沙过程;进一步分析可以看出随着雨强的增大,雨滴动能的溅蚀打击和径流冲刷剥蚀的作用均逐渐增强,这与Yang等[27]的研究结论较为一致,也说明雨强是影响水土流失的主要因素之一,大雨强赋予雨滴更强的动能以及超渗而形成的大径流,使得表层土壤结构更易崩解随径流迁移,从而加剧表层土壤流失。

表4不同雨强下的拟合参数a、b、R2、RMSE和NSE导出到EXCEL

为了进一步直观描述本文建立的模型与实测值之间的拟合关系,将其分别绘于图4。可以看出模型计算值在一定地程度上均可以很好地描述产流过程中径流携沙量随时间推移的变化趋势,与前述径流模拟过程相同的是30,45mm/h2种雨强下实测泥沙数据均与模拟曲线的匹配性较差,这可能的原因是雨强较小时,雨强和径流对泥沙的动力学搬运过程产生较大滞后作用,并且泥沙运输过程本身便存在概率因素,单一的考虑动力学作用也可能是造成模拟结果较差的原因,同时加之径流过程模拟的误差传递也会影响模拟计算结果。而随着雨强的增大,在60,75,90mm/h下,模型计算结果与实测结果的拟合度则好于其他雨强。

进一步对实测产沙过程进行分析比较,可以看出随着降雨强度的增大,径流产沙起始时刻开始逐渐提前;产沙率在径流产沙初始时刻均较小,而后迅速增大至某一水平后放慢增大速度,且稳定产沙时期的产沙率随雨强的增大而呈增大趋势;雨强为90mm/h的情况下的产沙率较其他雨强依次分别增大3.45,2.2,1.56,1.25倍,说明雨强的增大可以显著增加径流携沙量,进而加剧了坡耕地农田表层土壤的流失。

图4不同雨强下的产沙率实测值与计算值

3.3参数敏感性分析

为了进一步从量化角度来分析径流和泥沙数据对模型中参数的敏感性,在60,90mm/h下,改变不同参数以及将对应变化参数情况下模拟计算值和实测值的均方根误差值随时间的变化过程绘于图5和图6,用Dif(计算值与实测值差值平方的均值)来判断模型预测的准确性,其计算方法见公式34。

由图5可知,在参数c从0变化到0.45时,2种雨强下的Dif值分别从0.63,1.03cm2/min增大到3.17,5.79cm2/min,说明雨强越大,参数c的大幅变化对Dif的影响越大;进一步分析可以看出,当参数c<0.05时,Dif值的变化是非常微小的,这也说明地表径流对参数c的敏感程度很小,换言之,当雨强<90mm/h时,参数c取小于0.05的值对模型模拟结果影响较小。

图5不同雨强下随参数c变化的径流观测值与预测值的均方根误差值

图6为2种雨强下变量a和b的值对Dif的影响;图6(a)和图6(b)表明在变量b为相应的定值时(表3),a值的变化对2种雨强下Dif的影响。当a从0.04变化到0.45时,Dif值均呈先减小后增大的趋势,且在0.14左右的时候取得最小值;图6(c)和图6(d)表明在变量a为不同雨强下所对应的定值时(表3),b值的变化对2种雨强下Dif的影响。当b从0.01变化到0.42的时候,2种雨强下Dif值的变化在前期均呈减小趋势,且在b>0.11以后开始呈现出差异,60mm/h雨强下的上升趋势较为缓慢,而90mm/h雨强下的上升趋势则变得相对明显;进一步对比图6(a)～图6(d),不难看出在相同雨强下,参数a和b变化对Dif值的影响均表现为较一致的规律,即Dif值对参数a的变化较参数b更为敏感。

图6不同雨强下随参数a、b变化的泥沙观测值与预测值的均方根误差值

4、结论

(1)本文建立的近似解析模型可以较好地模拟不同雨强下的产流产沙过程,且当雨强>60mm/h时,其模拟的结果优于30,45mm/h;对实测数据进行分析表明,开始产流时间与雨强呈负相关关系,雨强越大,径流贡献率和产沙率也随之增大;曼宁糙率系数均分布在0.00026～0.00028min/cm1/3,而入渗率参数c随雨强的增大呈减小趋势,可以用指数函数c=9×10-4p-0.644对其进行描述,相关系数R2为0.9765;产沙模型中的参数a和b均分别分布于0.13～0.15和0.10～0.12,由于泥沙运动的偶然性因素,产沙过程的模型模拟精度低于径流过程。

(2)通过对模型参数进行敏感性分析表明,入渗率参数c取值为0～0.05时,Dif值变化较小,即地表径流对c不敏感,当c>0.05时,Dif值开始逐渐增大;径流产沙过程中,Dif值随参数a和b的变化均有不同程度的变化。总而言之,径流产沙率对参数a比参数b更敏感。

参考文献:

[1]郑粉莉,徐锡蒙,覃超.沟蚀过程研究进展[J].农业机械学报,2016,47(8):48-59,116.

[7]曹文洪,张启舜,姜乃森.黄土地区一次暴雨产沙数学模型的研究[J].泥沙研究,1994(1):1-13.

[8]贾媛媛,郑粉莉,杨勤科.黄土高原小流域分布式水蚀预报模型[J].水利学报,2005,36(3):328-332.

[9]张玉斌,郑粉莉,贾媛媛.WEPP模型概述[J].水土保持研究,2004,11(4):146-149.

[10]李文杰,王洁,张帅帅.一维坡面土壤侵蚀数学模型研究[J].重庆交通大学学报(自然科学版),2015,34(1):72-77.

[11]戚隆溪,黄兴法.坡面降雨径流和土壤侵蚀的数值模拟[J].力学学报,1997,29(3):343-348.

[13]陶汪海,吴军虎.坡地产流产沙规律数值模拟研究[J].水土保持学报,2016,30(1):54-57.

[14]邹宏荣,王烨,陆国宾,等.淮河上游产沙模拟研究[J].水资源保护,2016,32(4):80-83,99.

[16]范文涛,牛文全,张振华,等.降水头垂直入渗Green-Ampt模型显式近似解研究[J].灌溉排水学报,2012,31(1):50-53.

[17]刘继龙,马孝义,张振华.不同条件下Green-Ampt模型累积入渗量显函数的适用性[J].应用基础与工程科学学报,2010,18(1):11-19.

[18]王全九,来剑斌,李毅.Green-Ampt模型与Philip入渗模型的对比分析[J].农业工程学报,2002,18(2):13-16.

[19]张振华,潘英华,蔡焕杰,等.Green-Ampt模型入渗率显式近似解研究[J].农业系统科学与综合研究,2006,22(4):308-311.

[29]吴军虎,邵凡凡,刘侠.蚯蚓粪对土壤团聚体组成和入渗过程水分运移的影响[J].水土保持学报,2019,33(3):81-87.

邵凡凡,吴军虎.黄土区空闲坡耕地地表径流产沙近似解析模型[J].水土保持学报,2020,34(03):120-128.

基金:国家自然科学基金重点项目(51239009);西北旱区生态水利工程国家重点实验室科研项目(2016ZZKT-9);陕西省教育厅重点实验室科研计划项目(17JS096).