李文杰，王 潔，張帥帥

(1. 重慶交通大學(xué) 國家內(nèi)河航道整治工程技術(shù)研究中心，重慶 400074；2. 重慶交通大學(xué) 省部共建水利水運(yùn)工程教育部重點(diǎn)實驗室，重慶 400074)

?

一維坡面土壤侵蝕數(shù)學(xué)模型研究

李文杰1，王 潔2，張帥帥2

(1. 重慶交通大學(xué) 國家內(nèi)河航道整治工程技術(shù)研究中心，重慶 400074；2. 重慶交通大學(xué) 省部共建水利水運(yùn)工程教育部重點(diǎn)實驗室，重慶 400074)

依據(jù)水沙運(yùn)動過程的力學(xué)機(jī)理，將坡面水文模型和侵蝕產(chǎn)沙模型嵌套，構(gòu)建了一維坡面侵蝕產(chǎn)沙數(shù)學(xué)模型，模型的水文模塊涵蓋了植被截留、蒸散發(fā)、入滲、地下水和地表水運(yùn)動過程，侵蝕產(chǎn)沙模塊則模擬了雨滴濺蝕、坡面流侵蝕以及坡面徑流的輸沙過程。利用徑流小區(qū)徑流和泥沙實測資料對模型進(jìn)行檢驗，結(jié)果表明：水文模塊的模擬精度較高，侵蝕產(chǎn)沙模塊也能合理地計算坡面產(chǎn)沙量，整體效果較好。

水利工程；坡面；水文過程；侵蝕產(chǎn)沙過程；一維數(shù)學(xué)模型；徑流小區(qū)模擬

防止坡地沖蝕和水土流失與社會經(jīng)濟(jì)發(fā)展息息相關(guān)[1]，因此對水土流失過程模擬意義重大。對流域水土流失進(jìn)行模擬時，其下墊面及氣候條件一般呈非均勻性，常用的處理方法是將流域離散為坡面和溝道。降雨降落至坡面并匯集，從而侵蝕坡面上的泥沙并攜帶其進(jìn)入溝道，因此坡面是構(gòu)成流域的核心，也是流域水土流失的根源。坡面一方面為流域的產(chǎn)沙提供可輸移侵蝕土壤顆粒，另一方面為流域的泥沙輸移提供水力條件，故坡面侵蝕產(chǎn)沙過程模擬是流域模擬的基礎(chǔ)，也是水土流失治理的關(guān)鍵。

坡面侵蝕產(chǎn)沙數(shù)學(xué)模型可以分為經(jīng)驗?zāi)Ｐ秃蜋C(jī)理模型，如國外的Usle、Rusle模型等，國內(nèi)的江忠善模型[2]、范瑞瑜模型[3]、金爭平模型[4]、李鉅章模型[5]、張繼生模型[6]等?？傮w來看，國內(nèi)經(jīng)驗型模型居多、自主開發(fā)的通用模型較少，而引進(jìn)國外模型會存在一些地域適應(yīng)性或經(jīng)濟(jì)合理性等的限制。

筆者基于水沙運(yùn)動的物理過程構(gòu)建具有物理機(jī)制的一維坡面土壤侵蝕模型，先建立水文模塊，在此基礎(chǔ)上開發(fā)侵蝕產(chǎn)沙模塊，包括雨滴濺蝕、坡面流侵蝕以及坡面徑流的輸沙過程，然后將模型應(yīng)用于遂寧水土保持試驗站徑流小區(qū)，對模型的兩個模塊進(jìn)行檢驗，以期為坡面土壤流失的治理及流域模型的開發(fā)奠定基礎(chǔ)。

1 坡面水文模型

首先構(gòu)建坡面水文模型，水文模型模擬包括降雨、植被截留、蒸散發(fā)、入滲、地表徑流、壤中流和地下水過程。模型對此進(jìn)行了如下概化：

1)降雨首先被植被冠層截留，滿足冠層截留能力后的雨量繼續(xù)下降至地表。植被冠層截留采用覆被相關(guān)法計算[7]，是植被葉面積指數(shù)的函數(shù)，相關(guān)數(shù)據(jù)可根據(jù)NDVI獲得。

2)蒸散發(fā)的模擬先采用Penman公式[8]計算潛在蒸發(fā)能力，實際的蒸散發(fā)量則由3個部分組成：植被冠層的截留雨量首先開始蒸發(fā)，當(dāng)其不能滿足潛在蒸發(fā)能力時，植被根系所在土壤層的水量從植被葉面發(fā)生蒸騰；若地表無植被覆蓋，蒸發(fā)則從地表土壤開始。

3)經(jīng)截留和蒸發(fā)之后到達(dá)地表的降雨，如果降雨強(qiáng)度小于土壤入滲能力，則降雨全部入滲至土壤非飽和帶，在滿足土壤飽和度之后產(chǎn)生壤中流或繼續(xù)下滲補(bǔ)充地下水。

地下水層與非飽和帶的分界面通常處于不斷波動的狀態(tài)，數(shù)值計算中很難作為邊界條件確定下來，故模型統(tǒng)一考慮地表以下水流的垂向運(yùn)動，將潛水面包含在計算區(qū)域內(nèi)，計算時將兩個方向的水流運(yùn)動解耦，首先確定初始地下水深度以及非飽和帶含水量分布；然后根據(jù)非飽和帶中含水量的分布以及地下水位分別計算沿坡面方向運(yùn)動的壤中流和地下水流量；最后統(tǒng)一計算壤中流和地下水出流后地表以下水流的垂向運(yùn)動，得到非飽和帶含水量以及地下水位的重分布，作為下一時層計算壤中流和地下水流量的初始條件。飽和-非飽和土壤水垂向運(yùn)動用以水頭Ψ為因變量的一維Richards方程描述[9]，采用隱式差分格式離散，得到三對角方程組采用追趕法求解。

4)壤中流和地下水簡化為平行坡面的一維流動，模型中地下水只作為河道的基流補(bǔ)給，用達(dá)西定律和質(zhì)量守恒描述；而水分在非飽和帶中可能會越過大部分土壤的體積沿著優(yōu)先途徑流動產(chǎn)生優(yōu)先流，或者局部土層飽和后發(fā)生側(cè)向流動而產(chǎn)生壤中流。

5)到達(dá)地表的降雨經(jīng)過入滲及填洼后，多余雨量則形成坡面徑流，采用一維運(yùn)動波模型[10]描述：

(1)

式中：x為沿坡面向下方向；t為時間，s；h為水深，m；q為單寬流量，m2/s；p為凈降雨強(qiáng)度，m/s；θ為坡面傾角，(°)；i為入滲率，m/s；n為Manning糙率系數(shù)；S0為坡面坡度，S0=sinθ。

采用Preissmann四點(diǎn)偏心格式離散運(yùn)動波模型，用Newton迭代法求解離散后的方程組。其中，凈降雨強(qiáng)度p為扣除植被冠層截流、蒸散發(fā)之后的凈降雨強(qiáng)度。

2 坡面侵蝕產(chǎn)沙模型

在坡面水文模型的基礎(chǔ)上，增加侵蝕產(chǎn)沙模塊，將其概化為雨滴濺蝕、坡面流侵蝕以及泥沙輸移過程，降雨擊濺地表形成雨滴濺蝕，坡面流沖刷地表產(chǎn)生坡面流侵蝕，雨滴濺蝕和坡面流侵蝕的泥沙發(fā)生沉積或者通過坡面徑流向坡腳搬運(yùn)。

2.1 雨滴濺蝕

雨滴濺蝕的原動力是雨滴的沖擊力，雨滴降落至地表，由于其沖擊力的擊打，使土粒以躍移的形式搬運(yùn)，產(chǎn)生濺蝕，故濺蝕與雨滴的物理性質(zhì)有關(guān)。此外，雨滴擊打在下墊面上才能發(fā)生濺蝕，所以下墊面的土壤性質(zhì)、坡度、坡面流水深、地表覆被率也是重要影響因素。綜合考慮以上因素，筆者采用了李文杰模型[11]，將濺蝕產(chǎn)生的主要矛盾歸結(jié)為降雨和土壤，分別用雨滴動能和土壤分離指數(shù)表示，坡度對濺蝕的影響通過引入無量綱的坡度影響因子來反映且包含了臨界坡度的影響，坡面流和植被覆蓋對濺蝕的屏蔽效應(yīng)則通過定義無量綱的屏蔽系數(shù)來確定，雨滴濺蝕模型如式(2)：

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

(2)

式中：Di為雨滴濺蝕率，g/(m2·h)；FW為坡面流水深屏蔽系數(shù)；k為土壤分離指數(shù)，g/J；KE為雨滴動能，J/(m2·h)；f(J)為坡度影響因子；Cv為植被覆蓋度，%。

2.2 坡面流侵蝕

坡面流侵蝕是降雨在坡面形成的水流對坡面土壤的分散和輸移，模型認(rèn)為坡面流造成的土壤侵蝕率與坡面流的挾沙力和含沙量的差值成正比：

(3)

式中：Df為坡面水流的土壤侵蝕率，kg/(m2·s)；Cs為水流中的泥沙濃度，kg/m3；β為水流侵蝕效率系數(shù)；Tc為坡面流挾沙能力，kg/m3；ωs為泥沙顆粒沉降速度，m/s。

由于坡面流的水深一般較小，極易受到降雨擊濺的擾動，從而影響坡面水流的輸沙能力。Li Wenjie，等[12]通過量綱分析，得到了一個基于水流功率的坡面流挾沙力公式，降雨的影響通過臨界水流功率來反映，無量綱的泥沙通量計算如式(4)：

(4)

求得泥沙通量后則可得到坡面流的挾沙力。

2.3 坡面產(chǎn)輸沙一維模型

雨滴濺蝕和坡面流侵蝕作為泥沙輸移的源項，根據(jù)質(zhì)量守恒原理，坡面侵蝕產(chǎn)沙的一維連續(xù)方程為：

(5)

式中：符號意義同前。

采用Preissmann四點(diǎn)偏心差分格式離散求解，離散后的方程為：

(6)

令M=2Δt/Δx，則系數(shù)A1，A2，A3，A4分別為：

(7)

給定初始條件和邊界條件均為0求解。

3 模型應(yīng)用

3.1 徑流小區(qū)物理參數(shù)

徑流小區(qū)位于遂寧水土保持試驗站，地處嘉陵江中下游重點(diǎn)水土流失區(qū)的遂寧市安居鎮(zhèn)，東經(jīng)105°28′51″，北緯30°21′51″，海拔300 m。徑流小區(qū)共有5°，10°，15°，20°，25°這5個相鄰的露天坡面，下端的觀測房通過管道與坡面相連，用量水池和三角堰進(jìn)行徑流和泥沙的觀測。試驗站設(shè)有一個雨量計，認(rèn)為其測得的降雨可以覆蓋徑流小區(qū)。徑流小區(qū)坡面為矩形，水平投影長度9.52 m，坡面的長度、寬度、坡度以及土壤深度等物理參數(shù)詳見表1。

表1 徑流小區(qū)參數(shù)

3.2 土壤屬性及降雨侵蝕數(shù)據(jù)

在小區(qū)現(xiàn)場采集土樣，測得泥沙顆粒級配曲線見圖1。依據(jù)國際制土壤質(zhì)地分類三角圖，該土壤屬于砂質(zhì)土壤。侵蝕產(chǎn)沙模塊未對粒徑進(jìn)行分組，采用中值粒徑計算，其值為0.008 mm。

圖1 徑流小區(qū)泥沙顆粒級配曲線

收集了1994，1999，2005，2007以及2008年的8場降雨過程，詳見表2。其中2008年兩場降雨測有徑流過程，其余只有徑流系數(shù)。對應(yīng)的泥沙數(shù)據(jù)中，只測有降雨的侵蝕總量。

表 2 降雨信息統(tǒng)計

3.3 模擬結(jié)果及分析

3.3.1 徑流模擬結(jié)果分析

收集的8場降雨中，2008年8月9日和2008年8月24—25日兩次降雨測有相應(yīng)的流量過程，模型選擇這兩場降雨徑流過程率定，采用的時間步長為60 s，用人工試錯法調(diào)整徑流參數(shù)使得模擬結(jié)果與實測結(jié)果吻合，見圖2。

圖2 率定期流量過程模擬結(jié)果

2008年兩次降雨-徑流過程的率定結(jié)果表明，模擬流量過程與實測流量過程基本吻合，尤其對降雨的主峰值以及峰值時間模擬準(zhǔn)確，次峰值偏差稍大，整體模擬效果令人滿意。保持率定后的參數(shù)不變，以1994，1999，2005和2007年的6次降雨作為驗證期，考慮到?jīng)]有流量過程的資料，采用徑流系數(shù)來驗證，驗證結(jié)果和誤差分布分別見圖3和圖4。

圖3 徑流系數(shù)驗證結(jié)果

圖4 徑流系數(shù)驗證結(jié)果相對誤差

驗證結(jié)果表明，徑流系數(shù)的模擬值與實測值基本吻合，個別坡面的模擬偏差稍大，定義(│計算值-實測值│/實測值)為相對誤差，可得徑流系數(shù)模擬的相對誤差主要分布在1%～10%，部分位于10%～20%之間，個別誤差較大。經(jīng)計算，模擬值與實測值的相關(guān)系數(shù)和確定性系數(shù)分別為0.93和0.87，認(rèn)為坡面侵蝕產(chǎn)沙模型的水文模塊精度可靠，可進(jìn)行下一步侵蝕產(chǎn)沙模塊的模擬。

3.3.2 泥沙模擬結(jié)果

保持率定的徑流相關(guān)參數(shù)不變，進(jìn)行泥沙模塊的率定和驗證。對應(yīng)每場降雨的泥沙資料無含沙量過程，只有侵蝕總量，因此根據(jù)侵蝕量進(jìn)行泥沙參數(shù)的率定。同樣選用2008年8月9日和2008年8月24—25日作為泥沙率定期，計算步長采用60 s，調(diào)整泥沙模塊相關(guān)參數(shù)使模擬侵蝕量與實測侵蝕量吻合，率定結(jié)果見圖5。

圖5 泥沙率定結(jié)果

驗證期為1994，1999，2005和2007年的6場降雨，侵蝕量模擬結(jié)果見圖6和圖7。

圖6 侵蝕量驗證結(jié)果

圖7 侵蝕量驗證結(jié)果相對誤差

驗證結(jié)果表明，泥沙模塊對侵蝕量的模擬結(jié)果與實測結(jié)果總體符合較好。侵蝕量模擬的相對誤差主要分布在30%以下及附近。個別誤差較大甚至超過50%，原因可能是侵蝕產(chǎn)沙模塊的本身的精度不夠，同時也是徑流模擬誤差的放大。經(jīng)計算，模擬值與實測值的相關(guān)系數(shù)和確定性系數(shù)分別為0.97和0.92，認(rèn)為泥沙模塊總體比較可靠。但模型對多峰過程模擬結(jié)果誤差較大，這種結(jié)果受多種因素的影響，比如對水文計算引入的誤差、地表以下水流運(yùn)動的一維化簡化等。例如，土壤在一次降雨峰值之后達(dá)到飽和，從而在二次峰值來臨時我們將不再考慮該因素。而事實上，土壤在兩次降雨峰值之間可能會與其他形式的水分流動有一定的交換，并在二次峰值時未處于飽和狀態(tài)，因此可能導(dǎo)致較大的誤差；又或者由于人類活動的影響，對地表形態(tài)造成了較大的破壞，使得模型對其的模擬產(chǎn)生了一定的偏差。由此可知，該模型仍有待進(jìn)一步改進(jìn)。

4 結(jié) 語

基于地表水沙運(yùn)動的動力學(xué)機(jī)理，在坡面水文模型的基礎(chǔ)上構(gòu)建侵蝕產(chǎn)沙模塊，通過模擬植被截留、蒸散發(fā)、入滲、地表徑流、壤中流、地下水、雨滴濺蝕、坡面流侵蝕以及泥沙輸移等物理過程，建立了一維坡面侵蝕產(chǎn)沙動力學(xué)模型。水文模塊將地表以下水流的垂向運(yùn)動和坡面方向運(yùn)動解耦以簡化計算，采用四點(diǎn)偏心格式離散求解運(yùn)動波方程得到地表徑流模型。泥沙侵蝕模塊將雨滴濺蝕和坡面流侵蝕作為泥沙輸移的源項，輸移過程用質(zhì)量守恒方程描述。用遂寧水土保持試驗站徑流小區(qū)的實測數(shù)據(jù)對模型的合理性以及精確度進(jìn)行了檢驗，結(jié)果表明模型結(jié)構(gòu)合理，整體模擬效果令人滿意。

本坡面侵蝕產(chǎn)沙模型的研究為流域模型的開發(fā)奠定了基礎(chǔ)，也為水土流失治理提供了理論依據(jù)和技術(shù)支持。模型的不足之處在于對泥沙侵蝕模塊的概化，應(yīng)繼續(xù)加強(qiáng)對水文泥沙過程機(jī)理的深入研究并模擬不同形式的侵蝕產(chǎn)沙方式，以提高模型產(chǎn)沙計算的精度。

[1] 王平義.森林與洪水及土壤沖蝕相關(guān)性研究[J].重慶交通學(xué)院學(xué)報,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] 江忠善,王志強(qiáng),劉志.黃土丘陵區(qū)小流域土壤侵蝕空間變化定量研究[J].土壤侵蝕與水土保持學(xué)報,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] 范瑞瑜.黃河中游地區(qū)小流域土壤流失量計算方程的研究[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] 金爭平,趙煥勛，何泰，等.皇甫川區(qū)小流域土壤侵蝕量預(yù)報方程研究[J].水土保持學(xué)報,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] 李鉅章,景可,李風(fēng)新.黃土高原多沙粗沙區(qū)侵蝕模型探討[J].地理科學(xué)進(jìn)展,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] 張繼生,王平義,劉亞輝.三峽庫區(qū)中小流域產(chǎn)沙數(shù)學(xué)模型[J].重慶交通學(xué)院學(xué)報,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].土壤侵蝕與水土保持學(xué)報,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].南京:河海大學(xué),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] 劉青泉,李家春,陳力，等.坡面流及土壤侵蝕動力學(xué)(I)—坡面流[J].力學(xué)進(jìn)展,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].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報,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

重慶市教委科學(xué)技術(shù)研究項目(KJ120406)；重慶市科委自然科學(xué)基金項目(cstc2013jcyjA30003);交通運(yùn)輸部應(yīng)用基礎(chǔ)研究項目(2014329814310)

李文杰(1984—)，男，河北秦皇島人，副教授，博士，主要從事泥沙運(yùn)動力學(xué)、流域水沙過程模擬及航道整治等方面的研究。E-mail:li_wj1984@163.com。

TV143.4

A

1674-0696(2015)01-072-06