李斌兵，肖培青，余叔同

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

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

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

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

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

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

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

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

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

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

式中:A 為徑流橫截面面積，m2;Q 為清水流量，m3/s;x 為沿坡距離，m;t 為降雨時(shí)間，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。

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

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

1.2 模型求解

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

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

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

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

圖1 徑流-時(shí)間離散圖Fig.1 Runoff-time discrete graph

式(2)求解過(guò)程為:

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

設(shè)a=Kgρh

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

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

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

文獻(xiàn)［9］中給出了在安塞水土保持綜合試驗(yàn)站實(shí)體模型上進(jìn)行的模擬降雨實(shí)驗(yàn)結(jié)果，實(shí)體模型是根據(jù)黃土丘陵溝壑區(qū)坡溝空間分異特征和野外調(diào)查量化的典型坡溝系統(tǒng)的參數(shù)特征制作的，模型土槽高6.01 m、寬3.05 m、投影長(zhǎng)12 m，投影面積40 m2，見(jiàn)圖2。待降雨結(jié)束后，快速用LIDAR 三維激光掃描儀獲取地面地形信息，將測(cè)量數(shù)據(jù)導(dǎo)入ArcGIS 平臺(tái)進(jìn)行處理，經(jīng)過(guò)插值生成模擬降雨試驗(yàn)后的DEM，將切溝初始形成時(shí)刻(溝底開(kāi)始出現(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，切溝出口處的徑流量和侵蝕量觀測(cè)計(jì)算方法見(jiàn)文獻(xiàn)［16-17］。

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

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

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

續(xù)表1

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

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

3 結(jié)論與討論

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

2) 從模型計(jì)算與試驗(yàn)觀測(cè)的比較結(jié)果看，坡面上部溝深模擬值相對(duì)誤差較大，而在坡面下部相對(duì)誤差較小。這是因?yàn)槠旅嫦虏肯啾壬喜?，切溝侵蝕劇烈。總的來(lái)說(shuō)，本文模型模擬能夠反映切溝溝深的發(fā)育情況。

3) 盡管在坡溝系統(tǒng)下，動(dòng)態(tài)切溝模型取得了較好的模擬效果，下一步擬開(kāi)展在流域尺度下的驗(yàn)證和測(cè)試。其難度在于:切溝的空間尺度大，分布位置上地形多險(xiǎn)峻，在我國(guó)黃土丘陵溝壑區(qū)地形則更是破碎，切溝侵蝕過(guò)程很復(fù)雜;因此，時(shí)至今日，國(guó)內(nèi)外還沒(méi)有布設(shè)過(guò)野外切溝侵蝕過(guò)程的水沙觀測(cè)。因此，難以為模型驗(yàn)證提供必需的水流水力學(xué)參數(shù)和切溝形態(tài)演化數(shù)據(jù)。實(shí)體模型仍是目前有效的手段之一，研究具有嚴(yán)格比尺意義的實(shí)體模型，滿足幾何相似、降雨相似、水力相似等條件是建立小流域侵蝕模型的基礎(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ì)溝侵蝕能力的水動(dòng)力學(xué)試驗(yàn)研究［J］.土壤學(xué)報(bào)，2000，37(1):9-15

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

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

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

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

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

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

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

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

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

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

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

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

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

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