胡 鵬，余明輝

（武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢 430072）

1 研究背景

坡面流是指降雨或融雪在扣除地面截留、填洼與下滲等損失后在重力作用下沿坡面表層流動的一種特殊水流，有時也包括雨水在坡面上部下滲后，經(jīng)過表層土壤，以壤中流的形式在坡面下部復而流出地面，再度形成的一種片狀水流。坡面流水力特性取決于許多因素，如降水強度和歷時、土壤質(zhì)地或種類、前期水分條件、植被密度和類型、地貌特性包括洼坑和小丘數(shù)量和大小、坡度和坡長、邊界穩(wěn)定性條件［1］。作為水流污染物遷移、土壤侵蝕及泥沙輸移的主要動力因素，不同邊界條件下坡面流的形成、發(fā)展過程是該領(lǐng)域的研究基礎(chǔ)。

隨著計算機技術(shù)和計算方法的迅速發(fā)展，坡面流的數(shù)值模擬求解越來越受到重視。作為坡面徑流模型中重要的參數(shù)，阻力系數(shù)是影響坡面流計算精度的關(guān)鍵參數(shù)［2-3］。對坡面流阻力特性研究主要是通過野外和實驗室內(nèi)的放水或降雨模擬試驗，根據(jù)實驗結(jié)果應用回歸分析等方法建立阻力系數(shù)與雷諾數(shù)、降雨強度、土壤粒徑、坡度等因子的表達式，如姚文藝［4］、Abrahams等［5］、Hirsch［6］、Hu等［7］和Lawrence［8］建立了形式各異的阻力系數(shù)公式。Zhang 等［9］將坡面上阻力的處理分為三種類型：阻力系數(shù)為常數(shù)；不同土地類型、地表覆蓋情況、耕作措施而采用不同阻力系數(shù)；考慮阻力系數(shù)的空間非均勻性，即其在每個網(wǎng)格的變化。由于坡面流阻力受土壤類型、植被覆蓋、床面粗糙度及降雨影響等眾多因素影響，難以從理論上描述，大多數(shù)水文模型中，仍將阻力系數(shù)作為常數(shù)處理。盡管阻力系數(shù)為常數(shù)使用較廣泛，但其選取具有較強的經(jīng)驗性，越來越多的研究關(guān)注阻力的空間非均勻性對坡面匯流的影響［10-11］。

綜上所述，現(xiàn)有成果盡管在坡面流阻力特性研究上取得了一些成果，但由于試驗條件的差異，簡化處理方法不同，建立的阻力公式仍有很大局限性，不便于應用。本文選取三種有代表性的阻力模式，即阻力系數(shù)為常數(shù)、基于阻力線性疊加原理的阻力分割模型、以淹沒度為變量的Lawrence 模型，對其在裸坡、有礫石覆蓋坡面、有植被覆蓋坡面等三種典型坡面的坡面徑流計算精度的影響規(guī)律展開研究，并對植被覆蓋坡面，探究降雨強度對不同阻力模式模擬結(jié)果間差異的影響。

2 研究方法和計算模式

為研究不同阻力模式在坡面上的適用性，本文建立了坡面降雨徑流數(shù)學模型，對形成坡面徑流的兩個主要過程降雨入滲和坡面匯流進行模擬。降雨入滲模型采用考慮坡度影響的Green-Ampt 模型，可完整模擬降雨入滲過程，得到入滲率隨時間變化曲線。坡面徑流模型采用擴散波模型，可模擬坡面匯流過程，得到坡面各水力要素的時空變化過程。模型的離散求解采用具有二階精度的Maccormark 有限差分格式，網(wǎng)格劃分與時間步長選取滿足柯朗數(shù)小于1 的條件。模型中選取了三種不同的阻力模式，對不同類型坡面進行了模擬，并將不同阻力模式得到的模擬值與實測值比較分析。

2.1 坡面降雨入滲模型Green-Ampt 模型是最早基于物理過程和毛管理論的入滲模型，運用Darcy定律得到入滲率和入滲量之間的關(guān)系及其隨時間的變化過程，能較好地考慮土壤飽和導水率、有效孔隙率、初始含水量和累計入滲量對入滲過程的影響。經(jīng)典Green-Ampt 模型僅適用于平坡入滲，由于坡面坡度較大，用其求解坡面入滲過程會引入一定的誤差。因此采用Chen 和Young 提出的考慮坡度對于土壤入滲的影響的改進模型［12］。

式中：i為入滲率，m/s；Ks為土壤飽和導水率，m/s；S為土壤吸力，m；I為累計入滲量，m；θs為土壤飽和含水率；θi為土壤初始含水率；γ為坡面傾角，°。

整個降雨入滲過程的入滲率可描述為：

式中：p為降雨強度，mm/h；tp為開始積水時間，tp=Ip p；Ip為累計入滲量臨界值，當i=pcosγ時，由下式推出Ip計算公式：

I為積水開始后的累計入滲量，也包括未積水時段的入滲量在內(nèi)。由于不是從t=0時刻開始積水，I的計算須采用修正后的公式：

式中：ts為假設由t=0開始積水到入滲率i=p時所需的時間，可以理解為一個虛擬時間，s，計算公式如下：

2.2 坡面徑流模型坡面徑流可由圣維南方程組進行模擬。但實際的坡面水流運動因邊界條件復雜，同時，由于坡面流水深很淺，在實際坡面流動中受微地貌影響很大，完整的圣維南方程組并不一定能夠很好地描述這種特殊的流動。簡化的圣維南方程即運動波與擴散波模型被廣泛的應用于坡面流模擬中，擴散波模型包括連續(xù)方程與動量方程，可表示如下［13］：

式中：h為水深，m；t為時間，s；p為降雨強度，m/s；i為入滲率，m/s；Sf為能坡；S0為坡面比降；x為沿坡面向下坐標；q為單寬流量，q=hu，m2/s；u為流速，m/s。

u可由Darcy-Weisbach公式計算：

式中：f為Darcy-Weisbach阻力系數(shù)；R為水力半徑。

2.3 阻力模式Darcy-Weisbach 阻力系數(shù)f是無量綱參數(shù)且有一定的理論基礎(chǔ)，大部分學者選擇f對坡面流阻力特性研究進行研究，蔣昌波等［14］和Abrahams 等［15］也認為f比Chezy 系數(shù)和Manning 糙率系數(shù)更適合描述坡面地表水流的阻力。本文選取了下面三種阻力模式：

阻力系數(shù)為常數(shù)：

Lawrence 模型：Lawrence 以淹沒度為變量，提出可將坡面流分為完全淹沒、臨界淹沒和部分淹沒三種狀態(tài)。完全淹沒時，基于紊流理論計算阻力；臨界淹沒時，假定摻長與特征粗糙度成正比關(guān)系；部分淹沒時，僅考慮粗糙單元的拖曳力，分別推導建立了三種狀態(tài)下阻力系數(shù)表達式［8］：

完全淹沒：

臨界淹沒：

部分淹沒：

式中：Λ為淹沒度，Λ=h Dr；h為水深；Dr為特征粗糙尺度；λ為床面粗糙單元的覆蓋率，%；CD為系數(shù)。

阻力分割模型：在坡面流阻力研究特性研究中，阻力線性疊加原理得到了廣泛應用，Abrahams提出根據(jù)地表特征差異將坡面流阻力分為4部分，降雨阻力、顆粒阻力、形態(tài)阻力和波阻力，通過阻力的線性疊加得到綜合阻力。根據(jù)Abrahams，Hirsch 的系列試驗結(jié)果，Hu 等人運用量綱分析與回歸分析得到了部分阻力計算公式［5，7］。降雨阻力公式采用Shen的試驗結(jié)果［16］：

降雨阻力：

顆粒阻力：

形態(tài)阻力：

波阻力：

綜合阻力：

式中：p為降雨強 度，mm/min；Re為 雷諾數(shù)，Re=uh υ，υ為運 動 黏度，m2/s；Fr為 弗 勞德數(shù)，g為重力加速度，m/s2。

3 計算結(jié)果與分析

裸坡、有礫石覆蓋坡面、有植被覆蓋坡面為三種常見坡面類型，用上述三種不同阻力模式的降雨徑流模型對各類型坡面進行模擬，將數(shù)值模擬結(jié)果與實測值進行比較分析。裸坡上試驗數(shù)據(jù)分別采用中國科學院水土保持研究所和Abban 分別進行的室內(nèi)人工降雨—入滲—徑流試驗數(shù)據(jù)，試驗中土樣分別為黃土和粉質(zhì)壤土［17-18］。礫石是指平均粒徑大于2 mm 小于64 mm 巖石或礦物碎屑物，礫石覆蓋坡面常見于沙漠地帶，有礫石覆蓋坡面采用Jomaa 進行的室內(nèi)實驗數(shù)據(jù)，礫石的覆蓋度為20%，粒徑5 ～7 cm［19］。植被有明顯的阻水滯流攔沙作用，因此植被過濾帶常被用于坡面防護。植被冠層對降雨有截留能力，因此需考慮植被蓄水能力，有植被覆蓋坡面采用Neibling進行的場地實驗，距坡腳2.5 m 范圍內(nèi)有植被覆蓋，植被覆蓋率為40%。上述試驗中土樣均采用當?shù)鼐哂写硇酝寥?，?jīng)測定其性質(zhì)后進行實驗，降雨均采用人工降雨模擬裝置，保持降雨強度均勻，降雨強度均取當?shù)亟涤甑牡湫痛碇?，坡面入口處無來流，坡面出口自由出流，能很好的反映天然實際情況。各試驗坡度均在緩坡范圍內(nèi)，對不同長度的坡長模擬對比結(jié)果分析也表明不同坡長下不同阻力模型對模擬結(jié)果的影響規(guī)律相一致，因此可認為試驗坡長具有一定的代表性。試驗過程中流量由尾部集流裝置定時取樣測量，流速測量采用染色劑法測得表面流速乘以修正系數(shù)即可得斷面平均流速。從整體上看，可認為各試驗具有一定的代表性與可比性。各坡面基本參數(shù)見表1。

表1 各坡面基本參數(shù)

圖1 裸坡出口實測流量與計算流量對比

圖1—3 為各類型坡面流模擬出口流量過程與實測值對比圖。由圖可知，在裸坡上，無論是黃土還是粉質(zhì)壤土，各阻力模式計算得到的流量過程并無差異，且均與實測資料吻合較好。在有礫石覆蓋坡面上與有植被覆蓋坡面上，漲水階段各模式模擬結(jié)果有較大差別，阻力不僅僅會因為邊界條件改變而變化，也會隨水深，流速等水力要素變化而改變，故阻力與水力要素間是互相影響的，且是一個復雜的動態(tài)過程，考慮了阻力時空變化的阻力分割模式計算結(jié)果與實測資料最吻合，由于Lawrence 模型在水深較小，即淹沒度較小的時候，由于其只考慮繞流阻力，會低估阻力系數(shù)，所以模擬流量值高于實測值。在坡面匯流的穩(wěn)定階段，可以發(fā)現(xiàn)各模式計算結(jié)果差別不大，此時，坡面流流量主要由降雨強度與入滲量決定，受阻力系數(shù)影響不大，這是由于達到穩(wěn)定階段后，坡面各水力要素隨時間變化很小，一個時間步長內(nèi)相差不大，因此由水力要素計算得到的阻力系數(shù)也幾乎不變，由坡面徑流模型中的連續(xù)性方程與動量方程可知，此時決定坡面流流量的關(guān)鍵主要是源項，即降雨量與入滲量。Huang 等［10］研究空間變化的阻力系數(shù)對流量過程影響時也得到了阻力系數(shù)對坡面徑流流量的影響在坡面匯流的漲水與退水階段較大，而在穩(wěn)定階段影響較小這一結(jié)論。對有植被覆蓋坡面，模擬了其退水階段，由圖3 可知，常數(shù)模型及阻力分割模型和實測值吻合較好，但由于Lawrence 模型只考慮淹沒度這個單一變量，在退水階段會低估退水時間。

圖2 礫石覆蓋坡面出口流量與計算流量對比

圖3 植被覆蓋坡面出口流量與計算流量對比

阻力模式的選取不僅影響坡面徑流流量過程，也會影響坡面徑流流速的模擬結(jié)果，坡面徑流流速與污染物遷移，泥沙侵蝕過程密切相關(guān)，因此有必要對不同阻力模式下的流速模擬結(jié)果進行討論。裸坡上模擬的出口流速過程與實測值對比見圖4，從總體上看，阻力分割模型模擬效果最好，在初始漲水階段，常數(shù)模型與實測值更為接近，而達到穩(wěn)定階段后，流速模擬值偏低，這是由于隨著水深的增大，阻力系數(shù)應變小，而常數(shù)模型并不能體現(xiàn)這一特性，這與Mügler 等［20］在研究坡面徑流模擬與污染物遷移時得出的阻力系數(shù)為常數(shù)的數(shù)學模型會低估最大流速結(jié)論一致。盡管在裸坡上Law?rence 模型模擬流量過程并無差別，但是流速模擬值較實測值總體偏高。

圖4 出口流速與實測值對比

圖5 流速沿程變化

圖5 為有植被覆蓋坡面在40 min 時坡面流速沿程分布。常數(shù)模型模擬流速沿程增大，最大流速在坡腳，阻力分割模型與Lawrence 模型模擬流速最大值并不是在坡腳，而是出現(xiàn)在植被覆蓋區(qū)域上部，很好地體現(xiàn)植被的滯流效應，這與楊春霞等［21］和孫佳美等［22］研究植被覆蓋對坡面產(chǎn)流產(chǎn)沙影響時得到的實驗結(jié)果相符合。

4 討論

坡面水文過程的發(fā)生主要動力來自于降雨，降雨強度會影響坡面的下滲過程，產(chǎn)流的時間，坡面徑流流量以及達到穩(wěn)定階段時間。對于同一坡面，隨著降雨強度的增大，其漲水階段時間會減少，徑流流量增大，從而導致不同阻力計算模式下坡面徑流流量模擬結(jié)果間差異減小，但當差異減少至一定程度趨于穩(wěn)定后，即使降雨強度增大，差異值也不再改變，將這一分界點定義為降雨強度閾值。對于不同坡面，降雨強度閾值會隨坡長、坡度以及植被覆蓋度變化而改變，因此有必要研究其變化規(guī)律。在有植被覆蓋的坡面，不同的阻力模式對坡面流流量模擬結(jié)果有較大影響，考慮時空變化的阻力系數(shù)的阻力分割模式與實測值吻合最好，為探究降雨強度對不同阻力模式模擬結(jié)果的影響，將阻力分割模式計算得到的結(jié)果作為標準值，通過分析不同降雨強度下阻力系數(shù)常數(shù)模型與阻力分割模型計算結(jié)果間的Nash-Sutcliffe 效率系數(shù)，確定降雨強度閾值。大于閾值時，不同阻力模式對計算結(jié)果影響甚微；反之，不同阻力模式計算差異較大，需選擇合適的阻力計算模式。

Nash-Sutcliffe效率系數(shù)計算公式為［23］：

圖6 Ef 與降雨強度關(guān)系

式中：Ef為Nash-Sutcliffe效率系數(shù)；Oj為測量值；Mj為模擬值；Oˉ為均值；n為時間點個數(shù)。

Ef值越接近1代表模擬值與實測值越吻合。

圖6（a）為固定坡長與植被覆蓋率條件下，坡度分別為7%、12%和25%時，根據(jù)模擬結(jié)果計算得到的Nash-Sutcliffe 效率系數(shù)與降雨強度間關(guān)系圖。圖中存在明顯轉(zhuǎn)折點，隨著降雨強度的增大，E f值不發(fā)生變化，且E f值接近1，說明模擬結(jié)果相差很小，轉(zhuǎn)折點的降雨強度即為降雨強度閾值。Pa?panicolaou 等通過研究坡面徑流功率與降雨強度關(guān)系，也發(fā)現(xiàn)坡面上存在臨界降雨強度，大于臨界降雨強度時，坡面阻力特性對坡面徑流模擬影響較小，幾乎可忽略不計［11］。在坡度為7%時，降雨強度閾值為119 mm/h，在坡度為12%時，降雨強度閾值為74 mm/h，在坡度為25%時，降雨強度閾值為67.5 mm/h。在坡長一定條件下，坡度越緩，其降雨強度閾值越大。比較同一降雨強度下不同坡度情況下的Nash-Sutcliffe 效率系數(shù)，也可發(fā)現(xiàn)坡度越大，Nash-Sutcliffe 效率系數(shù)也越高，即不同阻力模型對計算結(jié)果影響更小。

圖6（b）為固定坡長與坡度條件下，植被覆蓋度為20%和40%時，根據(jù)模擬結(jié)果計算得到的Nash-Sutcliffe 效率系數(shù)與降雨強度間關(guān)系圖。由圖可知，不同植被覆蓋度對降雨強度閾值影響不大，均在119 mm/h 左右。在同一降雨強度下，植被覆蓋度越低，Nash-Sutcliffe 效率系數(shù)略高，不同阻力模型對計算結(jié)果影響略小。

圖6（c）為固定坡度與植被覆蓋度條件下，坡長分別為3.05 m、6.1 m 和12.2 m 時，根據(jù)模擬結(jié)果計算得到的Nash-Sutcliffe 效率系數(shù)與降雨強度間關(guān)系圖。在坡長為3.05 m 時，其降雨強度閾值為67.5 mm/h；在坡長為6.1 m時，其降雨強度閾值為119 mm/h；在坡長為12.2 m時，其降雨強度閾值為164 mm/h。在坡度一定條件下，坡長越大，其降雨強度閾值越大。比較同一降雨強度下不同坡長情況下的Nash-Sutcliffe 效率系數(shù)，也可發(fā)現(xiàn)坡長越長，Nash-Sutcliffe 效率系數(shù)越小，即阻力模式的選取對計算結(jié)果影響越大。

5 結(jié)論

基于降雨徑流數(shù)學模型，研究了不同阻力模式對坡面匯流計算結(jié)果的影響，可得出以下結(jié)論：（1）在裸坡上，不同阻力模式坡面流流量模擬結(jié)果差別很小，流速模擬結(jié)果表明阻力分割模型與實際更吻合，阻力系數(shù)為常數(shù)模型會低估最大流速，而Lawrence 模型模擬流速值較實測值總體偏高；在有礫石覆蓋和植被覆蓋的坡面上，考慮了阻力系數(shù)時空變化的阻力分割模型模擬精度最高，且能體現(xiàn)植被的滯流效應。從整體上看，阻力系數(shù)對坡面徑流流量模擬的影響在坡面匯流的漲水與退水階段較大，而在穩(wěn)定階段很小。（2）對有植被覆蓋的坡面，存在降雨強度閾值，大于閾值時，不同阻力模式對計算結(jié)果影響甚微；反之，需選擇合適的阻力計算模式。并且坡度坡長是影響降雨強度閾值的關(guān)鍵因子，坡長越長，坡度越緩，其降雨強度閾值越大。