土壤滲透系數各向異性對我國西南山區小流域產匯流過程的影響研究

陳蒼乙，冉啟華，劉 琳，潘海龍，葉 盛

（浙江大學建筑工程學院水文與水資源工程研究所，杭州310058）

1 研究背景

山區流域地勢陡峭，暴雨頻發，容易在短時間內形成較大洪水，可能導致山洪、泥石流等地質災害，嚴重威脅著人民群眾的生命財產安全。山區流域降雨產流水文響應過程是山洪形成的基礎，也是目前山洪相關研究的核心科學問題之一。鑒于山洪成災快、破壞力強等特點，也為更全面地認識這些因素對山洪的影響，數值模擬常作為一種有效手段來研究山洪發生過程。土壤飽和滲透系數表征了土壤的下滲能力，是影響降雨徑流過程非常重要的物理量，通常被作為數值模擬的一個關鍵參數。在自然界中，土壤滲透系數通常具有各向異性，一般用各向異性比（Kr）來表示，即水平方向（Kh）和垂直方向（Kv）飽和滲透系數的比值（Kh/Kv），可能是由土層層狀或片狀構造［1］、沉積過程和復雜的應力加載［2］等原因所導致的。關于土壤滲透系數各向異性的具體大小，目前學者的研究尚未形成一致。有研究表明，土壤的滲透系數在垂直方向上高于水平方向［3-5］；但也有研究得出相反的結論，水平方向的滲透系數高于垂直方向［6，7］。Terzaghi和Scott發現，重復層狀土壤的滲透系數比較大，可達到100［8，9］；C Germa 'n Soracco 等［10］也發現土壤樣品頂層15 cm 存在各向異性，水平方向的滲透系數約為垂直方向的5倍。

針對土壤滲透系數各向異性帶來的影響，現有已有一些學者對此進行了研究。胡瑞等［11］發現在各向異性假設下，隨著水平滲透系數的增加，單寬流量明顯增大。周鵬鵬等［12］建立了我國南方某地下水水流數值模型，結果表明，隨著各向異性比Kr的增大，洞室涌水量增加。曾文西等［13］發現在同一雨型作用下，Kr越小，雨水入滲深度越大。Yu等［14］指出Kr增加會導致邊坡入滲能力和安全系數的降低。然而由于流域性質和水文過程的復雜性和異質性［15］，較多水文模擬案例中為簡化模型分析計算過程均將土壤滲透系數設定為各向同性（Kh/Kv=1）。而上述相關研究表明，這種簡化可能會導致無法準確估計流量、洞內涌水量和邊坡安全系數等結果，進而可能造成嚴重的社會、經濟后果。

為進一步明確土壤滲透系數各向異性在降雨徑流過程中的影響，本文運用一種基于物理概念的水文數值模型，對我國西南山區堿坪溝流域的降雨徑流過程進行了模擬。設置了不同土壤滲透系數各向異性的情景，系統性分析了滲透系數各向異性對流域土壤水分運動以及產匯流過程的影響。

2 研究區域

考慮到山洪災害的頻率、汶川地震影響等因素，本研究選用位于四川省都江堰市西北部的龍溪河流域作為研究區域，并選擇龍溪河流域內的支溝堿坪溝流域進行詳細研究。

龍溪河屬于長江流域岷江水系（圖1），堿坪溝為龍溪河中游左岸的次級支溝［16］。堿坪溝流域面積約3.5 km2，海拔為1 050～2 199 m，坡面的平均坡度為31.7%。整體自東北向西南方向延伸，主溝長1.7 km，整個流域被森林覆蓋，平均縱比降36.3%［17］。堿坪溝地質構造復雜，斷裂較為發育，出露的地層主要有震旦系灰綠色安山巖、凝灰巖及安山玄武巖和第四系地層，其中第四系地層主要為第四系殘坡積層，巖性為含碎石粉質黏土［18］。氣候為亞熱帶濕潤氣候，都江堰市年降水量約1 134.8 mm，主要集中在5-9月，占到全年降水量的80%以上，產流方式以蓄滿產流和地下潛流為主［16，19］。

圖1 研究區域地理位置及水系分布Fig.1 Location and river network of study area

3 數據與研究方法

3.1 數 據

本研究使用的水文數據包括堿坪溝流域2018年6月24日至7月21日雨季實測降雨和徑流數據。降雨數據由JBD-2 翻斗式雨量計每1 min 采集一次。流量數據利用基于大尺度粒子圖像測速（簡稱LSPIV）和立體視覺（簡稱SI）的流速面積法計算得到，測算精度為45 min/次，在本研究中將該方法所得的流量值作為實測流量數據使用。該方法將大尺度粒子圖像測速技術與立體視覺技術相結合，組成一套完整的測流計算系統（SILSPIV），旨在實時非接觸式獲取河道斷面的水位、表面流場和斷面流量等要素［20］。SI-LSPIV 系統用于流量測算的準確性已在之前的研究中得到了驗證［21，22］。

為了解堿坪溝流域內土壤滲透系數實際情況，在綜合考慮流域地形及道路等客觀因素后，研究人員從上游到流域出口共選了5 處較有代表性的地點采集土壤樣品，沿深度方向在每個點的10～50 cm 和50～100 cm 兩個范圍內使用環刀分別從豎直向和水平向采集用于測量土壤飽和滲透系數和孔隙度的原狀土樣，密封保存后帶回實驗室。采用恒定水頭法測定土壤飽和滲透系數，用飽和浸泡法測定土壤孔隙率。受地質條件限制，土壤采樣難度較大，僅在其中4 處地點同時采集到了用于縱向滲透系數（Kv）和橫向滲透系數（Kh）測試的原狀土樣。實測結果（表1）表明每組樣品的縱向和橫向滲透系數均有差異，由此可以得出，堿坪溝流域土壤滲透系數確實存在各向異性。

表1 各采樣點土壤縱向和橫向滲透系數對比Tab.1 Comparison between vertical and horizontal hydraulic conductivity of sampling points

3.2 InHM模型介紹

本文所用的模型為基于物理概念的分布式水文模型InHM（Integrated Hydrology Model），該模型最初在1999年由Waterloo大學的Vander-Kwaak 開發［23］。模型能夠在不預設產流機制的前提下，用基本物理定律來描述產流的自然過程。可模擬地下水在飽和和非飽和區的三維流動、地下水在大孔隙中的三維流動、地表水二維坡面流動及二維河道流動等［24］。目前已被成功應用于泥沙運動模擬［25-28］、流域尺度下的地貌演變研究［29］、小尺度下水文響應模擬［30-34］等領域。

地表的瞬時流（包括坡面流和明渠流）是通過二維的淺水方程的擴散波近似計算。該二維地表水流被概化為第二種連續介質，通過一層厚度為as的水體覆蓋在地下土體上與地下飽和/非飽和多孔介質互動。地表淺水深度為ψs，則地表水流的流動可以表述為：

式中：qb為各種邊界上的源/匯速率，s-1；qs為地表水流速度m/s；qe為地表和地下的水量交換速率，s-1；t為時間，s；hs為網格中未顯示的地表微地形的平均高度，m；Sws為地表飽和度；ψsmobile和ψsstore為貯存水深和流動水深，m。

地下水在地下孔隙介質飽和/非飽和土體和大孔隙（兩者均為地下連續體）中的流動在InHM 模型計算中通過Richards 方程表示：

式中：φ為孔隙率；Sw為土壤飽和度；fa為各介質所對應的面積百分比，%；fv為各介質所對應的體積百分比，%；q?為達西流量，m/s。

關于InHM 的詳細描述，請參閱Vander Kwaak［23］和Heppner［25］。

3.3 模型設定

基于堿坪溝流域的5 m 精度DEM 數據，提取出流域的河道，并在此基礎上生成模型所需的三維網格（圖2）。流域邊界的網格水平精度為100 m，河道的網格水平精度為10 m［35］。根據對流域的實地觀測，網格垂向上總共分為三層：表層為第四紀地層土壤（0～0.5 m），網格垂向精度為0.1 m；中間為透水碎石層（0.5～3 m），網格垂向精度為0.5 m；底部為不透水基巖層（3～13 m），網格垂向精度為2 m。

圖2 堿坪溝流域模型網格Fig.2 Finite-element meshes for the Jianpinggou catchment

對龍溪流域的土壤調查表明，土壤性質與植被類型密切相關［36］。堿坪溝流域中的植被分布與海拔高度相關，因此將表層土壤按照海拔高度進行分區：1 300 m 以上、1 180～1 300 m 和1 180 m 以下。模型中表層和碎石層土壤的飽和滲透系數和孔隙率從土壤樣品實測值中取平均值，基巖層的飽和滲透系數和孔隙率分別設定為1×10-8m/s和0.3（表2）。

表2 模型設定中的土壤飽和滲透系數和孔隙率Tab.2 The porosity and saturated hydraulic conductivity used in the model

土壤含水率與壓強水頭及相對滲透率的關系基于van Genuchten 方法［37］由模型率定獲得。地表的Manning 糙率系數由下式［38］計算：

式中：Cv為地表植被覆蓋率，以小數形式表示。

本研究中，堿坪溝流域的植被覆蓋率假定為90%，經計算Manning系數為0.6。

為了使正式模擬開始時流域內的含水量盡可能接近流域實際情況，模擬前先對全流域進行20 d 的預模擬，選取流域出口處流量與實測流量近似的時刻作為正式模擬研究的初始時刻。

3.4 模型率定驗證

考慮到像堿坪溝這樣的山區流域存在山洪暴發的風險，該模型按洪水事件進行模擬，以便更好地了解對山洪的影響。為此，選取了7 場洪水場次，其中4 場用于率定，3 場作為驗證工況。洪峰流量作為山洪預報的關鍵特征之一，以此為根據進行了率定和驗證。在實測徑流特征數據的基礎上，通過調整土壤特征曲線參數，模擬了流域出口處的洪峰流量。圖3展示了所選事件的觀測值和模擬值的對比，Qobs為觀測值，Qmod為模擬值。可以看出，該模型在水文過程模擬中表現較好，NSE值分別達到了0.88 和0.89，模擬結果基本能夠代表實際水文響應結果，可以此為基礎進行進一步的分析研究。率定參數α和β分別為0.13和2.2。

圖3 洪峰流量觀測值與模擬值的對比Fig.3 Comparison between the observed and simulated peak runoff for the selected events

3.5 工況設置

為了研究土壤滲透系數各向異性對產匯流的影響，我們對各向異性比進行了敏感性分析。由表2可知，中間層的滲透系數相對較大，對土壤水分運動的影響也相對更突出，因此本文重點研究中間透水碎石層土壤滲透系數各向異性的影響。根據表1實測數據各向異性比大小，本研究對碎石層設定了6 種不同各向異性比的工況以覆蓋實際可能出現的情況，以及1 種滲透系數各向同性（Kh/Kv=1）的基礎工況（表3）。每種工況均采用相同的降雨事件。

表3 不同各向異性比工況Tab.3 The simulated scenarios with different anisotropy

4 結果與討論

4.1 滲透系數各向異性對匯流過程的影響

圖4為不同Kr值下流域出口處的匯流過程流量曲線。從圖中可以看出，不同Kr值對流域出流有很大的影響，且隨著Kr的增大，峰值流量和總流量均增大。在滲透系數各向同性情況下（基礎工況），峰值流量為5 m3/s，總流量為5 138.6 m3。當Kr為10 時，峰值流量可達8.5 m3/s，比基礎工況高出約70%；總流量可達14 350.3 m3，是基礎工況的近3 倍。當Kr很小（即Kr=1/5）時，總流量為2 413.4 m3，不及各向同性工況下總流量的1/2。由此可見，在不改變其他土壤參數的情況下，匯流過程對滲透系數各向異性非常敏感。這一匯流特性可認為與該流域下墊面條件有關。改變土壤橫向滲透系數可能會使得地表徑流和地下徑流交互水量改變，也可以改變土壤水分運動的流速，峰值流量和流域出口處總徑流量因此而改變。

圖4 不同Kr下流域出口處匯流過程線Fig.4 Simulated hydrograph at the outlet under different scenarios

此外，隨著Kr的增加，初始水文響應過程對相同降雨事件變得更加敏感。當Kr小于1 時，初始徑流量波動不大；當Kr大于1 時，初始水文響應出現一個較小的流量峰值。這種敏感性的上升可能是由于橫向土壤滲透系數的增加，橫向導水性越大，水流越容易在土壤中水平流動，土壤含水層中儲存的水分也就會越多地出流，使得起始時刻流量迅速增加。因此，若不考慮各向異性，容易對初始水文響應產生不準確的估計。

圖5為Kr與洪水事件特征參數之間的散點圖。隨著Kr的增加，到達峰值流量的時間逐漸減少。對比工況1（Kr=10）和工況6（Kr=1/10），峰現時間相差1 586 s，約為26 min。這段時間對于山洪的早期預警是至關重要的，在山區洪水監測預報的實際應用中，如果不考慮土壤滲透系數的各向異性，可能導致預警延遲，對人身財產安全造成嚴重的危害。

圖5 不同Kr下的峰現時間、峰值流量及總流量Fig.5 Scatter plots between the anisotropic ratio（Kr）and time to peak discharge，peak discharge and total discharge

4.2 滲透系數各向異性對產流過程的影響

4.2.1 滲透系數各向異性對地表水深空間分布的影響

由圖6所示的地表水深空間分布可知，地表水深D也隨Kr變化。隨著Kr的增大，流域源頭地表水深逐漸降低，因為水流橫向流動速度增大時，坡面上的水下滲進入土壤中后加快匯聚至河道。同時，從圖6也可以看出，雖然地表產流面積較小，但河道水深較大。在滲透系數各向同性工況下，流域出口處水深為0.89 m；當Kr為10 時，水深1.43 m；當Kr為1/10 時，水深0.77 m。由此可見，在不同Kr值情況下，流域出口處水深差值相差較大。考慮到堿坪溝陡峭的地形條件，河道周邊區域通常是海拔相對較低處，坡面上的降雨入滲后，形成的地下暴雨徑流自坡頂匯流至河道及其周邊并出滲，匯集到河道中，使得河道及其周邊地表水深較大。

4.2.2 滲透系數各向異性對坡面土壤水分運動的影響

為了直觀展示土壤滲透系數各向異性帶來的影響，圖7給出了峰現時刻剖面AA'處［圖1（c）］的流速和流向分布情況。剖面AA'位于上游兩條支溝匯合處河道南岸坡面上，與下游主河道垂直。由于斷面處于相對平坦的區域，位置水頭的影響有限，水流運動主要由壓力水頭驅動。圖中標有箭頭的曲線表示水流方向，流速用不同的顏色體現；橫向實線為土壤表層、中間透水碎石層和不透水基巖層的分界線。

由圖7可看出，各工況下流速均在中間透水碎石層與底部不透水層交界處達到最大，由地表至交界處流速逐漸增大，然后往更深層逐漸減小。滲透系數各向異性比的不同對流速有著顯著影響，水流速度隨著Kr的增大而增大，最大和最小Kr的工況各自對應的最大流速相差近一個數量級。

除了流速之外，水流運動方向也隨Kr變化。在基礎工況中［圖7（d）］，土壤表層和中間碎石層中的水下滲至不透水基巖層的交界面，此時流速達到最大，隨后沿著地形坡度流動。但隨著Kr的增大，中間碎石層與不透水基巖層交界面處的流速明顯增大，流向逐漸向偏向于平行地表的方向，因而土壤水分可能會被限制在一個相對較淺的深度［39］［圖7（a），（b）］。而在Kr較小的情況下，水流基本不參與水平流動，大部分水體均垂直于地表入滲至更深的地方［圖7（e）-（g）］。

圖7很好地解釋了圖6的現象，當Kr增加時，在坡面上水流橫向流動速度加快，流速方向也逐漸變為水平方向，增加了坡面排水量，進而促進了地表水的入滲，使地表水深進一步降低，大量地表水以及土壤中儲存的水分出滲到河道形成徑流，最終抬高了河道的水深。

4.2.3 滲透系數各向異性對河道土壤水分運動的影響

為了研究河道河床處的土壤水分運動，本文進一步研究了主河道縱剖面的流速分布［圖1（c）中的BB'］，圖8為主河道BB'縱剖面流速示意圖。從圖中可以看出，河道地形呈階梯狀，與圖7中相對平緩的山坡不同，沿主河道的流速變化幅度較明顯，在陡坡處流速較大，在平緩臺面處流速較小。這種變化主要發生在表層和中間透水碎石層，基巖層的速度變化不大。隨著Kr的增大，透水碎石層的流速逐漸增大。各向同性工況下，河道坡面地下暴雨徑流的最大流速為2×10-4m/s；當Kr為10時，最大速度可以達到2×10-3m/s；而當Kr為0.1 時，最大流速可僅為9×10-5m/s。這種現象不論是在坡面上還是平臺處均出現，說明了河道沿程的流速整體均隨Kr增加而增加。

圖7 AA′縱剖面流速流向示意圖Fig.7 Longitudinal profile of flow velocity and flow direction under different Kr scenarios

圖8（h）是基礎工況下［圖8（d）］的局部放大圖，選取了一段包含了連續坡面及平臺的地形，圖上箭頭表示為水流流速方向。結果表明，在河道中出滲和入滲呈交替分布。出滲區域位于陡峭坡面處，由于地勢相對較高，其總水頭較高，土壤中的水容易排出，從而流速較大；而在地勢平緩處，水流匯聚，使得這些區域水流由于水頭影響更容易產生入滲。這種土壤水分運動規律也符合山區河流中常見的典型階梯式深潭系統。因此，隨著Kr的增大，坡面土壤中的水分會出滲到河道中，使得河道中的流量增大。這也進一步證實了前文對于圖6的解釋，即Kr較大時主河道水深較高。

圖8 主河道BB′縱剖面流速示意圖Fig.8 Longitudinal profile of flow velocity along the stream channel BB′under different Kr scenarios

5 結 論

土壤滲透系數各向異性在自然界中普遍存在，但其在數值模擬中往往被簡化為各向同性情形。本文將基于物理概念的分布式水文模型（InHM）應用于中國西南部山區流域，研究了滲透系數各向異性對降雨產匯流過程造成的影響。

結果表明，土壤透水碎石層的滲透系數各向異性對研究區域的產匯流過程有較大的影響。洪峰流量和峰現時間均隨Kr而變化，當Kr為10 時，峰值流量比各向同性工況下高出約70%，峰現時間比Kr=0.1 的工況減少近26 min。隨后的模擬進一步表明，這些影響是由土壤水分運動的流速和水流方向變化導致的。隨著Kr的增大，流速顯著增大，尤其是坡面橫向水流運動，更多的土壤水分以地下暴雨徑流的形式水平流動，從而促進了降雨自地表下滲。當地下暴雨徑流到達河道時，較大的Kr也會加速其出滲，使得流量增加。土壤水分運動速度的提高也是縮短峰現時間的原因。

由此可見，在不考慮各向異性的情況下進行數值模擬可能會嚴重低估洪澇災害的危害程度，對峰現時間的預警也不夠準確，從而影響疏散調度。因此，土壤滲透系數各向異性的特性不能被忽略，合理考慮各向異性比能夠提高山洪模擬的精確度，在山洪暴發危及生命財產安全的山區流域，研究土壤滲透系數各向異性對防洪預警具有重要意義。□