一维坡面土壤侵蚀数学模型研究

李文杰，王 洁，张帅帅

(1. 重庆交通大学 国家内河航道整治工程技术研究中心，重庆 400074；2. 重庆交通大学 省部共建水利水运工程教育部重点实验室，重庆 400074)



一维坡面土壤侵蚀数学模型研究

李文杰1，王 洁2，张帅帅2

(1. 重庆交通大学 国家内河航道整治工程技术研究中心，重庆 400074；2. 重庆交通大学 省部共建水利水运工程教育部重点实验室，重庆 400074)

依据水沙运动过程的力学机理，将坡面水文模型和侵蚀产沙模型嵌套，构建了一维坡面侵蚀产沙数学模型，模型的水文模块涵盖了植被截留、蒸散发、入渗、地下水和地表水运动过程，侵蚀产沙模块则模拟了雨滴溅蚀、坡面流侵蚀以及坡面径流的输沙过程。利用径流小区径流和泥沙实测资料对模型进行检验，结果表明：水文模块的模拟精度较高，侵蚀产沙模块也能合理地计算坡面产沙量，整体效果较好。

水利工程；坡面；水文过程；侵蚀产沙过程；一维数学模型；径流小区模拟

防止坡地冲蚀和水土流失与社会经济发展息息相关[1]，因此对水土流失过程模拟意义重大。对流域水土流失进行模拟时，其下垫面及气候条件一般呈非均匀性，常用的处理方法是将流域离散为坡面和沟道。降雨降落至坡面并汇集，从而侵蚀坡面上的泥沙并携带其进入沟道，因此坡面是构成流域的核心，也是流域水土流失的根源。坡面一方面为流域的产沙提供可输移侵蚀土壤颗粒，另一方面为流域的泥沙输移提供水力条件，故坡面侵蚀产沙过程模拟是流域模拟的基础，也是水土流失治理的关键。

坡面侵蚀产沙数学模型可以分为经验模型和机理模型，如国外的Usle、Rusle模型等，国内的江忠善模型[2]、范瑞瑜模型[3]、金争平模型[4]、李钜章模型[5]、张继生模型[6]等。总体来看，国内经验型模型居多、自主开发的通用模型较少，而引进国外模型会存在一些地域适应性或经济合理性等的限制。

笔者基于水沙运动的物理过程构建具有物理机制的一维坡面土壤侵蚀模型，先建立水文模块，在此基础上开发侵蚀产沙模块，包括雨滴溅蚀、坡面流侵蚀以及坡面径流的输沙过程，然后将模型应用于遂宁水土保持试验站径流小区，对模型的两个模块进行检验，以期为坡面土壤流失的治理及流域模型的开发奠定基础。

1 坡面水文模型

首先构建坡面水文模型，水文模型模拟包括降雨、植被截留、蒸散发、入渗、地表径流、壤中流和地下水过程。模型对此进行了如下概化：

1)降雨首先被植被冠层截留，满足冠层截留能力后的雨量继续下降至地表。植被冠层截留采用覆被相关法计算[7]，是植被叶面积指数的函数，相关数据可根据NDVI获得。

2)蒸散发的模拟先采用Penman公式[8]计算潜在蒸发能力，实际的蒸散发量则由3个部分组成：植被冠层的截留雨量首先开始蒸发，当其不能满足潜在蒸发能力时，植被根系所在土壤层的水量从植被叶面发生蒸腾；若地表无植被覆盖，蒸发则从地表土壤开始。

3)经截留和蒸发之后到达地表的降雨，如果降雨强度小于土壤入渗能力，则降雨全部入渗至土壤非饱和带，在满足土壤饱和度之后产生壤中流或继续下渗补充地下水。

地下水层与非饱和带的分界面通常处于不断波动的状态，数值计算中很难作为边界条件确定下来，故模型统一考虑地表以下水流的垂向运动，将潜水面包含在计算区域内，计算时将两个方向的水流运动解耦，首先确定初始地下水深度以及非饱和带含水量分布；然后根据非饱和带中含水量的分布以及地下水位分别计算沿坡面方向运动的壤中流和地下水流量；最后统一计算壤中流和地下水出流后地表以下水流的垂向运动，得到非饱和带含水量以及地下水位的重分布，作为下一时层计算壤中流和地下水流量的初始条件。饱和-非饱和土壤水垂向运动用以水头Ψ为因变量的一维Richards方程描述[9]，采用隐式差分格式离散，得到三对角方程组采用追赶法求解。

4)壤中流和地下水简化为平行坡面的一维流动，模型中地下水只作为河道的基流补给，用达西定律和质量守恒描述；而水分在非饱和带中可能会越过大部分土壤的体积沿着优先途径流动产生优先流，或者局部土层饱和后发生侧向流动而产生壤中流。

5)到达地表的降雨经过入渗及填洼后，多余雨量则形成坡面径流，采用一维运动波模型[10]描述：

(1)

式中：x为沿坡面向下方向；t为时间，s；h为水深，m；q为单宽流量，m2/s；p为净降雨强度，m/s；θ为坡面倾角，(°)；i为入渗率，m/s；n为Manning糙率系数；S0为坡面坡度，S0=sinθ。

采用Preissmann四点偏心格式离散运动波模型，用Newton迭代法求解离散后的方程组。其中，净降雨强度p为扣除植被冠层截流、蒸散发之后的净降雨强度。

2 坡面侵蚀产沙模型

在坡面水文模型的基础上，增加侵蚀产沙模块，将其概化为雨滴溅蚀、坡面流侵蚀以及泥沙输移过程，降雨击溅地表形成雨滴溅蚀，坡面流冲刷地表产生坡面流侵蚀，雨滴溅蚀和坡面流侵蚀的泥沙发生沉积或者通过坡面径流向坡脚搬运。

2.1 雨滴溅蚀

雨滴溅蚀的原动力是雨滴的冲击力，雨滴降落至地表，由于其冲击力的击打，使土粒以跃移的形式搬运，产生溅蚀，故溅蚀与雨滴的物理性质有关。此外，雨滴击打在下垫面上才能发生溅蚀，所以下垫面的土壤性质、坡度、坡面流水深、地表覆被率也是重要影响因素。综合考虑以上因素，笔者采用了李文杰模型[11]，将溅蚀产生的主要矛盾归结为降雨和土壤，分别用雨滴动能和土壤分离指数表示，坡度对溅蚀的影响通过引入无量纲的坡度影响因子来反映且包含了临界坡度的影响，坡面流和植被覆盖对溅蚀的屏蔽效应则通过定义无量纲的屏蔽系数来确定，雨滴溅蚀模型如式(2)：

Di=FW·k·KE·f(J)·e-0.01Cv

(2)

式中：Di为雨滴溅蚀率，g/(m2·h)；FW为坡面流水深屏蔽系数；k为土壤分离指数，g/J；KE为雨滴动能，J/(m2·h)；f(J)为坡度影响因子；Cv为植被覆盖度，%。

2.2 坡面流侵蚀

坡面流侵蚀是降雨在坡面形成的水流对坡面土壤的分散和输移，模型认为坡面流造成的土壤侵蚀率与坡面流的挟沙力和含沙量的差值成正比：

(3)

式中：Df为坡面水流的土壤侵蚀率，kg/(m2·s)；Cs为水流中的泥沙浓度，kg/m3；β为水流侵蚀效率系数；Tc为坡面流挟沙能力，kg/m3；ωs为泥沙颗粒沉降速度，m/s。

由于坡面流的水深一般较小，极易受到降雨击溅的扰动，从而影响坡面水流的输沙能力。Li Wenjie，等[12]通过量纲分析，得到了一个基于水流功率的坡面流挟沙力公式，降雨的影响通过临界水流功率来反映，无量纲的泥沙通量计算如式(4)：

(4)

求得泥沙通量后则可得到坡面流的挟沙力。

2.3 坡面产输沙一维模型

雨滴溅蚀和坡面流侵蚀作为泥沙输移的源项，根据质量守恒原理，坡面侵蚀产沙的一维连续方程为：

(5)

式中：符号意义同前。

采用Preissmann四点偏心差分格式离散求解，离散后的方程为：

(6)

令M=2Δt/Δx，则系数A1，A2，A3，A4分别为：

(7)

给定初始条件和边界条件均为0求解。

3 模型应用

3.1 径流小区物理参数

径流小区位于遂宁水土保持试验站，地处嘉陵江中下游重点水土流失区的遂宁市安居镇，东经105°28′51″，北纬30°21′51″，海拔300 m。径流小区共有5°，10°，15°，20°，25°这5个相邻的露天坡面，下端的观测房通过管道与坡面相连，用量水池和三角堰进行径流和泥沙的观测。试验站设有一个雨量计，认为其测得的降雨可以覆盖径流小区。径流小区坡面为矩形，水平投影长度9.52 m，坡面的长度、宽度、坡度以及土壤深度等物理参数详见表1。

表1 径流小区参数

3.2 土壤属性及降雨侵蚀数据

在小区现场采集土样，测得泥沙颗粒级配曲线见图1。依据国际制土壤质地分类三角图，该土壤属于砂质土壤。侵蚀产沙模块未对粒径进行分组，采用中值粒径计算，其值为0.008 mm。

图1 径流小区泥沙颗粒级配曲线

收集了1994，1999，2005，2007以及2008年的8场降雨过程，详见表2。其中2008年两场降雨测有径流过程，其余只有径流系数。对应的泥沙数据中，只测有降雨的侵蚀总量。

表 2 降雨信息统计

3.3 模拟结果及分析

3.3.1 径流模拟结果分析

收集的8场降雨中，2008年8月9日和2008年8月24—25日两次降雨测有相应的流量过程，模型选择这两场降雨径流过程率定，采用的时间步长为60 s，用人工试错法调整径流参数使得模拟结果与实测结果吻合，见图2。

图2 率定期流量过程模拟结果

2008年两次降雨-径流过程的率定结果表明，模拟流量过程与实测流量过程基本吻合，尤其对降雨的主峰值以及峰值时间模拟准确，次峰值偏差稍大，整体模拟效果令人满意。保持率定后的参数不变，以1994，1999，2005和2007年的6次降雨作为验证期，考虑到没有流量过程的资料，采用径流系数来验证，验证结果和误差分布分别见图3和图4。

图3 径流系数验证结果

图4 径流系数验证结果相对误差

验证结果表明，径流系数的模拟值与实测值基本吻合，个别坡面的模拟偏差稍大，定义(│计算值-实测值│/实测值)为相对误差，可得径流系数模拟的相对误差主要分布在1%～10%，部分位于10%～20%之间，个别误差较大。经计算，模拟值与实测值的相关系数和确定性系数分别为0.93和0.87，认为坡面侵蚀产沙模型的水文模块精度可靠，可进行下一步侵蚀产沙模块的模拟。

3.3.2 泥沙模拟结果

保持率定的径流相关参数不变，进行泥沙模块的率定和验证。对应每场降雨的泥沙资料无含沙量过程，只有侵蚀总量，因此根据侵蚀量进行泥沙参数的率定。同样选用2008年8月9日和2008年8月24—25日作为泥沙率定期，计算步长采用60 s，调整泥沙模块相关参数使模拟侵蚀量与实测侵蚀量吻合，率定结果见图5。

图5 泥沙率定结果

验证期为1994，1999，2005和2007年的6场降雨，侵蚀量模拟结果见图6和图7。

图6 侵蚀量验证结果

图7 侵蚀量验证结果相对误差

验证结果表明，泥沙模块对侵蚀量的模拟结果与实测结果总体符合较好。侵蚀量模拟的相对误差主要分布在30%以下及附近。个别误差较大甚至超过50%，原因可能是侵蚀产沙模块的本身的精度不够，同时也是径流模拟误差的放大。经计算，模拟值与实测值的相关系数和确定性系数分别为0.97和0.92，认为泥沙模块总体比较可靠。但模型对多峰过程模拟结果误差较大，这种结果受多种因素的影响，比如对水文计算引入的误差、地表以下水流运动的一维化简化等。例如，土壤在一次降雨峰值之后达到饱和，从而在二次峰值来临时我们将不再考虑该因素。而事实上，土壤在两次降雨峰值之间可能会与其他形式的水分流动有一定的交换，并在二次峰值时未处于饱和状态，因此可能导致较大的误差；又或者由于人类活动的影响，对地表形态造成了较大的破坏，使得模型对其的模拟产生了一定的偏差。由此可知，该模型仍有待进一步改进。

4 结 语

基于地表水沙运动的动力学机理，在坡面水文模型的基础上构建侵蚀产沙模块，通过模拟植被截留、蒸散发、入渗、地表径流、壤中流、地下水、雨滴溅蚀、坡面流侵蚀以及泥沙输移等物理过程，建立了一维坡面侵蚀产沙动力学模型。水文模块将地表以下水流的垂向运动和坡面方向运动解耦以简化计算，采用四点偏心格式离散求解运动波方程得到地表径流模型。泥沙侵蚀模块将雨滴溅蚀和坡面流侵蚀作为泥沙输移的源项，输移过程用质量守恒方程描述。用遂宁水土保持试验站径流小区的实测数据对模型的合理性以及精确度进行了检验，结果表明模型结构合理，整体模拟效果令人满意。

本坡面侵蚀产沙模型的研究为流域模型的开发奠定了基础，也为水土流失治理提供了理论依据和技术支持。模型的不足之处在于对泥沙侵蚀模块的概化，应继续加强对水文泥沙过程机理的深入研究并模拟不同形式的侵蚀产沙方式，以提高模型产沙计算的精度。

[1] 王平义.森林与洪水及土壤冲蚀相关性研究[J].重庆交通学院学报,2003,22(6):121-124. Wang Pingyi.Investigation on the co-relations between forest and flood/soil erosion [J].Journal of Chongqing Jiaotong University,2003,22(6):121-124.

[2] 江忠善,王志强,刘志.黄土丘陵区小流域土壤侵蚀空间变化定量研究[J].土壤侵蚀与水土保持学报,1996,1(2):1-9. Jiang Zhongshan,Wang Zhiqiang,Liu Zhi.Quantitative study on spatial variation of soil erosion in a small watershed in the loess hilly region [J].Journal of Soil Erosion and Soil and Water Conservation,1996,1(2):1-9.

[3] 范瑞瑜.黄河中游地区小流域土壤流失量计算方程的研究[J].中国水土保持,1985,2:12-18. Fan Ruiyu.Study on formula for computation of soil loss from small watershed on the middle reaches of the Huanghe River [J].Soil and Water Conservation in China,1985,2:12-18.

[4] 金争平,赵焕勋，何泰，等.皇甫川区小流域土壤侵蚀量预报方程研究[J].水土保持学报,1991,5(1):8-18. Jin Zhengping,Zhao Huanxun,He Tai,et al.An approach to the equation predicting soil erosion of small watershed in Huangfuchuan area [J].Journal of Soil and Water Conservation,1991,5(1):8-18.

[5] 李钜章,景可,李风新.黄土高原多沙粗沙区侵蚀模型探讨[J].地理科学进展,1999,18(1):46-53. Li Juzhang,Jing Ke,Li Fengxin.A study on the erosion model in areas with high and coarse sediment yield [J].Progress in Geography,1999,18(1):46-53.

[6] 张继生,王平义,刘亚辉.三峡库区中小流域产沙数学模型[J].重庆交通学院学报,2004,23(4):117-120. Zhang Jisheng,Wang Pingyi,Liu Yahui.Mathematical model on sediment yield of middle small watershed at the area of the Three Goreges Reservior [J].Journal of Chongqing Jiaotong University,2004,23(4):117-120.

[7] 仪垂祥,刘开瑜.植被截留降水量公式的建立[J].土壤侵蚀与水土保持学报,1996,2(3):47-49. Yi Chuixiang,Liu Kaiyu.Research on a formula of rainfull interception by vegetation [J].Journal of Soil Erosion and Soil and Water Conservation,1996,2(3):47-49.

[8] 王加虎.分布式水文模型理论与方法研究[D].南京:河海大学,2006. Wang Jiahu.The Theory and Method of Distributed Hydrologic Model [D].Nanjing:Hohai University,2006.

[9] G R Foster,L J Lane,J D Nowlin,et al.Estimating erosion and sediment yield on field-sized areas [J].Transactions of the ASAE,1981,24(5):1253-1262.

[10] 刘青泉,李家春,陈力，等.坡面流及土壤侵蚀动力学(I)—坡面流[J].力学进展,2004,34(3):360-372. Liu Qingquan,Li Jiachun,Chen Li,et al.Overland flow and soil erosion dynamics (I):Overland flow [J].Advances in Mechanics,2004,34(3):360-372.

[11] 李文杰,李丹勋,王兴奎.雨滴溅蚀模型研究[J].应用基础与工程科学学报,2011,19(5):689-699. Li Wenjie,Li Danxun,Wang Xingkui.Study on raindrop splash detachment model [J].Journal of Basic Science and Engineering,2011,19(5):689-699.

[12] Li Wenjie,Li Danxun,Wang Xingkui.An approach to estimating sediment transport capacity of overland flow [J].Science China:Technological Sciences,2011,54(10):2649-2656.

Study on One-Dimensional Mathematical Model of Slope Erosion

Li Wenjie1, Wang Jie2, Zhang Shuaishuai2

(1. National Engineering Research Center for Inland Waterway Regulation, Chongqing Jiaotong University, Chongqing 400074, China;2. Key Laboratory of Ministry of Education for Hydraulic & Water Transport Engineering, Chongqing Jiaotong University, Chongqing 400074, China)

A one dimensional model of slope erosion and sediment yield was established by combining the hydrological module and soil erosion and sediment yield module, according to the dynamical mechanism of hydrological and soil erosion processes. The processes of interception, evapotranspiration, infiltration, groundwater and overland flow were simulated in the hydrology module, and the sediment module simulated the raindrop splash detachment, overland flow erosion and sediment transport processes. The whole model was verified by the measured data of runoff plot and sediment. The simulation results indicate that the hydrological module is relatively precise and the soil erosion and sediment yield module simulates the sediment yield reasonably, so the overall effect is better.

hydraulic engineering; slope; hydrological process; soil erosion and sediment yield process; one dimensional mathematical model; simulation of runoff plot

10.3969/j.issn.1674-0696.2015.01.16

2013-03-25；

2013-06-21

重庆市教委科学技术研究项目(KJ120406)；重庆市科委自然科学基金项目(cstc2013jcyjA30003);交通运输部应用基础研究项目(2014329814310)

李文杰(1984—)，男，河北秦皇岛人，副教授，博士，主要从事泥沙运动力学、流域水沙过程模拟及航道整治等方面的研究。E-mail:li_wj1984@163.com。

TV143.4

A

1674-0696(2015)01-072-06