黃土丘陵溝壑區(qū)坡溝系統(tǒng)切溝侵蝕數(shù)值模擬

李斌兵，肖培青，余叔同

(1.武警工程大學(xué)信息技術(shù)教研室，710086，西安;2.黃河水利科學(xué)研究院水利部黃河泥沙重點實驗室，450003，鄭州;3.西北農(nóng)林科技大學(xué)黃土高原土壤侵蝕與旱地農(nóng)業(yè)國家重點實驗室，712100，陜西楊凌)

切溝侵蝕包括溝底下切、溝壁擴(kuò)展、溝頭溯源侵蝕及其相伴的輸沙沉積等過程，侵蝕過程既受水力作用又受重力作用，既有必然性又有隨機(jī)性，行為非常復(fù)雜，研究難度較大，相關(guān)的研究成果較少，切溝重力侵蝕試驗和野外觀測資料缺乏，對認(rèn)識切溝侵蝕力學(xué)機(jī)制并對其侵蝕量進(jìn)行預(yù)報非常困難。溝蝕侵蝕過程及預(yù)報的問題嚴(yán)重影響了流域侵蝕產(chǎn)沙預(yù)報的結(jié)果，成為制約其發(fā)展的“瓶頸”。在切溝侵蝕演變過程與侵蝕產(chǎn)沙研究中，由于自然裸露坡面或陡坡農(nóng)地細(xì)溝侵蝕向淺溝和切溝演變是一個相對漫長的過程;因此很難觀測到切溝形成和演化過程中某位置和某時刻的瞬時水力學(xué)參數(shù)和土壤剝離及沉積量，目前只能通過概化模型進(jìn)行模擬。在切溝侵蝕水動力學(xué)機(jī)制研究中，國內(nèi)外學(xué)者［1-8］對坡面侵蝕的動力學(xué)規(guī)律進(jìn)行了大量的研究，目前國外幾個代表性的模型在徑流水力學(xué)基礎(chǔ)方面還存在較大的差異。例如，WEPP 模型采用了徑流剪切力，EUROSEM和LISEM 模型采用了單位水流功率，而GUEST 模型則采用了水流功率。剪切力、水流功率和單位水流功率間存在明顯的差異，哪一個參數(shù)更能準(zhǔn)確描述土壤侵蝕過程，仍需要進(jìn)一步深入研究。國內(nèi)坡面流水力特性試驗研究多集中于細(xì)溝流研究，對于坡溝系統(tǒng)徑流水力學(xué)特性的研究較少。王文龍等［4，7］采用室內(nèi)人工模擬降雨試驗的方法，初步探討了坡溝系統(tǒng)各種侵蝕方式的水力特性及其開始發(fā)生侵蝕的動力臨界條件。肖培青等［9］采用變坡度坡溝系統(tǒng)概化模型和人工模擬降雨試驗相結(jié)合的方法，研究了坡溝系統(tǒng)水流水力學(xué)參數(shù)變化特征。但是目前切溝流水動力學(xué)過程研究的深度和廣度還遠(yuǎn)不能滿足黃土高原水土流失規(guī)律研究和數(shù)學(xué)模型建設(shè)的需要，大大制約了具有物理成因的水動力學(xué)模型的進(jìn)展。在切溝形態(tài)研究中，許多學(xué)者［10-14］提出了包含野外動態(tài)監(jiān)測、高精度GPS 測量、三維激光掃描和GIS 技術(shù)相結(jié)合的溝蝕演變過程研究方法，估算了切溝侵蝕的體積和侵蝕量，但沒有給出基于物理過程的切溝侵蝕模型。

綜上所述，國內(nèi)外切溝侵蝕模型研究還處于起步階段，面臨著很多挑戰(zhàn)，現(xiàn)有的模型還不能模擬局部土壤剝離和泥沙沉積量變化，因而不能反映切溝瞬時形態(tài)的特征與時空變化。筆者通過對切溝侵蝕過程和機(jī)制的實驗與數(shù)值模擬研究，揭示在水力和重力復(fù)合作用下，切溝侵蝕下切和輸沙沉積過程機(jī)制，建立切溝下切方程，進(jìn)行數(shù)值模擬，在坡溝系統(tǒng)上進(jìn)行了驗證，以揭示切溝侵蝕過程機(jī)制，建立切溝侵蝕模型，推動侵蝕預(yù)報模型發(fā)展。

1 坡溝系統(tǒng)切溝侵蝕模擬

1.1 坡溝系統(tǒng)切溝侵蝕模型

在切溝由起初發(fā)育到進(jìn)入穩(wěn)定狀態(tài)過程中，起始階段發(fā)展變化最快，切溝的長度、深度、寬度、表面積和體積變化劇烈，遠(yuǎn)離穩(wěn)定，隨著切溝侵蝕過程的進(jìn)行，切溝發(fā)育演變趨緩，直至穩(wěn)定。其發(fā)展過程為:如果水流剪切力大于土壤臨界剪切力，在表土層或溝的底部形成一個矩形溝道，隨著降雨過程的持續(xù)進(jìn)行，切溝溝頭不斷溯源侵蝕使其加長、溝槽不斷下切使其溝槽加深、溝壁崩塌使其加寬，切溝逐步由切溝發(fā)育的早期向中期、后期演變，最終進(jìn)入切溝發(fā)展的穩(wěn)定期。切溝是具有地理空間的微地貌，在切溝發(fā)育演變過程中，切溝的各種特征不僅隨時間變化，而且在其上、中、下游等不同的部位也有空間上的差異;因此，切溝發(fā)育過程是一個時空動態(tài)的演化過程，切溝下切和溝壁坍塌等發(fā)育的速度與水流速度、水流剪切力、水深、流量、水流含沙量、泥沙沉速、土壤質(zhì)地與土壤可蝕性等諸多水沙及土壤參數(shù)相關(guān)。

質(zhì)量守恒原理是水沙運動遵循的普遍規(guī)律。在相距為dx 的坡段，在dt 時距內(nèi)輸入和輸出該坡段的沙量之差與該坡段內(nèi)水流中攜帶的沙量變化量之和，必然等于dt 時距dx 距離內(nèi)水流下伏土壤表面的凈侵蝕或沉積量;因而，可基于質(zhì)量平衡原理，建立切溝侵蝕輸沙微分方程，描述不同時間點沿程侵蝕輸沙量的時空動態(tài)變化，方程見式(1)。式(1)的左邊定義了dx 坡段dt 時距內(nèi)輸入輸出沙量變化量與水流攜帶沙量變化量之和，方程的右邊則定義了該坡段該時距內(nèi)的凈侵蝕或沉積。

侵蝕輸沙過程的推移，會影響沿溝長剝離量與沉積量的變化，因此，在侵蝕過程中，坡面局部的土壤剝離和泥沙沉積會引起局部坡度的瞬態(tài)改變。其原因是溝床隨時間是非均勻侵蝕，在時段dt 內(nèi)地表剝離或沉積的凈泥沙量可用式(2)表示，即為切溝侵蝕下切方程，定義了溝底高程隨著分離的土壤顆粒和沉積的泥沙量變化而變化的過程。

坡溝實驗結(jié)果［9］表明，切溝底部土壤分離速率主要與水流有效剪切力成正比，見式(3)。

式(1)和(2)的初始條件是:

式中:A 為徑流橫截面面積，m2;Q 為清水流量，m3/s;x 為沿坡距離，m;t 為降雨時間，s;C 為平均的體積比含沙量，m3/m3;Sd為切溝底部土壤顆粒分離率，m/s;Z 為溝底的高度，m;w 為徑流寬度，m;vf為泥沙顆粒在湍流下的沉降速度，m/s;ε 為土壤孔隙度，m3/m3;K 為土壤可蝕性系數(shù)，s/m;τ 為溝底剪切力，m2/s2;τcr為臨界剪切力，m2/s2;h 為徑流的水深，m;ρ 為水的密度，kg/m3;g 為重力加速度，m/s2;S 為水力坡度，m/m。

對于湍流中的泥沙沉降速度，一般小于在靜止水中的速度，選擇張瑞瑾［15］提出的泥沙沉降速度計算公式

式中:vf為泥沙沉降速度，m/s;ν 為黏滯系數(shù)，cm2/s;ds為泥沙粒徑，mm;ρs為泥沙密度，kg/m3，通常(ρs－ρ)/ρ=1.65。

1.2 模型求解

根據(jù)產(chǎn)匯流方程［16-17］和切溝侵蝕模型，利用坡溝系統(tǒng)的DEM，按照時空變化進(jìn)行推演計算，時間推演以次降雨開始時間為起點，降雨結(jié)束為終點，將時間劃分為若干時刻，空間從坡溝徑流源點開始沿每條徑流線一直演算到坡溝出口，在每一時刻上遍歷每一條徑流的全部網(wǎng)格，求解出每一網(wǎng)格的流量Q(i，j)和泥沙含量C(i，j)。

式(1)中，C 為待求的，C 隨時間t 和沿坡距離x 而變化，下式中用C(x，t) 表示C 的動態(tài)變化，用差分方法對其進(jìn)行求解，首先用Δx 與Δt 分割(0，0)－(x，t)平面，見圖1，x 為沿溝長的距離(m)，t 為降雨結(jié)束時間(s)，采用下列近似公式:

式中θ 為Pressmann 隱式差分格式中的權(quán)重系數(shù)。

將偏微分方程(1)離散為下列差分方程(6)，進(jìn)而求解。

圖1 徑流-時間離散圖Fig.1 Runoff-time discrete graph

式(2)求解過程為:

把方程(3)代入方程(2)得:

設(shè)a=Kgρh

因S 為坡度，即高程隨位置的變化率，所以

式(7)可以用Lax-Wendroff 顯式方程求解，即用有限差分代替方程中的微分，再利用迭代算法求解。其推導(dǎo)過程為:

2 模擬結(jié)果與檢驗

文獻(xiàn)［9］中給出了在安塞水土保持綜合試驗站實體模型上進(jìn)行的模擬降雨實驗結(jié)果，實體模型是根據(jù)黃土丘陵溝壑區(qū)坡溝空間分異特征和野外調(diào)查量化的典型坡溝系統(tǒng)的參數(shù)特征制作的，模型土槽高6.01 m、寬3.05 m、投影長12 m，投影面積40 m2，見圖2。待降雨結(jié)束后，快速用LIDAR 三維激光掃描儀獲取地面地形信息，將測量數(shù)據(jù)導(dǎo)入ArcGIS 平臺進(jìn)行處理，經(jīng)過插值生成模擬降雨試驗后的DEM，將切溝初始形成時刻(溝底開始出現(xiàn)下切)生成的DEM 作為切溝模型的起算DEM 數(shù)據(jù)，網(wǎng)格大小為20 cm，采用徑向基插值方法生成DEM。降雨強(qiáng)度分別為0.83、1.66 和3.60 mm/min，土壤吸力參數(shù)sf取572.52 mm，土壤初始飽和含水率為0.155，飽和含水率為0.465，曼寧糙率系數(shù)n=0.108(s/m1/3)，飽和導(dǎo)水率Ks=0.010 07 mm/s，切溝出口處的徑流量和侵蝕量觀測計算方法見文獻(xiàn)［16-17］。

圖2 坡溝系統(tǒng)實體模型示意圖Fig.2 Schematic diagram of solid model of hillslope-gully system

不同降雨強(qiáng)度降雨試驗中經(jīng)測量得到的溝深度值與本文模型計算得到的模擬值見表1、2、3。0.83 mm/min 降雨強(qiáng)度降雨試驗中溝深模擬值與實測值的相對誤差基本上在±25%以內(nèi)，1.66 mm/min 降雨強(qiáng)度降雨試驗中，溝深模擬值與實測值相對誤差基本上在±15%以內(nèi)，3.60 mm/min 降雨強(qiáng)度降雨試驗中，溝深模擬值與實測值相對誤差基本上在±16%以內(nèi)。大、中、小3 種降雨強(qiáng)度下的檢驗結(jié)果說明在不同降降雨強(qiáng)度度下采用本文模型計算出的溝深模擬值與實測值較接近。進(jìn)一步分析溝深相對誤差的分布可知坡面上部溝深模擬值相對誤差較大，而在坡面下部相對誤差較小，這是因為坡面上部溝蝕發(fā)育沒有坡面下部劇烈，坡面上部雛形侵蝕溝的溝深較小，給測量精度帶來影響，測量數(shù)據(jù)可能與實際值有所偏差。總的來說，采用本文切溝模型模擬能夠反映溝蝕溝深的發(fā)育情況。

表1 0.83 mm/min 降雨強(qiáng)度溝深模擬相對誤差Tab.1 Relative errors of simulated gully depth at 0.83 mm/min rainfall intensity

續(xù)表1

表2 1.66 mm/min 降雨強(qiáng)度下溝深模擬的相對誤差Tab.2 Relative errors of simulated gully depth at 1.66 mm/min rainfall intensity

表3 3.60 mm/min 降雨強(qiáng)度下溝深模擬的相對誤差Tab.3 Relative errors of simulated gully depth at 3.60 mm/min rainfall intensity

3 結(jié)論與討論

1) 依據(jù)質(zhì)量守恒原理，考慮了影響切溝水沙運動的各物理要素，建立了坡溝系統(tǒng)切溝侵蝕模型，給出了有限差分法的推導(dǎo)過程。將建立的切溝侵蝕模型與降雨入滲、產(chǎn)匯流模型進(jìn)行耦合計算，利用迭代算法編寫程序求解了不同時刻和不同坡長下的切溝深度數(shù)據(jù)。經(jīng)過與試驗觀測結(jié)果比較檢驗，平均相對誤差在15%以內(nèi)，取得了較高的模擬精度。

2) 從模型計算與試驗觀測的比較結(jié)果看，坡面上部溝深模擬值相對誤差較大，而在坡面下部相對誤差較小。這是因為坡面下部相比上部，切溝侵蝕劇烈。總的來說，本文模型模擬能夠反映切溝溝深的發(fā)育情況。

3) 盡管在坡溝系統(tǒng)下，動態(tài)切溝模型取得了較好的模擬效果，下一步擬開展在流域尺度下的驗證和測試。其難度在于:切溝的空間尺度大，分布位置上地形多險峻，在我國黃土丘陵溝壑區(qū)地形則更是破碎，切溝侵蝕過程很復(fù)雜;因此，時至今日，國內(nèi)外還沒有布設(shè)過野外切溝侵蝕過程的水沙觀測。因此，難以為模型驗證提供必需的水流水力學(xué)參數(shù)和切溝形態(tài)演化數(shù)據(jù)。實體模型仍是目前有效的手段之一，研究具有嚴(yán)格比尺意義的實體模型，滿足幾何相似、降雨相似、水力相似等條件是建立小流域侵蝕模型的基礎(chǔ)和重要工作。

［1］ Nearing M A，Norton L D，Bulgakav G A，et al.Hydraulics and erosion in eroding rills［J］.Water Resources Research，1997，33(4):865-876

［2］ Foster G R，Huggins L F，Meyer L D.A laboratory study of rill hydraulics:I.velocity relationship［J］.Trans of ASAE，1984，27(3):790-796

［3］ 張科利，唐克麗.黃土坡面細(xì)溝侵蝕能力的水動力學(xué)試驗研究［J］.土壤學(xué)報，2000，37(1):9-15

［4］ 王文龍，雷阿林，李占斌.土壤侵蝕鏈內(nèi)細(xì)溝淺溝切溝動力機(jī)制研究［J］.水科學(xué)進(jìn)展，2003，14(4):471-475

［5］ 丁文峰，李占斌，魯克新，等.坡面細(xì)溝發(fā)生臨界水動力條件初探［J］.土壤學(xué)報，2003，40(6):822-828

［6］ 張光輝，劉寶元，何小武.黃土區(qū)原狀土壤分離過程的水動力學(xué)機(jī)理研究［J］.水土保持學(xué)報，2005，19(4):48-52

［7］ 王文龍，王兆印，雷阿林，等.黃土丘陵溝壑區(qū)坡溝系統(tǒng)不同侵蝕方式的水力特性初步研究［J］.中國水土保持科學(xué)，2007，5(2):11-17

［8］ 李占斌，秦百順，亢偉，等.陡坡面發(fā)育的細(xì)溝水動力學(xué)特性室內(nèi)試驗研究［J］.農(nóng)業(yè)工程學(xué)報，2008，24(6):64-68

［9］ 肖培青，鄭粉莉，姚文藝.坡溝系統(tǒng)坡面徑流流態(tài)及水力學(xué)參數(shù)特征研究［J］.水科學(xué)進(jìn)展，2009，20(2):236-240

［10］伍永秋，劉寶元.切溝、切溝侵蝕與預(yù)報［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報，2000，8(2):134-142.

［11］胡剛，伍永秋，劉寶元，等.GPS 和GIS 進(jìn)行短期溝蝕研究初探:以東北漫川漫崗黑土區(qū)為例［J］.水土保持學(xué)報，2004，18(4):16-19.

［12］游智敏，伍永秋，劉寶元.利用GPS 進(jìn)行切溝侵蝕監(jiān)測研究［J］.水土保持學(xué)報，2004，18(5):91-94

［13］胡剛，伍永秋，劉寶元，等.東北漫崗黑土區(qū)切溝侵蝕發(fā)育特征［J］.地理學(xué)報，2007，62(11):1165-1173

［14］張鵬，鄭粉莉，王彬，等.高精度GPS，三維激光掃描和測針板三種測量技術(shù)監(jiān)測溝蝕過程的對比研究［J］.水土保持通報，2008，28(5):11-15

［15］張瑞瑾.河流泥沙動力學(xué)［M］.北京:中國水利水電出版社，1998:48-55

［16］李斌兵，鄭粉莉，王占禮.黃土丘陵區(qū)小流域分布式水文和侵蝕模型建立和模擬［J］.土壤通報，2010，41(5):1153-1160

［17］李斌兵，鄭粉莉.黃土坡面不同土地利用下的降雨入滲模擬與數(shù)值計算［J］.干旱地區(qū)農(nóng)業(yè)研究，2008，26(5):118-123