土壤滲透系數(shù)各向異性對我國西南山區(qū)小流域產(chǎn)匯流過程的影響研究

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

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

1 研究背景

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

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

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

2 研究區(qū)域

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

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

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

3 數(shù)據(jù)與研究方法

3.1 數(shù) 據(jù)

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

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

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

3.2 InHM模型介紹

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

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

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

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

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

關(guān)于InHM 的詳細(xì)描述，請參閱Vander Kwaak［23］和Heppner［25］。

3.3 模型設(shè)定

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

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

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

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

土壤含水率與壓強(qiáng)水頭及相對滲透率的關(guān)系基于van Genuchten 方法［37］由模型率定獲得。地表的Manning 糙率系數(shù)由下式［38］計(jì)算：

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

本研究中，堿坪溝流域的植被覆蓋率假定為90%，經(jīng)計(jì)算Manning系數(shù)為0.6。

為了使正式模擬開始時(shí)流域內(nèi)的含水量盡可能接近流域?qū)嶋H情況，模擬前先對全流域進(jìn)行20 d 的預(yù)模擬，選取流域出口處流量與實(shí)測流量近似的時(shí)刻作為正式模擬研究的初始時(shí)刻。

3.4 模型率定驗(yàn)證

考慮到像堿坪溝這樣的山區(qū)流域存在山洪暴發(fā)的風(fēng)險(xiǎn)，該模型按洪水事件進(jìn)行模擬，以便更好地了解對山洪的影響。為此，選取了7 場洪水場次，其中4 場用于率定，3 場作為驗(yàn)證工況。洪峰流量作為山洪預(yù)報(bào)的關(guān)鍵特征之一，以此為根據(jù)進(jìn)行了率定和驗(yàn)證。在實(shí)測徑流特征數(shù)據(jù)的基礎(chǔ)上，通過調(diào)整土壤特征曲線參數(shù)，模擬了流域出口處的洪峰流量。圖3展示了所選事件的觀測值和模擬值的對比，Qobs為觀測值，Qmod為模擬值?？梢钥闯?，該模型在水文過程模擬中表現(xiàn)較好，NSE值分別達(dá)到了0.88 和0.89，模擬結(jié)果基本能夠代表實(shí)際水文響應(yīng)結(jié)果，可以此為基礎(chǔ)進(jìn)行進(jìn)一步的分析研究。率定參數(shù)α和β分別為0.13和2.2。

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

3.5 工況設(shè)置

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

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

4 結(jié)果與討論

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

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

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

此外，隨著Kr的增加，初始水文響應(yīng)過程對相同降雨事件變得更加敏感。當(dāng)Kr小于1 時(shí)，初始徑流量波動(dòng)不大；當(dāng)Kr大于1 時(shí)，初始水文響應(yīng)出現(xiàn)一個(gè)較小的流量峰值。這種敏感性的上升可能是由于橫向土壤滲透系數(shù)的增加，橫向?qū)栽酱螅髟饺菀自谕寥乐兴搅鲃?dòng)，土壤含水層中儲(chǔ)存的水分也就會(huì)越多地出流，使得起始時(shí)刻流量迅速增加。因此，若不考慮各向異性，容易對初始水文響應(yīng)產(chǎn)生不準(zhǔn)確的估計(jì)。

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

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

4.2 滲透系數(shù)各向異性對產(chǎn)流過程的影響

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

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

4.2.2 滲透系數(shù)各向異性對坡面土壤水分運(yùn)動(dòng)的影響

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

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

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

圖7很好地解釋了圖6的現(xiàn)象，當(dāng)Kr增加時(shí)，在坡面上水流橫向流動(dòng)速度加快，流速方向也逐漸變?yōu)樗椒较颍黾恿似旅媾潘?，進(jìn)而促進(jìn)了地表水的入滲，使地表水深進(jìn)一步降低，大量地表水以及土壤中儲(chǔ)存的水分出滲到河道形成徑流，最終抬高了河道的水深。

4.2.3 滲透系數(shù)各向異性對河道土壤水分運(yùn)動(dòng)的影響

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

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

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

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

5 結(jié) 論

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

結(jié)果表明，土壤透水碎石層的滲透系數(shù)各向異性對研究區(qū)域的產(chǎn)匯流過程有較大的影響。洪峰流量和峰現(xiàn)時(shí)間均隨Kr而變化，當(dāng)Kr為10 時(shí)，峰值流量比各向同性工況下高出約70%，峰現(xiàn)時(shí)間比Kr=0.1 的工況減少近26 min。隨后的模擬進(jìn)一步表明，這些影響是由土壤水分運(yùn)動(dòng)的流速和水流方向變化導(dǎo)致的。隨著Kr的增大，流速顯著增大，尤其是坡面橫向水流運(yùn)動(dòng)，更多的土壤水分以地下暴雨徑流的形式水平流動(dòng)，從而促進(jìn)了降雨自地表下滲。當(dāng)?shù)叵卤┯陱搅鞯竭_(dá)河道時(shí)，較大的Kr也會(huì)加速其出滲，使得流量增加。土壤水分運(yùn)動(dòng)速度的提高也是縮短峰現(xiàn)時(shí)間的原因。

由此可見，在不考慮各向異性的情況下進(jìn)行數(shù)值模擬可能會(huì)嚴(yán)重低估洪澇災(zāi)害的危害程度，對峰現(xiàn)時(shí)間的預(yù)警也不夠準(zhǔn)確，從而影響疏散調(diào)度。因此，土壤滲透系數(shù)各向異性的特性不能被忽略，合理考慮各向異性比能夠提高山洪模擬的精確度，在山洪暴發(fā)危及生命財(cái)產(chǎn)安全的山區(qū)流域，研究土壤滲透系數(shù)各向異性對防洪預(yù)警具有重要意義?！?/p>