摘要： 為了揭示不同大孔隙參數(shù)取值條件下邊坡水分場的變化規(guī)律，利用2個耦合的Richards方程描述大孔隙流，并聯(lián)合運(yùn)動波方程，建立坡面徑流與大孔隙流耦合模型；基于有限元軟件COMSOL的偏微分方程接口，實(shí)現(xiàn)所建立耦合模型的數(shù)值求解，并自行設(shè)計室內(nèi)模型試驗(yàn)，驗(yàn)證數(shù)值結(jié)果，分析不同大孔隙參數(shù)對體積含水率和坡面積水深度的影響。結(jié)果表明：COMSOL軟件的偏微分方程接口可以實(shí)現(xiàn)所建立耦合模型的數(shù)值求解，相比于無大孔隙坡面，考慮大孔隙時水分入滲深度明顯更大；邊坡飽和區(qū)深度與濕潤鋒深度均隨大孔隙域占比的增大而增大，均隨大孔隙域與基質(zhì)域?qū)禂?shù)之比與經(jīng)驗(yàn)參數(shù)的增大而減小；當(dāng)降雨歷時為30 min時，坡面積水深度區(qū)分度最大，表現(xiàn)在隨著大孔隙域占比的增大而減小，均隨大孔隙域與基質(zhì)域?qū)禂?shù)之比和經(jīng)驗(yàn)參數(shù)的增大而增大； 3種大孔隙參數(shù)按照對邊坡水分場影響程度由大到小的順序?yàn)榇罂紫队蛘急?、?jīng)驗(yàn)參數(shù)、大孔隙域與基質(zhì)域?qū)禂?shù)之比。

關(guān)鍵詞： 坡面徑流； 水分場； 有限元法； 大孔隙流； 濕潤鋒

中圖分類號： U416.1

文獻(xiàn)標(biāo)志碼： A

Abstract： To reveal change law of slope water fields under the condition of different macropore parameters， two coupled Richards equations were used to describe macropore flows， and a coupled model of slope runoffs and macropore flows was established by combining kinematic wave equations. Partial differential equation interface of finite element software COMSOLwas used to realize numerical solving of the established coupled model， and an innovative laboratory test was conducted to verify numerical results. Influences of different macropore parameters on water volume fraction and ponding depth of slope area were analyzed. The results show that the partial differential equation interface of COMSOL software can realize the numerical solving of the established coupled model. Compared with slopes without macropores， water infiltrationdepth is significantly greater when the macropores are considered. Depths of saturation zones and wetting fronts increase with the increase of macropore domain ratio， and decrease with the increase of transmissivity ratio of macropore

收稿日期： 2022-03-25 網(wǎng)絡(luò)首發(fā)時間：2023-03-10T11∶37∶32

基金項(xiàng)目： 國家自然科學(xué)基金項(xiàng)目（41772297）； 福建省中青年教師教育科研項(xiàng)目（JAT210542）

第一作者簡介： 李尚輝（1989—），男，湖北黃石人。講師，博士研究生，研究方向?yàn)檫吰聺B流與穩(wěn)定。Email： 451134454@qq.com。

通信作者簡介： 闕云（1980—），男，江西撫州人。教授，博士，博士生導(dǎo)師，研究方向?yàn)閹r土工程防災(zāi)減災(zāi)。Email： Queyun_2001@fzu.edu.cn。

中國東南沿海省份福建、 浙江一帶是世界上臺風(fēng)暴雨最活躍的地區(qū)之一，每年受臺風(fēng)暴雨影響而誘發(fā)的邊坡失穩(wěn)現(xiàn)象屢見不鮮。如2016年臺風(fēng)“鲇魚”登錄福建省泉州市， 最大降雨強(qiáng)度高達(dá)62 mm/h， 多地發(fā)生了以滑坡為主的地質(zhì)災(zāi)害， 造成32人死亡［1］。 降雨誘發(fā)邊坡失穩(wěn)在時序上屬于降雨—坡體水分運(yùn)移—水力交互作用—邊坡失穩(wěn)的過程，因此準(zhǔn)確刻畫坡內(nèi)水分運(yùn)移的過程是認(rèn)識邊坡失穩(wěn)的關(guān)鍵。 許多學(xué)者基于傳統(tǒng)的非飽和流理論對水分運(yùn)移過程進(jìn)行了大量研究， 并取得了許多有益成果［2］。

土壤結(jié)構(gòu)的變異性最為復(fù)雜。福建省廣泛分布的花崗巖殘積土是一種典型的大孔隙材料［3］。另外，受植物根孔、動物活動和自然氣候變化條件下土壤干濕交替和凍融循環(huán)的影響，實(shí)際邊坡表面土壤內(nèi)也分布著許多大孔隙，進(jìn)而影響土壤孔隙結(jié)構(gòu)特性［4］。大孔隙的存在使得降雨作用下土壤大孔隙中的水流可繞過周圍大部分土壤基質(zhì)而在短時間內(nèi)到達(dá)土壤深層，即產(chǎn)生大孔隙流［5］。大孔隙流是一種優(yōu)先流，相比于傳統(tǒng)的非飽和流，大孔隙流的水分運(yùn)移速率更大，并且內(nèi)部存在水分交換現(xiàn)象，不適宜采用單一達(dá)西定律進(jìn)行描述［6］。目前對大孔隙流的研究主要集中在土壤學(xué)科［7］，關(guān)于巖土地質(zhì)災(zāi)害中的邊坡失穩(wěn)問題，很少研究考慮到邊坡中的大孔隙結(jié)構(gòu)。

臺風(fēng)暴雨強(qiáng)度遠(yuǎn)大于土壤飽和導(dǎo)水系數(shù)，此時降雨量部分入滲，坡面短時間內(nèi)即出現(xiàn)徑流［8］，根據(jù)飽和-非飽和滲流邊界條件，降雨入滲速率趨于土壤飽和導(dǎo)水系數(shù)。實(shí)際上，當(dāng)降雨強(qiáng)度大于飽和導(dǎo)水系數(shù)時，入滲與徑流同時發(fā)生，兩者互相影響，坡內(nèi)滲流與坡面徑流的物理問題是一個耦合過程。王樂等［9］利用Richards方程描述邊坡飽和-非飽和滲流，利用運(yùn)動波方程描述坡面徑流，建立了滲流-徑流聯(lián)合求解的三維有限元模型，通過與Geo-Seep法所得結(jié)果進(jìn)行對比，驗(yàn)證了聯(lián)合求解的有限元法所得累積入滲量符合Geo-Seep模擬值。Tian等［10］提出了耦合地下水滲流Richards方程與地表徑流運(yùn)動波方程的二維聯(lián)合數(shù)值模擬方法，結(jié)果表明，聯(lián)合數(shù)值模擬方法計算精度較高，結(jié)果可靠。韓同春等［11］基于坡面徑流控制方程和Moore雙層入滲理論，建立耦合模型，并利用Python語言編制相應(yīng)的計算程序，分析耦合條件下雙層結(jié)構(gòu)邊坡降雨產(chǎn)流過程。已有耦合模型僅限于非飽和流，并且通過流量交換的迭代計算實(shí)現(xiàn)二者間耦合的數(shù)值模擬，穩(wěn)定性差且計算效率低［12］。有限元軟件COMSOL在求解非線性偏微分方程（PDE）時收斂性好且計算效率高，內(nèi)置的Richards物理場接口在非飽和地下水入滲方程中得到廣泛應(yīng)用［13］，對于沒有專門內(nèi)置接口的物理問題，也可通過COMSOL軟件中的PDE接口進(jìn)行建模，從而實(shí)現(xiàn)耦合方程的數(shù)值解，為求解強(qiáng)非線性PDE組提供了便捷途徑［14-15］。

本文中針對中國沿海省份福建省臺風(fēng)強(qiáng)降雨特性及土壤中的大孔隙結(jié)構(gòu)，建立坡面徑流與大孔隙流耦合模型，基于COMSOL軟件的PDE接口，實(shí)現(xiàn)所建立耦合模型的數(shù)值求解，分析不同影響因素對水分場分布特性的影響。

1 坡面徑流與大孔隙流耦合模型

1.1 理論方程

1.1.1 坡面徑流控制方程

坡面徑流控制方程采用一維圣維南方程的簡化形式，即運(yùn)動波方程［16］

1.1.2 大孔隙流控制方程

在目前描述大孔隙流的理論模型中，雙重滲透模型最具有代表性。在雙重滲透模型中，假定土壤由基質(zhì)域和大孔隙域構(gòu)成，并假設(shè)水分在這2個域中均能流動［17］。圖1所示為傳統(tǒng)非飽和流與大孔隙流入滲示意圖。在傳統(tǒng)的非飽和流模式下，土壤只含基質(zhì)域，水流從土壤上部均勻入滲，濕潤鋒推進(jìn)速度相同。在大孔隙流模式下，當(dāng)水分向下流動時，基質(zhì)域與大孔隙域之間入滲速率明顯不同，并且2個域中水分相互交換，濕潤鋒不在同一界面，該現(xiàn)象體現(xiàn)了大孔隙流的主要特點(diǎn)。

1.2 坡面徑流與大孔隙流耦合過程

在降雨初始時刻，坡面處于干燥狀態(tài)，土表入滲率理論上趨于無窮大，降雨強(qiáng)度遠(yuǎn)小于入滲率，來水全部入滲，邊界條件為第二類邊界條件，即流量邊界。隨著降雨的進(jìn)行，尤其是在強(qiáng)降雨作用下，土表體積含水率迅速增大，入滲率快速減小，當(dāng)入滲率減至降雨強(qiáng)度時開始產(chǎn)流，坡面邊界轉(zhuǎn)化為第一類邊界，即水頭邊界，此時方程（1）、（2）構(gòu)成坡面徑流與大孔隙流耦合模型。

耦合過程如下：在同一時間步內(nèi)，坡面徑流與入滲是通過共同的壓力水頭和坡面入滲率進(jìn)行耦合的，并且耦合過程通過邊界發(fā)生，方程（2）的邊界條件由方程（1）確定；而方程（1）中的入滲率由方程（2）確定，具體方法是根據(jù)方程（2）中的孔隙水壓力和非飽和導(dǎo)水系數(shù)，利用達(dá)西定律進(jìn)行求解［19］。大孔隙流入滲率由基質(zhì)域入滲率與大孔隙域入滲率共同決定，2個域的入滲率均利用達(dá)西定律進(jìn)行計算，再根據(jù)各自占比，利用加權(quán)求和的方法進(jìn)行求解。大孔隙流入滲率求解公式為

坡面徑流與大孔隙流耦合流程如圖2所示。

1.3 坡面徑流與大孔隙流耦合模型求解

基于有限元理論的多物理場數(shù)值仿真軟件COMSOL數(shù)學(xué)分支下的PDE接口具有強(qiáng)大的非線性PDE求解能力，PDE接口提供了3種基本形式的PDE，即系數(shù)型、廣義型、弱解型。可以根據(jù)需要自定義PDE，并設(shè)置參數(shù)、變量和邊界條件。

選取系數(shù)型PDE實(shí)現(xiàn)所建立坡面徑流與大孔隙流耦合模型的數(shù)值求解， 系數(shù)型PDE在整個計算區(qū)域Ω中的主控制方程［15］帶有方向， 鑒于本文中所求因變量均為標(biāo)量， 在不考慮方向的基礎(chǔ)上， COMSOL軟件中系數(shù)型PDE在Ω中的主控制方程形式為

與其他有限元軟件相比，COMSOL軟件可以靈活地自定義變量和參數(shù)，不同物理場之間耦合性強(qiáng)，對復(fù)雜方程的求解具有一定的優(yōu)越性。運(yùn)用PDE接口必須預(yù)先定義方程中的各種系數(shù)和表達(dá)式，這些系數(shù)和表達(dá)式既可以是自定義參數(shù)，也可以是關(guān)于時空坐標(biāo)和因變量的函數(shù)，編寫域方程和邊界條件應(yīng)匹配COMSOL軟件的PDE界面的特殊形式。方程（1）、（2）轉(zhuǎn)換為系數(shù)型PDE后，有限元軟件COMSOL的PDE接口中系數(shù)和表達(dá)式的設(shè)置如表1所示。

2 大孔隙邊坡室內(nèi)試驗(yàn)

2.1 試驗(yàn)邊坡幾何尺寸

試驗(yàn)用土樣源自廈蓉高速公路漳州段擴(kuò)建工程里程樁樁號為K103+160的路段，經(jīng)檢測，土樣的天然體積含水率為10.2%。設(shè)計室內(nèi)大孔隙邊坡，對比試驗(yàn)監(jiān)測數(shù)據(jù)與數(shù)值計算結(jié)果。試驗(yàn)邊坡總長度為120 cm，其中坡頂和坡趾長度為20 cm，邊坡高度為60 cm，坡趾高度為20 cm，坡比值為1∶2，布設(shè)2個體積含水率監(jiān)測點(diǎn)A1、 A2，如圖3所示。大孔隙布置方式是填土之前在土壤內(nèi)部插入一根直徑為2 cm的不銹鋼空心管，填土完成后將粒徑為2 mm的石英砂灌入管中，采用螺旋上升的方式緩慢拔出鋼管，近似模擬邊坡大孔隙［20］。

2.2 試驗(yàn)裝置

圖4所示為降雨作用下大孔隙邊坡試驗(yàn)裝置布置。 試驗(yàn)裝置主要由有機(jī)玻璃模型箱主體、 底座、 降雨系統(tǒng)、 數(shù)據(jù)監(jiān)測裝置及拍攝裝置五大部分組成。降雨系統(tǒng)采用南京南林電子科技有限公司研發(fā)的NLJY-10型人工模擬降雨系統(tǒng)， 主要由噴頭、 水管支架、 壓力控制系統(tǒng)、 供水箱、 雨量筒及降雨控制系統(tǒng)組成。 數(shù)據(jù)監(jiān)測裝置包括土壤水分監(jiān)測儀、 徑流速率和累積徑流量收集裝置。 采用PC-2SQ型控制儀連接TDR-3型水分傳感器（量程為0～100%，精度為±2%）制作土壤水分監(jiān)測儀。 采用經(jīng)濟(jì)且常用的染色劑法測定徑流速率。 因?yàn)樽霞t色的溶液區(qū)分度高， 所以染色溶劑選用高錳酸鉀溶液。 測定原理是當(dāng)坡面水流從開始形成至穩(wěn)定狀態(tài)時，根據(jù)紫紅色溶液流過一定坡面距離所需的時間計算徑流速率［21］。利用自制的三角玻璃引流槽將坡面累積徑流收集于容器中，利用玻璃膠將三角玻璃引流槽粘接在有機(jī)玻璃模型箱出水口處。

2.3 坡面徑流與大孔隙流耦合模型驗(yàn)證

基于2.2節(jié)中的試驗(yàn)條件， 開展室內(nèi)邊坡降雨試驗(yàn)，其中降雨強(qiáng)度為30 mm/h， 降雨歷時為500 min。 圖5所示為監(jiān)測點(diǎn)A1、 A2體積含水率試驗(yàn)值與計算值對比。 由圖可知： 2個監(jiān)測點(diǎn)體積含水率的試驗(yàn)值與計算值變化趨勢大致吻合， 均表現(xiàn)為先保持不變，隨后迅速增大，臨近最大值時增長放緩，最終趨于穩(wěn)定。監(jiān)測點(diǎn)A1體積含水率響應(yīng)時間計算值與試驗(yàn)值分別為76、 81 min，計算值較試驗(yàn)值提前了6.2%；監(jiān)測點(diǎn)A2體積含水率響應(yīng)時間計算值與試驗(yàn)值分別為236、 221 min，計算值比試驗(yàn)值落后6.8%，2個監(jiān)測點(diǎn)體積含水率響應(yīng)時間相差較小。超過響應(yīng)時間后，體積含水率均迅速增大，臨近飽和狀態(tài)，同一時刻體積含水率計算值大于試驗(yàn)值，說明數(shù)值計算條件下更容易達(dá)到飽和狀態(tài)，主要原因是試驗(yàn)條件下土壤內(nèi)部存在氣體，對水分入滲存在阻滯作用?？傮w上，體積含水率試驗(yàn)結(jié)果與計算結(jié)果相互驗(yàn)證良好。

圖6所示為坡面累積徑流量、徑流速率試驗(yàn)值與計算值對比。由圖可知：累積徑流量試驗(yàn)值與計算值均隨著降雨歷時的增加而增大， 在降雨歷時為0～30 min時，累積徑流量試驗(yàn)值與計算值曲線幾乎重合，此后出現(xiàn)拐點(diǎn)，增長速率均顯著增大，并且計算值略大于試驗(yàn)值；當(dāng)降雨歷時為220 min時，累積徑流量試驗(yàn)值與計算值幾乎相等，分別為16 294、 16 456 cm3。徑流速率試驗(yàn)值與計算值曲線變化趨勢也大致相同，降雨歷時超過30 min后計算值趨于平穩(wěn)，約為0.011 2 m/s，而試驗(yàn)值則出現(xiàn)上下波動現(xiàn)象，流速波動界限約為0.010 9 m/s，因此累積徑流量、徑流速率的試驗(yàn)值與計算值可以相互驗(yàn)證。

由監(jiān)測點(diǎn)體積含水率、累積徑流量和徑流速率的對比可知，采用粒徑為2 mm的粗石英砂人工制造大孔隙的方法在一定程度上可以近似模擬大孔隙，試驗(yàn)結(jié)果與數(shù)值計算結(jié)果可以相互驗(yàn)證，表明利用COMSOL軟件的PDE接口求解坡面徑流與大孔隙流耦合方程具有可行性。

3 邊坡水分場參數(shù)敏感性分析

選擇福建省長樂-永定公路S203線某殘積土

路堤邊坡為例， 邊坡高度為6 m， 坡底寬度為12 m， 坡頂、 坡趾寬度均為2 m， 坡趾高度為2 m， 坡比值為1∶2， 主要分析邊坡坡腳剖面。 采用映射的方式劃分網(wǎng)格， 網(wǎng)格尺寸極細(xì)化， 以盡量提高計算精度， 整個邊坡幾何模型劃分網(wǎng)格個數(shù)為5 229， 邊界網(wǎng)格個數(shù)為459。 邊坡幾何模型和網(wǎng)格劃分如圖7所示。

計算邊界條件設(shè)置如下：降雨入滲作用下邊坡幾何模型兩側(cè)及底部為無流量邊界，當(dāng)降雨歷時小于產(chǎn)生徑流時間時，坡面采用流量邊界；當(dāng)降雨歷時大于產(chǎn)生徑流時間時，根據(jù)圖2中的迭代步驟，以積水深度作為壓力水頭邊界條件，初始條件是固定體積含水率為0.1。徑流方程邊界條件為坡頂處流速和積水深度為0，初始條件為開始降雨時坡面上不出現(xiàn)徑流，徑流速率和積水深度為0。

通過定量濾紙法測得體積含水率與吸力之間的關(guān)系，利用MATLAB軟件中Isqcurvefit函數(shù)擬合VG模型參數(shù)，得到土水特征曲線參數(shù)。根據(jù)行業(yè)標(biāo)準(zhǔn)JTG 3430—2020《公路土工試驗(yàn)規(guī)程》［22］，采用變水頭試驗(yàn)測量殘積土基質(zhì)域飽和導(dǎo)水系數(shù)，通過粗砂填制大孔隙域，采用常水頭滲透試驗(yàn)測量飽和導(dǎo)水系數(shù)。采用COMSOL軟件中PDE接口建立邊坡有限元模型，計算參數(shù)如表2所示。

3.1 大孔隙域占比

圖8所示為不同大孔隙域占比ωf時體積含水率隨分析剖面深度的變化。由圖可知：不同ωf時分析剖面一定深度處體積含水率均達(dá)到飽和值，該深度對應(yīng)的區(qū)域?yàn)轱柡蛥^(qū)，飽和區(qū)往下出現(xiàn)過渡區(qū)，在較小深度范圍內(nèi)體積含水率迅速減小至初始值，即達(dá)到濕潤鋒深度。當(dāng)ωf為0.02時，分析剖面深度為0～0.15 m的土壤處于飽和狀態(tài)，該區(qū)間即為飽和區(qū)深度；分析剖面深度為0.15～0.33 m土壤的體積含水率迅速減小，該深度區(qū)段即為過渡區(qū)；分析剖面深度大于0.33 m時體積含水率無響應(yīng)，即濕潤鋒深度為0.33 m。當(dāng)ωf由0.02增至0.05、 0.1、 0.2時，飽和區(qū)深度從0.15 m增至0.33、 0.77、 1.66 m，依次增大120.3%、 133.3%、 115.6%；濕潤鋒深度從0.33m增至0.49、 0.90、 1.81 m，依次增大48.5%、 83.7%、 101.1%。說明隨著ωf的增大，水分入滲深度增大，這是由于大孔隙域?qū)禂?shù)遠(yuǎn)大于基質(zhì)域?qū)禂?shù)，ωf越大，水分快速運(yùn)移的空間也越大，因此能夠在相同時間內(nèi)運(yùn)移至更深處。

圖9所示為積水深度隨大孔隙域占比ωf的變化。由圖可知： 當(dāng)降雨歷時為10 min時， ωf曲線幾乎水平， 此時不同ωf時積水深度基本相同。 當(dāng)降雨歷時為20 min時， 積水深度隨著ωf的增大而略微減小，當(dāng)ωf從0.02增至0.2時， 積水深度從4.515 mm減至4.502 mm， 最大降幅僅為0.289%。當(dāng)降雨歷時為30 min時，ωf取值為0.02、 0.05、 0.1、 0.2對應(yīng)的積水深度為6.382、 6.295、 6.275、 6.216 mm，最大降幅達(dá)2.601%，近乎是降雨歷時為20 min時的10倍。當(dāng)降雨歷時為60 min時， 積水深度變化曲線恢復(fù)水平。 總體上， 隨著ωf的增大， 坡面積水深度呈現(xiàn)下降趨勢， 但是下降幅度隨著降雨歷時的增加而先增大后減小， 在降雨歷時約為30 min時區(qū)分度最顯著。

3.2 大孔隙域與基質(zhì)域?qū)禂?shù)之比

圖10所示為不同大孔隙域與基質(zhì)域?qū)禂?shù)之比u時體積含水率隨分析剖面深度的變化。由圖可知：水分入滲深度隨著u的增大而減小，當(dāng)u為100、 200、 400、 800時，對應(yīng)的飽和區(qū)深度為0.47、 0.43、 0.35、 0.31 m，依次降低8.5%、 18.6%、 11.4%；濕潤鋒深度為0.66、 0.62、 0.55、 0.51 m，依次降低6.1%、 11.3%、 7.3%。這與大孔隙域占比ωf產(chǎn)生的影響相反，主要原因是根據(jù)K=ωmKsm+ωfKsf［17］（其中Ksm、 Ksf為基質(zhì)域、 大孔隙域的飽和導(dǎo)水系數(shù)），當(dāng)u增大時，基質(zhì)域飽和導(dǎo)水系數(shù)減小，體積含水率隨分析剖面深度的變化

致使水分入滲深度減小。當(dāng)u從200增至400時，飽和區(qū)和濕潤鋒深度降幅均最大。

圖11所示為積水深度隨u的變化。 由圖可以看出： 當(dāng)降雨歷時為10、 20 min時， 隨著u的增大， 積水深度無明顯變化。 當(dāng)降雨歷時為30 min時， 積水深度隨著u的增大而增大， 當(dāng)u為100、 200、 400、 800時， 對應(yīng)的積水深度為6.208、 6.234、 6.275、6.289 mm， 最大增長幅度為1.305%。當(dāng)降雨歷時為60 min時，積水深度恢復(fù)至無顯著變化狀態(tài)?？傮w上，降雨歷時為30 min時積水深度區(qū)分度最大。

3.3 大孔隙經(jīng)驗(yàn)參數(shù)

圖12所示為不同大孔隙經(jīng)驗(yàn)參數(shù)rw時體積含水率隨分析剖面深度的變化。由圖可以看出， 當(dāng)rw取值為0.1、 0.2、 0.4、 0.8時， 飽和區(qū)深度為0.97、0.88、 0.76、 0.65 m，" 依次降低8.9%、 14.3%、 14.0%， 濕潤鋒深度為1.13、 1.05、 0.95、 0.84 m， 依次降低7.6%、 9.5%、 11.9%。 這說明入滲深度隨著經(jīng)驗(yàn)參數(shù)rw的增大而減小， 與入滲深度隨u的變化趨勢是一致的。

圖13所示為積水深度隨rw的變化。由圖可知：當(dāng)降雨歷時為10 min時，rw對積水深度幾乎無影響。當(dāng)降雨歷時為20 min時，rw由0.1增至0.8，積水深度由4.482 mm增至4.501 mm，最大增幅為0.424%。當(dāng)降雨歷時為30 min且rw為0.1、 0.2、 0.4、 0.8時，對應(yīng)的積水深度依次為6.205、 6.255、 6.275、 6.319 mm，最大增幅為1.837%。當(dāng)降雨歷時為60 min時，積水深度無顯著變化。與ωf、 u對積水深度影響最大時對應(yīng)的降雨歷時相同，rw對積水深度的影響也是在降雨歷時為30 min時最大。

圖8—13中均給出了無大孔隙即相應(yīng)的大孔隙參數(shù)為0時體積含水率和積水深度值。通過與大孔隙參數(shù)不為0時對比發(fā)現(xiàn)，無大孔隙時，剖面表層土飽和區(qū)深度較小，濕潤鋒深度約為0.15 m，入滲深度遠(yuǎn)小于大孔隙參數(shù)不為0時的入滲深度，說明大孔隙域存在與否對水分入滲的影響很大。除了降雨歷時為10 min，不同大孔隙參數(shù)取值時，相對于無大孔隙時，積水深度稍小，說明無論大孔隙存在與否，坡面積水深度都具有一定差異。

4 結(jié)論

本文中采用2個添加水分交換項(xiàng)的Richards方程，描述土壤大孔隙域和基質(zhì)域水分運(yùn)移過程，結(jié)合運(yùn)動波方程，建立了坡面徑流與大孔隙流耦合模型?；谟邢拊浖﨏OMSOL中的PDE接口，實(shí)現(xiàn)耦合模型的數(shù)值求解，并分析了不同大孔隙參數(shù)對水分場和徑流特性的影響，得到以下主要結(jié)論：

1）水分入滲深度隨著大孔隙域占比ωf的增大而增大，其中ωf從0.1增至0.2時濕潤鋒深度增長幅度最大， 達(dá)到101.1%； 大孔隙域與基質(zhì)域?qū)禂?shù)之比u越大， 水分入滲深度越小， 其中當(dāng)u由200增至400時， 濕潤鋒深度減小幅度最大， 達(dá)到11.3%； 大孔隙經(jīng)驗(yàn)參數(shù)rw越大， 水分入滲深度越小， 其中當(dāng)rw由0.4增至0.8時， 濕潤鋒深度減小幅度最大， 達(dá)到11.9%。

2）積水深度總體上隨著ωf的增大而減小，隨著u、 rw的增大而增大，當(dāng)降雨歷時為30 min時，積水深度區(qū)分度最大。當(dāng)ωf從初始值增至最大值時，積水深度減小2.601%，當(dāng)u、 rw從初始值增至最大值時，積水深度增大1.305%、1.837%。3種大孔隙參數(shù)按照對邊坡水分場影響程度由大到小的順序?yàn)棣豧、 rw、 u。

3）COMSOL軟件的PDE接口可以實(shí)現(xiàn)坡面徑流與大孔隙流耦合模型數(shù)值求解，并且與試驗(yàn)結(jié)果驗(yàn)證良好。相比于無大孔隙結(jié)構(gòu)，考慮大孔隙時水分入滲深度明顯更大，坡面積水深度稍小，無論大孔隙存在與否，邊坡水分場的分布均具有較大差異。

本文中的研究成果為深刻認(rèn)識邊坡力學(xué)場的變化提供了重要依據(jù)，進(jìn)而揭示了邊坡水力交互作用機(jī)制， 可以豐富降雨型大孔隙邊坡失穩(wěn)機(jī)制并為災(zāi)害防治提供技術(shù)支撐。 所建立的大孔隙流與坡面徑流耦合作用通過邊界發(fā)生， 此外， 描述大孔隙流的雙重滲透方程組內(nèi)部還含有耦合項(xiàng)Γw， 因此實(shí)際求解過程十分復(fù)雜。 鑒于本文中主要研究臺風(fēng)強(qiáng)降雨時的水分運(yùn)移問題，" 強(qiáng)降雨作用下坡面短時間內(nèi)即發(fā)生積水， 為了提高計算效率， 本文中忽略了坡面產(chǎn)流臨界時間的計算， 下一步將結(jié)合力學(xué)場的分析實(shí)現(xiàn)更全面的研究。

參考文獻(xiàn)：

［1］ 楊寅， 包紅軍， 彭濤. 臺風(fēng)“鲇魚”強(qiáng)降水引發(fā)的地質(zhì)災(zāi)害氣象風(fēng)險預(yù)警檢驗(yàn)與分析［J］. 暴雨災(zāi)害， 2019， 38（3）： 221.

［2］ 靳紅華， 王林峰， 任青陽， 等. 降雨循環(huán)條件下高切坡穩(wěn)定性演化過程及預(yù)測方法［J］. 土木與環(huán)境工程學(xué)報（中英文）， 2021， 43（4）： 12.

［3］ 許旭堂. 非飽和殘積土坡對降雨入滲的響應(yīng)及其失穩(wěn)演變過程研究［D］. 福州： 福州大學(xué)， 2015.

［4］ 徐則民， 黃潤秋. 山區(qū)流域高蓋度斜坡對極端降雨事件的地下水響應(yīng)［J］. 地球科學(xué)進(jìn)展， 2011， 26（6）： 598.

［5］ 徐宗恒， 徐則民， 曹軍尉， 等. 土壤優(yōu)先流研究現(xiàn)狀與發(fā)展趨勢［J］. 土壤， 2012， 44（6）： 905.

［6］ BOGAARD T A， GRECO R. Landslide hydrology： from hydrology to pore pressure［J］. Wiley Interdisciplinary Reviews： Water， 2016， 3（3）：439.

［7］ JARVIS N， KOESTEL J， LARSBO M. Understanding preferential flow in the vadose zone： recent advances and future prospects［J］. Vadose Zone Journal， 2016， 15（12）： 1.

［8］ 柯云斌， 符元帥， 闕云， 等. 降雨誘發(fā)殘積土陡坡坡面沖刷破壞顆粒流模擬［J］. 濟(jì)南大學(xué)學(xué)報（自然科學(xué)版）， 2019， 33（5）：453.

［9］ 王樂， 田東方. 邊坡滲流與坡面徑流聯(lián)合求解三維有限元模型［J］. 人民珠江， 2019， 40（5）： 24.

［10］ TIAN D F， ZHENG H， LIU D F. A 2D integrated FEM model for surface watergroundwater flow of slopes under rainfall condition［J］. Landslides， 2017， 14（2）： 577.

［11］ 韓同春， 蘇鈺欽， 張宇. 雙層結(jié)構(gòu)邊坡降雨入滲與坡面徑流耦合分析［J］. 工程科學(xué)與技術(shù)， 2020， 52（6）： 145.

［12］ 童富果， 田斌， 劉德富. 改進(jìn)的斜坡降雨入滲與坡面徑流耦合算法研究［J］. 巖土力學(xué)， 2008， 29（4）： 1035.

［13］ SHIGORINA E， RDIGER F， TARTAKOVSKY A M， et al. Multiscale smoothed particle hydrodynamics model development for simulating preferential flow dynamics in fractured porous media［J］. Water Resources Research， 2020， 57（3）： 1.

［14］ SUN X， LUO H， KENICHI S. A coupled thermalhydraulicmechanicalchemical （THMC） model for methane hydrate bearing sediments using COMSOL Multiphysics［J］. Journal of Zhejiang University： Science A， 2018， 19（8）： 600.

［15］ 孟慶娟， 石志飛. 基于周期理論和COMSOL PDE的排樁減振特性研究［J］. 巖土力學(xué)， 2018， 39（11）： 4252.

［16］ 劉俊新， 劉育田， 胡啟軍. 非飽和地表徑流-滲流和流固體耦合條件下降雨入滲對路堤邊坡穩(wěn)定性研究［J］. 巖土力學(xué)， 2010， 31（3）： 904.

［17］ GERKE H H， VAN GENUCHTEN M T. A dualporosity model for simulating the preferential movement of water and solutes in structured porous media［J］. Water Resources Research， 1993， 29（2）： 305.

［18］ 雷志棟， 楊詩秀， 謝森傳. 土壤水動力學(xué)［M］. 北京： 清華大學(xué)出版社， 1988： 4.

［19］ 張培文， 劉德富， 鄭宏， 等. 降雨條件下坡面徑流和入滲耦合的數(shù)值模擬［J］. 巖土力學(xué)， 2004， 25（1）： 111.

［20］ KHNE J M， MOHANTY B P. Water flow processes in a soil column with a cylindrical macropore： experiment and hierarchical modeling［J］. Water Resources Research， 2005， 41（3）： 1.

［21］ NI S M， ZHANG D Q， FENG S Y. Quantitative relationship between hydraulics parameters and soil erosion rate on remolded soil slopes with different textures［J］. Acta Pedologica Sinica， 2019， 56（6）： 1336.

［22］ 交通運(yùn)輸部公路科學(xué)研究院. 公路土工試驗(yàn)規(guī)程： JTG 3430—2020［S］. 北京： 人民交通出版社， 2020： 8.

（責(zé)任編輯：王 耘）