一維坡面土壤侵蝕數學模型研究

李文杰，王 潔，張帥帥

(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