黃土高原北部水蝕風蝕交錯區產流條件及徑流系數

盧龍彬，付 強，黃金柏

（東北農業大學 水利與建筑學院，哈爾濱150030）

降雨徑流在整個流域系統中是最為活躍的水文因子，是流域水循環的基本環節，研究降雨徑流過程有助于更好地開發利用流域水資源［1］。對于季節性水資源缺乏問題相對突出的黃土高原北部水蝕風蝕交錯區而言，如何實現對可利用水資源的準確評估已成為一個亟待解決的課題［2］。到目前為止，對流域產流的研究較多，如高軍俠等［3］利用人工模擬降雨試驗，根據水力學理論，統計分析了黃土高原坡面超滲徑流特征，指出黃土高原坡面產流以超滲產流為主；周宏飛等［4］認為，在下墊面條件一定時，降雨能否產生徑流，很大程度上取決于降雨歷時和降雨強度；李慧敏等［5］指出徑流量的大小在很大程度上取決于降雨因素?？偨Y降雨徑流數值模型的有關研究可知，迄今為止的分布式流域降雨—徑流過程的數值模型，大多是針對某一選定的流域或規模較小的試驗區，普遍存在通用性不強的特點，很難在不同地區中小流域尺度之間推廣應用［6］。針對黃土高原北部水蝕風蝕交錯區，地表水資源的研究相對較少［7－8］，鑒于此，本研究選取黃土高原北部六道溝流域為研究區，基于對觀測水文數據的綜合分析，揭示研究區的產流機制?；谶\動波理論結合GIS技術構建適用于研究流域的降雨—徑流數值模型，通過對降雨—徑流過程數值計算結果的分析，推求研究區的徑流系數。研究結果可為黃土高原北部中小尺度流域地表徑流的準確推求提供實用的計算方法，并期待為該地區地表水資源的深入研究以及流域數字化建設提供基礎數據。

1 材料與方法

1.1 研究區概況

六道溝流域位于陜西神木縣以西14 k m處（北緯35°20′—40°10′，東經103°33′—113°53′），地處黃土高原北部的水蝕風蝕交錯區。多年平均降水量為437 mm，降雨季節性分布很不均勻，70%以上的降雨集中在7—9月（圖1），且多暴雨，年潛在蒸發量超過1 000 mm［9－10］。流域內地形非常復雜，溝道縱橫交錯，土地沙化嚴重。對土層進行實際調查發現：該區域的表層土厚度為10～20 c m，表層以下是厚度超過15 m的堅硬的黃土層，再下層是侏羅紀砂巖［11－12］。

圖1 研究區多年降雨量（1956－2009年）與降雨量月分布

1.2 水文觀測

水文觀測主要包括降雨、河道內徑流、土壤水分等。因為六道溝流域內自然條件相似，為減少外界因素破壞，觀測點選在人為因素干擾較少的流域上游。降雨的觀測采用翻倒式自動計數雨量計（型號：7852 M-L10，尺寸 Φ165×240 H（mm）），每發生0.2 mm降雨記錄數據1次?？衫肂ox Car Pro.4.3軟件對觀測得到的降雨數據在時間軸上用不同的時間間隔如3，5，10 min等進行分割，從而可以得到不同時間步長序列的雨量數據。由于只有在雨季集中降雨發生期間才產生地表徑流，且流量受降雨強度的影響較大，所以采用精度較高的自動記數水位計（型號及制造商：KADEC21-MZPT （KONA System co.ltd）對地表徑流進行觀測，其觀測精度為誤差小于1 mm，觀測的時間步長為5 min；為了使流量數據和降雨數據在時間步長上保持一致，以時間步長5 min對降雨數據進行分割，即文中所提及的雨強（平均降雨強度）為5 min降雨量的算術平均值。土壤水分的觀測采用多探頭感應式土壤水分計進行（制造商：Decagon Devices Inc，型號：EC-5），測量的指標為土壤的體積含水率，測量精度為誤差小于3%。觀測點位于地表徑流觀測斷面兩側的坡面上（圖2）［13］。

圖2 地表徑流及土壤水分觀測示意圖

2 結果與分析

2.1 長歷時低強度降雨條件下的產流機制

圖3所示為2008-06-15—2008-06-16降雨、觀測斷面土壤水分變化及地表徑流產生的過程，本次降雨事件歷時1 335 min，降雨量為86.2 mm，最大降雨強度為6.7×10－6m／s（2 mm／5 min）。由圖3可知，降雨開始后5 c m處的土壤水分在某一時刻開始增加，而兩側坡面5 c m處的土壤水分含量開始增加的時間點不同；20 c m處的土壤水分含量開始增加的時刻滯后于5 c m處；50 c m處土壤水分含量在徑流產生之前保持相對穩定。產流前，由于降雨強度低于表層土壤的不飽和滲透系數，加之表層土壤未達到飽和狀態，所以不滿足產流條件。當表層土壤（土壤水分觀測深度：5 c m）達到飽和時，開始產流，產流時的降雨強度為2×10－6m／s（0.6 mm／5 min）。產流時，20 c m深處的土壤水分含量有了較明顯的增加，而50 c m處的土壤水分含量依然保持相對穩定狀態，此為長歷時低強度降雨條件下的產流機制。

2.2 短歷時高強度降雨條件下的產流機制

圖4所示為觀測斷面2008-08-29降雨、土壤水分變化及地表徑流產生的過程。此次降雨發生在當日午后14：00—15：00，平均降雨強度超過1.67×10－5m／s（5 mm／5 min）。由圖4可知，降雨開始10 min后開始產流，此時降雨強度為1.4×10－5m／s（4.2 mm／5 min）。產流時，除觀測斷面右側坡面5 c m處土壤水分略有增加外，其它各觀測點土壤水分保持穩定狀態。對隨機選取的多次降雨—徑流事件的分析表明：在短歷時條件下，降雨強度達到0.52 mm／min（2.6 mm／5 min）時，表層土壤在不飽和時便可以產流，結合對圖4所示降雨—徑流事件的分析結果，可以得出試驗流域在短歷時高強度降雨條件下，表層土壤未達到飽和狀態時的產流條件是降雨強度達到或超過0.52 mm／min（2.6 mm／5 min）。

圖3 長歷時低強度降雨條件下的產流事件

圖4 短歷時高強度降雨條件下的產流事件

2.3 降雨因子與徑流系數的關系

2.3.1 平均降雨強度與徑流系數的相關性 降雨強度、降雨歷時和降雨量等降雨因子都是徑流過程和徑流量的影響因子［14－15］，為了評價主要降雨因子對徑流的影響，隨機選取多次降雨—徑流事件，其降雨歷時、降雨量、平均雨強及徑流系數如表1所示。因為試驗流域所在地區只有在雨季集中降雨期間且滿足產流條件時才能產生地表徑流，所以本文所述的隨機降雨—徑流事件指在觀測得到的降雨—徑流事件中抽樣選取的事件。其中既包含短歷時的降雨事件又包含長歷時降雨事件。根據表1中數據繪制平均降雨強度和徑流系數的關系曲線如圖5所示，由該圖可知：平均降雨強度與徑流系數之間存在著較高的正相關性（R2＞0.87），即隨著平均降雨強度的增加，徑流系數呈增大的趨勢。

2.3.2 降雨歷時及降雨量與徑流系數的關系 根據表1中數據繪制的降雨歷時與徑流系數、降雨量與徑流系數的關系曲線分別如圖6、圖7所示，由圖6、圖7可知，降雨歷時與徑流系數之間沒有顯著的相關性，但存在輕微的反比例趨勢，即隨著單次降雨事件歷時的增加，徑流系數呈減小的趨勢，但不明顯。降雨量與徑流系數之間亦不存在顯著的相關性，二者之間也存在著輕微的反比趨勢，即隨著單次降雨事件降雨量的增加，徑流系數有輕微減小的趨勢。由以上分析可知，在平均降雨強度、降雨歷時和降雨量3個降雨特征因子中，平均降雨強度是影響徑流系數的主要因子，二者存在著較高的正相關性。

圖5 平均降雨強度－徑流系數關系曲線

圖6 降雨歷時－徑流系數關系曲線

圖7 降雨量－徑流系數關系曲線

3 降雨－徑流數值模型

3.1 模型構建

3.1.1 河網構建 利用 Arc GIS工具中的 Archydr ol gy和Arcview平臺，對流域地形進行分析，包括基于DEM流域特征的提取，水系的提取、子流域的劃分、流域邊界線的提取、地形指數的計算、單元匯流路徑長度的提取等［16－17］。依據空間上各溝道連接的相對位置關系將試驗流域分割成若干個分布式小流域，從而可以通過各個小流域的空間結合來表達試驗流域的整體。按最大坡度方向決定地表匯流方向［18］。擬河道網是在數字標高模型和水流方向線生成的基礎上，利用Arcview生成的，按上述方法生成的地表徑流觀測點上游集水區（計算區域）的河網如圖8所示。為了實現基于運動波理論建立的模型的徑流計算，將流域以擬河道網為依據進行分布式分割，將流域劃分成若干個分布式小流域。

3.1.2 水文地質模型 根據對試驗流域土壤縱剖面實際調查的結果，基于研究區實際水文地質條件，構建水文地質模型如圖9所示。模型由坡面區間和河道區間構成，坡面區間在縱斷面由兩層構成，河道區間在縱剖面上由一層構成。

圖9 水文地質模型

3.1.3 算法建立

（1）連續方程式。圖10所示為任意流域的單元體水收支示意圖，對于地表徑流，在計算時間步長Δt內以質量守恒定律（水收支平衡法）建立如下數學模型：

式中：b——單元體寬度（m）；d x——流向上的單位長度（m）；Δt——時 間 步 長 （s）；r——降 雨量 （m／s）；f1——第一層土壤平均滲透速度（m／s）；A——地表徑流的斷面面積（m2）；Q——流入量（m3／s）；（Q＋?Q／?x）——流出量（m3／s）。

因為A＝b·h，Q＝A·v［h為地表徑流深度（m）；v為流速（m／s）］，其他因子同前，對方程式進行整理兩邊同時約去b，得

按照同樣方法，可以推得第一層地下水的連續方程式為：

式中：λ——第一層的有效孔隙率；ˉq——單寬流量（m2／s）；ˉh——水深（m）；f1——第一層土壤平均滲透速度（m／s）；f2——由第一層向第二層的滲透速度（m／s）；E——蒸散發量（m／s）；蒸散發量認為從第一層產生（土壤水分）；t——時間（s）；x——水流方向上的距離（m）；用同樣的方法，如果土壤在縱剖面上按照土壤性質（包括含水層）可分成多個層，則可以建立任意層的水流運動的連續方程式：

式中：i——層號；λi——第i層的有效孔隙率；ˉqi——第i層單寬流量（m2／s）；ˉhi——第i層水深；fi——由第i－1層向第i層的滲透速度；fi＋1——由第i層向第i＋1層的滲透速度，降雨期間，因為植物的蒸散發可以忽略，同時，考慮到研究區的地表徑流只有在雨季集中降雨發生期間才能產生，而且存在時間很短等特點，所以，采用一個損失系數評價降雨期間的蒸散發［19］。

圖10 流域單元體水收支示意圖

（2）運動方程式。運動波理論是基于河道抵抗法則基礎之上，水流基礎方程式之一的運動方程式可以用等流的河道抵抗法則置換，即一維水流計算的運動方程式可以用平均流速公式來代替，因此，地表徑流計算的運動方程式為曼寧平均流速公式：

式中：h——水深（m）；q——單寬流量（m2／s）；n——曼寧粗度系數（s m－1／3）；R——水力半徑（m）；I——水力坡降（計算時用河道坡度近似代替），其它因子與前述一致。地下水（滲透）連續方程式為：

3.1.4 差分計算 在時間上和空間上用有限差分法對地表流和滲透流連續方程式進行離散化，實現連續計算。因為數值計算要依賴于流域最末級子流域端部的邊界條件和初始條件，所以采用后退差分法進行差分，差分公式如下：

地表流連續式的差分公式：

滲透流連續式的差分公式：

式中：n——計算的時間次序；i——水流方向上計算格子的編號；其它因子與前述相同［21］。

3.1.5 計算參數 計算所用主要參數如表2所示。參數由實地調查、水準測量、參考同一研究區的有關研究成果以及利用模型進行海量運算，對部分參數進行尋優等方法確定的。

表2 主要計算參數

3.2 數值模擬

（1）一次“降雨—徑流”事件的數值模擬：以2008年9月23—24日觀測的典型降雨和地表徑流為例，數值模擬的結果如圖11所示。誤差分析結果表明，模擬結果與觀測流量的誤差小于3%。

（2）數次“降雨—徑流”事件的數值模擬：以2004年8月11日至22日觀測的降雨及產流結果為例，數值模擬結果如圖12所示。誤差分析結果表明，觀測流量與模擬流量之間的誤差均在誤差判斷基準允許范圍之內（誤差＜3%）。

圖11 一次降雨徑流事件模擬（2008-09-23－24）

圖12 數次降雨徑流模擬結果（2004-08-11－22）

3.3 數值計算

利用降雨—徑流數值模型和觀測的降雨量，對試驗流域2005—2009年（5 a）的降雨徑流過程分別進行各年全年計算，結果如表3所示。

表3 2005－2009年降雨及年徑流系數

由表3可知：2005—2009年的平均年降雨量為439 mm／a，與多年平均降雨量437 mm／a相近，5 a的年平均徑流系數為0.11，與在同一區域用長期觀測的方法而得到的結果幾乎相等，可以近似地推得研究流域所在地區的年徑流系數應該在0.10～0.15?；谝陨辖Y論可知：本研究所開發的地表徑流計算模型效率較高，可以為該地流域提供一種較為精確的數值計算方法。

4 結論

本研究基于對觀測水文數據的分析，揭示了研究區在長歷時低強度降雨以及短歷時高強度降雨條件下的產流機制；開發了適用于黃土高原北部水蝕風蝕交錯區降雨—徑流過程的數值計算模型，利用試驗流域2005—2009年的降雨數據進行了計算并對結果進行了分析，得出的結論如下：

（1）試驗流域在長歷時低強度條件降雨條件下，產流的必要條件是表層土壤達到飽和且降雨強度達到或超過0.12 mm／min（0.6 mm／5 min）；在短歷時高強度降雨條件下，產流的條件是降雨強度達到或者超過0.52 mm／min（2.6 mm／5 min）。

（2）在平均降雨強度、降雨歷時和降雨量3個降雨因子中，平均降雨強度是影響徑流系數的主要因子，其與徑流系數之間存在較好的正相關性；降雨歷時和降雨量與降雨徑流系數之間只存在著輕微的反比趨勢。

（3）試驗流域2005—2009年的年平均徑流系數為0.11，根據2005—2009年的年平均降雨量與試驗流域所在地區多年平均降雨量幾乎相等的事實，近似地推得試驗流域所在地區的多年徑流系數為0.1～0.15。

（4）研究開發的分布式流域降雨—徑流數值模型適用于試驗流域所在地區的水文地質條件，因此可在同一地區中小尺度流域的降雨徑流計算中推廣應用。

［1］ 孫穎娜，芮孝芳，付強，等.基于隨機微分方程的流域產匯流模型［J］.水利學報，2011，42（2）：187-191.

［2］ Hinokidani O，Huang J B，Yasuda H，et al.Annual water budget of a small basin in the norther n Loess Plateau in China［J］.Journal of Arid Land Studies，2010，20（3）：167-172.

［3］ 高軍俠，劉作新，黨宏斌.黃土高原坡面模擬降雨超滲徑流特征分析［J］.土壤通報，2004，36（6）：780-784.

［4］ 周宏飛，王大慶，馬健，等.天山山區草地覆被和雨強對產流和產沙的影響研究［J］.水土保持通報，2009，29（5）：26-29.

［5］ 李慧敏，張建軍，徐佳佳，等.黃土高原不同下墊面小流域徑流特征［J］.東北林業大學學報，2011，39（10）：90-113.

［6］ 芮孝芳，黃國如.分布式水文模型的現狀與未來［J］.水利水電科技進展，2004，24（2）：55-58.

［7］ Hinokidani O，Huang J B，Yasuda H，et al.Study on surface r unoff characteristics of a s mall ephemeral catchment in the norther n Loess Plateau，China［J］.Jour nal of Arid Land Studies，2010，20（3）：173-178.

［8］ 王幼奇，樊軍，邵明安，等.黃土高原水蝕風蝕交錯區三種植被蒸散特征 ［J］.生態學報，2009，29（10）：5384-5396.

［9］ 唐克麗，侯慶春，王斌科，等.黃土高原水蝕風蝕交錯帶和神木試區的環境背景及整治方向［J］.中國科學院水利部西北水土保持研究所集刊，1993，18：1-5.

［10］ 王萬忠，焦菊英，郝小品.黃土高原暴雨空間分布的不均勻性及點面關系［J］.水科學進展，1999，10（2）：165-169.

［11］ 李免，李占斌，劉普靈，等.黃土高原水蝕風蝕交錯帶土壤侵蝕坡向分異特征［J］.水土保持學報，2004，18（1）：63-65.

［12］ Wang Q X，Takahashi H.A land surface water deficit model for an arid and semiarid region：i mpact of desertification on the water deficit status in the Loess Plateau，China［J］.Journal of Cli mate，1999，12（1）：244-257.

［13］ 王云強，張興昌，從偉，等.黃土區不同土地利用方式坡面土壤含水率的空間變異研究［J］.農業工程學報，2005，19（1）：65-71.

［14］ 王貴作，任立良.基于柵格垂向混合產流機制的分布式水文模型［J］.河海大學學報：自然科學版，2009，37（4）：386-390.

［15］ 焦菊英，王萬忠，郝小品.黃土高原不同類型暴雨的降水侵蝕特征［J］.干旱區資源與環境，1999，13（1）：34-42.

［16］ Huang J B，Hinokidani O，Yasuda H，et al.Study on characteristics of the surface flow of the upstream region in Loess Plateau［J］.Annual Jour nal of Hydraulic Engineering，JSCE，2008，52：1-6.

［17］ 任立良，劉新仁.基于數字流域的水文過程模擬研究［J］.自然災害學報，2001，9（4）：45-52.

［18］ 李忠武，蔡強國，曾光明，等.基于GIS的黃土丘陵溝壑區土壤水分模型研究［J］.水利學報，2004（3）：123-128.

［19］ 黃金柏，檜谷治，梶川勇樹，等.分步型流域“降雨—流出”過程數值模擬方法的研究［J］.水土保持學報，2008，22（4）：52-55.

［20］ Kobayashi K，Takara K，Tachikawa Y.Parameter estimation of a distributed rainfall-runoff model by a levenber g-marquardt opti mization algorith m［J］.Annual Jour nal of Hydraulic Engineering，JSCE，2007，51：409-414.

［21］ 汪志榮，沈晉，王文焰，等.黃土區波涌畦灌條件下地表水流運動實驗與數值模擬［J］.農業工程學報，1994，10（1）：36-43.