薛 松，童富果，郝 霜，王 軍

（三峽大學(xué)水利與環(huán)境學(xué)院，湖北宜昌 443002）

土體初始飽和度對降雨入滲的影響

薛 松，童富果，郝 霜，王 軍

（三峽大學(xué)水利與環(huán)境學(xué)院，湖北宜昌 443002）

基于水、氣兩相流理論，采用有限單元法計(jì)算了壤土在不同飽和度下的降雨入滲過程，以研究初始飽和度變化對降雨入滲的影響。數(shù)值模擬結(jié)果表明：降雨入滲初期孔隙氣壓會逐漸增大，入滲強(qiáng)度會因孔隙氣的頂托作用逐漸減小并最終趨于穩(wěn)定。穩(wěn)定后的降雨入滲強(qiáng)度主要取決于土體基質(zhì)吸力和滲透特性，與土體初始飽和度大小直接相關(guān)；入滲強(qiáng)度存在最大、最小值，最大值通常會出現(xiàn)在初始飽和度為0或1時(shí)，而最小值所對應(yīng)的飽和度與土體非飽和特性有關(guān)；受土體基質(zhì)吸力和水、氣滲透性的雙重影響，入滲強(qiáng)度會先隨飽和度的增大而逐漸減小，當(dāng)降至最小值后，入滲強(qiáng)度又會隨飽和度的增大而增大。根據(jù)計(jì)算結(jié)果，提出了一個(gè)可計(jì)算相對穩(wěn)定入滲強(qiáng)度的公式。

壤土；兩相流；降雨入滲強(qiáng)度；初始飽和度；有限單元法

降雨入滲是一個(gè)涉及土中水、氣兩相流動的非飽和滲流過程，其對降雨誘發(fā)滑坡、堤防滲透破壞、污染物質(zhì)在水中的遷移等諸多工程問題會產(chǎn)生重要影響。降雨入滲通常與土體非飽和特性、土體初始飽和度、降雨強(qiáng)度等諸多因素有關(guān)［1-2］，其中土體初始飽和度是影響降雨入滲的關(guān)鍵因素之一［3］。

人們在研究非飽和滲流時(shí)通常采用試驗(yàn)方法［4-6］和數(shù)值方法［7-12］。試驗(yàn)方法多用于對入滲規(guī)律的定性研究；在定量研究方面，由于對多物理量控制困難，難以準(zhǔn)確分析單因素變化對入滲的影響。數(shù)值方法便于多物理量控制，適于研究單一變量與入滲的相關(guān)關(guān)系。根據(jù)降雨入滲數(shù)學(xué)模型對非飽和滲流對象的描述，非飽和滲流理論可分為單相流理論和多相流理論。傳統(tǒng)計(jì)算軟件和算法多數(shù)建立在描述水分單相流動的Richards方程上［7］，計(jì)算時(shí)入滲流量邊界需人為給定，而在降雨入滲的計(jì)算中入滲水量又是反映非飽和滲流規(guī)律的關(guān)鍵指標(biāo)，因此單相流理論存在一定的不足。此外，Richards方程假定孔隙氣壓力恒定，未考慮水、氣耦合的影響，忽略了氣相的作用。非飽和滲流涉及土中水、氣兩相的流動［13-14］，兩相流理論考慮了水、氣的相互影響［15-17］，能更準(zhǔn)確地反映土中水氣的運(yùn)動規(guī)律?；趦上嗔黢詈系姆匠讨恍杞o定降雨入滲的氣壓力邊界，降雨入滲流量通過水、氣耦合求解，因此更符合實(shí)際情況。以往對兩相流的研究大多集中在土中氣壓力變化對入滲的影響上，而在初始飽和度對降雨入滲的影響方面研究相對較少。

鑒于此，本文基于水、氣兩相流全耦合的微分控制方程，考慮土中水、氣的相互耦合作用，采用有限單元法計(jì)算不同初始飽和度下一維均質(zhì)土柱的降雨入滲過程，探討入滲過程中的水、氣耦合機(jī)理，分析初始飽和度變化對降雨入滲強(qiáng)度、飽和度沿深度分布的影響。

1 水、氣兩相流控制方程

基于連續(xù)介質(zhì)理論，將土體視為由固相、液相和氣相組成的多孔連續(xù)介質(zhì)，土中水、氣的運(yùn)動隨時(shí)間和空間不斷變化，可描述為關(guān)于時(shí)間和空間的偏微分方程組［18］。水在土體中的流動主要源于水壓力驅(qū)動，并受土體孔隙特征、含水狀態(tài)等因素影響?；谝合噘|(zhì)量守恒可導(dǎo)出土體內(nèi)水流動的控制微分方程：

式中：φ為土體孔隙率；Sr為水飽和度；μw為液相黏滯系數(shù)，Pa·s；krw為水相相對滲透系數(shù)；k為由孔隙特征決定的土體本征滲透系數(shù)，m2；pw為孔隙水壓力，Pa；Qw為液相內(nèi)源項(xiàng)，kg/（m3·s）；ρw為水相密度，kg/m3；g為重力加速度矢量。

氣相在土體中的流動主要由氣壓力驅(qū)動，并受土體孔隙特征、含水狀態(tài)、氣相黏滯系數(shù)等因素影響?；跉庀噘|(zhì)量守恒可導(dǎo)出土體內(nèi)氣流動的控制微分方程：

式中：μg為氣相黏滯系數(shù)，Pa·s；krg為氣相相對滲透系數(shù)；Qg為氣相內(nèi)源項(xiàng)，kg/（m3·s）；pg為孔隙氣壓力，Pa；ρg為氣相密度，kg/m3。

基質(zhì)吸力pc是水氣交界面承受氣壓力和水壓力不相等所產(chǎn)生的不平衡力，根據(jù)其定義［18］，可表達(dá)為

方程（1）（2）構(gòu)成基于時(shí)間和空間的非線性偏微分方程組，求解時(shí)需對其進(jìn)行時(shí)間和空間離散。本文空間離散采用Garlerkin有限單元法，時(shí)間離散采用差分法。以孔隙氣壓力pg、水飽和度Sr為未知量，通過對水、氣兩個(gè)方程的循環(huán)迭代計(jì)算實(shí)現(xiàn)每個(gè)時(shí)間步內(nèi)非線性方程組的求解。

2 計(jì)算模型

2.1 幾何模型

為計(jì)算通氣條件下邊坡的水分入滲過程，假設(shè)模擬區(qū)域?yàn)橐粶\層均質(zhì)各向同性土體，土柱高3m。根據(jù)模擬對象建立有限元模型，計(jì)算網(wǎng)格節(jié)點(diǎn)數(shù)40個(gè)，計(jì)算單元數(shù)20個(gè)，單元最小厚度0.1m。

為研究土體的入滲強(qiáng)度與初始飽和度相關(guān)關(guān)系，假定降雨強(qiáng)度遠(yuǎn)大于入滲強(qiáng)度，地表始終處于飽和狀態(tài)，地表徑流對入滲的影響相對較小，可忽略不計(jì)。土柱整體以同種初始飽和度起算，地表與底部均為大氣壓邊界，側(cè)邊為不透水不透氣邊界。

2.2 主要本構(gòu)關(guān)系及參數(shù)

由于方程（1）（2）求解時(shí)，具有5個(gè)未知參數(shù)Sr、pw、pg、krw、krg，還需引入3個(gè)本構(gòu)模型關(guān)系，即土水特征曲線、水相相對滲透系數(shù)曲線、氣相相對滲透系數(shù)曲線。本文選取長江流域最為常見的壤土作為計(jì)算對象分析其降雨入滲規(guī)律，通過分析相關(guān)文獻(xiàn)［19-24］得到其具體模型及參數(shù)如下：

a.土水特征曲線。反映土體基質(zhì)吸力大小與飽和度之間的關(guān)系（圖1），本文選取應(yīng)用較為廣泛且?guī)缀踹m用于所有土質(zhì)類型的van Genuchten模型［20］：

式中：Srw為殘余水飽和度；P0、m為與材料特性有關(guān)的參數(shù)。根據(jù)相關(guān)資料，確定壤土計(jì)算參數(shù)為P0= 0.278、m=0.359、Srw=0.1。

圖1 土體非飽和滲透特性曲線

b.水相相對滲透系數(shù)曲線。描述土體內(nèi)部水相相對流動性與飽和度關(guān)系的曲線（圖1），本文選取最為常用的van Genuchten-Mualem模型［20-21］：式中n為與材料特性有關(guān)的參數(shù)，本文計(jì)算取0.9。

c.氣相相對滲透系數(shù)曲線。描述土體內(nèi)部氣相相對流動性與飽和度關(guān)系的曲線，氣體的相對滲透系數(shù)采用Brooks-Corey模型［22］：

式中Srg為殘余氣飽和度。

根據(jù)方程（4）（5）中參數(shù)得到土體的土水特征曲線和水相相對滲透系數(shù)曲線，亦即圖1所示的土體的非飽和滲透特性曲線。方程（1）（2）中涉及的其他主要參數(shù)取值：ρw=1 t/m3，ρg=1.29 kg/m3，g= 9.8m/s2，φ=0.4，μw=1.0×10-3Pa·s，μg=1.84× 10-5Pa·s，k=2.88×10-13m2。

3 降雨入滲與初始飽和度相關(guān)性分析

降雨入滲是一個(gè)涉及水、氣兩相流動的過程，入滲過程受到土中水、氣兩相的相互作用及土體的非飽和滲透特性等多重因素的影響，而初始飽和度的改變對土體的水、氣含量和土體的非飽和特性均具有一定的影響。本節(jié)先分析降雨入滲過程中水、氣的耦合作用，而后研究降雨入滲強(qiáng)度與土體初始飽和度的相關(guān)關(guān)系。為便于分析比較，采用相對入滲強(qiáng)度來衡量降雨入滲的強(qiáng)弱，相對入滲強(qiáng)度為無量綱變量，定義為雨水入滲速率與土體飽和滲透系數(shù)之比。

3.1 降雨入滲中的水、氣耦合作用

降雨入滲是水、氣兩相相互影響的耦合過程，孔隙氣對雨水入滲的影響明顯。非飽和狀態(tài)下入滲強(qiáng)度的大小主要取決于水相相對滲透系數(shù)和孔隙水壓力梯度，而孔隙水壓力梯度受到基質(zhì)吸力、孔隙氣壓力及重力的影響。

降雨入滲初期，表層土體飽和度較小，基質(zhì)吸力較大，故水壓力梯度較大，入滲強(qiáng)度較大。此時(shí)表層土體氣相相對滲透系數(shù)較大，氣體流動性較好，孔隙氣可從表層溢出，壓力接近大氣壓力，因而氣體對入滲的影響可忽略。隨著表層土體逐漸達(dá)到飽和狀態(tài)，氣相相對滲透系數(shù)迅速減小并趨近于0，氣體難以從土體表層溢出，孔隙氣因受水分入滲的擠壓，壓力迅速上升并產(chǎn)生與入滲方向相反的頂托力，致使孔隙水壓力梯度迅速下降，故入滲強(qiáng)度迅速下降（圖2）。隨著孔隙氣壓力逐漸趨于穩(wěn)定，氣體對重力梯度的頂托也趨于穩(wěn)定，此時(shí)入滲主要受到濕潤鋒區(qū)域基質(zhì)吸力及水相滲透系數(shù)的共同作用，趨于相對穩(wěn)定的入滲狀態(tài)。

表層土體飽和后，受表層通氣條件改變和水分下滲擠壓作用的影響，孔隙氣壓力迅速上升，進(jìn)而導(dǎo)致入滲強(qiáng)度迅速下降。由于前期入滲強(qiáng)度較大，受擠壓而產(chǎn)生的孔隙氣體來不及流動，因此孔隙氣壓力會出現(xiàn)極大值（圖3）。隨著入滲強(qiáng)度的減小和土體內(nèi)部孔隙氣體的流動，孔隙氣壓力會緩慢下降并逐漸趨于穩(wěn)定，此時(shí)入滲強(qiáng)度也趨于穩(wěn)定。此外，隨著濕潤鋒的推移，土體內(nèi)部飽和度增大，氣相相對滲透系數(shù)減小，孔隙氣受擠壓作用更明顯，因而隨著飽和度的增大孔隙氣壓力呈增大趨勢。如土體深度0.6m處，孔隙氣壓力隨著雨水的入滲，而存在緩慢上升的趨勢。

圖2 降雨入滲強(qiáng)度隨時(shí)間變化曲線

圖3 孔隙氣壓力隨時(shí)間變化曲線（Sr=0.5）

土體表層飽和后，孔隙氣體難以溢出，表層土體因部分孔隙被氣體占據(jù)而難以達(dá)到完全飽和狀態(tài)，其飽和度大小與初始飽和度關(guān)系密切（圖4）。初始飽和度較大時(shí)，土體的氣相相對滲透系數(shù)較小，氣壓力梯度較大，孔隙氣壓力較大；而基質(zhì)吸力的大小與孔隙氣壓力直接相關(guān)，孔隙氣壓力越大，基質(zhì)吸力越大；受基質(zhì)吸力較大的影響（圖1），土體所能達(dá)到的飽和度相對較小。因此初始飽和度為0.6時(shí)的表層含水量明顯小于初始飽和度為0.15時(shí)的含水量，故初始飽和度越大，表層土體的飽和度越小，孔隙氣壓力對水分入滲影響越強(qiáng)。但當(dāng)土體初始飽和度接近完全飽和時(shí)，孔隙氣體含量極少，水分主要通過孔隙通道下滲，對氣體擠壓作用有限，孔隙氣壓力對水分入滲影響極弱。

隨著初始飽和度的增大，土體單位入滲斷面內(nèi)，氣體含量越少，水分入滲通道越多，水相相對滲透系數(shù)越大（圖1）。受基質(zhì)吸力與水相滲透性的影響，飽和度沿深度分布會因初始飽和度的不同而存在差異（圖4），初始飽和度較小時(shí)，土體的基質(zhì)吸力較大，單位區(qū)域內(nèi)水分的流入強(qiáng)度較大，而水相相對滲透系數(shù)較小，流出強(qiáng)度較小，受流入流出強(qiáng)度差的影響會呈現(xiàn)出較陡的濕潤鋒。飽和度較大時(shí)，基質(zhì)吸力較小，而水相相對滲透系數(shù)較大，入滲水量主要通過水分通道下滲，流入流出強(qiáng)度相差較小，故濕潤鋒較緩。

3.2 初始飽和度對降雨入滲的影響

降雨入滲初期，入滲強(qiáng)度變化劇烈，為便于比較分析，建立降雨入滲強(qiáng)度與初始飽和度的相關(guān)關(guān)系，本文主要研究穩(wěn)定后的入滲強(qiáng)度。入滲穩(wěn)定后，孔隙氣壓力趨于穩(wěn)定，其對重力梯度的頂托作用已經(jīng)完成，此時(shí)降雨入滲強(qiáng)度主要取決于基質(zhì)吸力和水相、氣相滲透特性，重力梯度因受氣壓力的頂托可忽略不計(jì)。

降雨初期，土體的入滲強(qiáng)度主要取決于土體表層基質(zhì)吸力，基質(zhì)吸力越大，孔隙水壓力梯度越大，故初始飽和度越小，初期降雨入滲強(qiáng)度越大（圖2）。隨著入滲逐漸趨于穩(wěn)定，入滲強(qiáng)度主要由濕潤鋒處的基質(zhì)吸力和水相相對滲透系數(shù)共同決定。入滲強(qiáng)度與初始飽和度呈非線性關(guān)系，穩(wěn)定入滲強(qiáng)度隨初始飽和度的增大先減小后增大，穩(wěn)定入滲強(qiáng)度的最大值會出現(xiàn)在殘余飽和度時(shí)（有效飽和度為0）或完全飽和狀態(tài)（有效飽和度為1），最小值會出現(xiàn)在非飽和狀態(tài)。初始飽和度較小時(shí)，土體的基質(zhì)吸力較大，而水相相對滲透系數(shù)較小，水分的入滲主要由基質(zhì)吸力驅(qū)動，故隨著初始飽和度的增大，基質(zhì)吸力減小，穩(wěn)定入滲強(qiáng)度減小。隨著初始飽和度的增大，雖然基質(zhì)吸力繼續(xù)減小，但水相相對滲透系數(shù)迅速增大（圖1），入滲強(qiáng)度會因基質(zhì)吸力與水相滲透性的雙重影響出現(xiàn)最小值（圖5）。入滲強(qiáng)度達(dá)到最小值后，隨著飽和度的增大，基質(zhì)吸力對入滲的影響繼續(xù)減小，水相相對滲透系數(shù)趨于1，水分主要通過孔隙通道下滲，此時(shí)入滲強(qiáng)度主要受到水相相對滲透系數(shù)的影響，故入滲強(qiáng)度隨初始飽和度的增大而增大。以壤土為例計(jì)算表明，壤土的穩(wěn)定入滲強(qiáng)度的最大值出現(xiàn)在接近殘余飽和度時(shí)，而最小值出現(xiàn)在飽和度為0.8～0.9時(shí)。

圖4 不同時(shí)刻土體飽和度沿深度分布

圖5 相對穩(wěn)定入滲強(qiáng)度與初始飽和度關(guān)系曲線

基于應(yīng)用的需要，對本文計(jì)算所得的壤土的降雨入滲關(guān)系數(shù)據(jù)進(jìn)行了簡單的回歸擬合，分析表明3次多項(xiàng)式能夠較好地?cái)M合出初始飽和度與穩(wěn)定入滲強(qiáng)度的相關(guān)關(guān)系：

式中：V（Sr0）為相對穩(wěn)定入滲強(qiáng)度；Sr0為初始水飽和度。

4 結(jié) 論

a.降雨導(dǎo)致表層土體趨近于飽和后，孔隙氣體受水分入滲的擠壓作用，會產(chǎn)生反向的頂托力，降雨入滲強(qiáng)度會因氣體頂托逐漸減小并趨于穩(wěn)定，穩(wěn)定后的入滲強(qiáng)度主要受基質(zhì)吸力與水相相對滲透系數(shù)的影響。飽和度較小時(shí)，穩(wěn)定入滲強(qiáng)度受基質(zhì)吸力的影響較大，隨飽和度的增大而減小，接近飽和時(shí)，穩(wěn)定入滲強(qiáng)度受水相相對滲透系數(shù)的影響較大，隨飽和度的增大而增大。

b.穩(wěn)定入滲強(qiáng)度的最大值出現(xiàn)在飽和度最小或完全飽和狀態(tài)，最小值與土體非飽和滲透特性相關(guān)，受土水特征關(guān)系和水相、氣相滲透特性的綜合影響。根據(jù)計(jì)算所得的相對穩(wěn)定入滲強(qiáng)度與飽和度的相關(guān)關(guān)系，提出了一個(gè)能估算相對穩(wěn)定入滲強(qiáng)度的入滲公式，計(jì)算表明壤土穩(wěn)定入滲強(qiáng)度的最大值出現(xiàn)在接近殘余飽和度時(shí)，最小值出現(xiàn)在飽和度為0.8～0.9時(shí)。

c.入滲過程中，由于孔隙氣的存在，土體表層難以達(dá)到完全飽和狀態(tài)，初始飽和度越大，孔隙氣壓越大，氣體對土體飽和度的影響越強(qiáng)，土體表層所能達(dá)到的飽和度相對越小。此外，受水相滲透性的影響，土體含氣率大時(shí)，水分主要依靠基質(zhì)吸力驅(qū)動下滲；含氣率小時(shí)，水分主要通過孔隙通道下滲。

［1］朱偉，陳學(xué)東，鐘小春.降雨入滲規(guī)律的實(shí)測與分析［J］.巖土力學(xué)，2006，27（11）：1873-1879.（ZHU Wei，CHEN Xuedong，ZHONG Xiaochun.Observation and analysis of rainfall infiltration［J］.Rock and Soil Mechanics，2006，27（11）：1873-1879.（in Chinese））

［2］呂剛，吳祥云.土壤入滲特性影響因素研究綜述［J］.中國農(nóng)學(xué)通報(bào)，2008，24（7）：494-499.（L¨U Gang，WU Xiangyun.Review on influential factors of soil infiltration characteristics［J］.Chinese Agricultural Science Bulletin，2008，24（7）：494-499.（in Chinese））

［3］朱偉，程南軍，陳學(xué)東，等.淺談非飽和滲流的幾個(gè)基本問題［J］.巖土工程學(xué)報(bào)，2006，28（2）：235-240.（ZHU Wei，CHENG Nanjun，CHEN Xuedong，et al.Some fundamental problems of unsaturated seepage［J］.Chinese Journal of Geotechnical Engineering，2006，28（2）：235-240.（in Chinese））

［4］劉目興，聶艷，于婧.不同初始含水率下黏質(zhì)土壤的入滲過程［J］.生態(tài)學(xué)報(bào)，2012，32（3）：871-878.（LIU Muxing，NIE Yan，YU Jing.The infiltration process of clay under different initial soil water contents［J］.Acta Ecologica Sinca，2012，32（3）：871-878.（in Chinese））

［5］曾辰，王全九，樊軍.初始含水率對土壤垂直線源入滲特征的影響［J］.農(nóng)業(yè)工程學(xué)報(bào)，2010，26（1）：24-30.（ZENG Chen，WANG Quanjiu，F(xiàn)AN Jun.Effect of initial water content on vertical line-source infiltration characteristics of soil［J］.Transactions of the CSAE，2010，26（1）：24-30.（in Chinese））

［6］張振華，謝恒星，劉繼龍，等.氣相阻力與土壤容重對一維垂直入滲影響的定量分析［J］.水土保持學(xué)報(bào)，2005，19（4）：36-39.（ZHANG Zhenhua，XIE Hengxing，LIU Jilong，et al.Quantitative analysis of influence of air resistance and soil bulk density on one-dimensional vertical infiltration［J］.Journal of Soil and Water Conservation，2005，19（4）：36-39.（in Chinese））

［7］雷志棟，楊詩秀.非飽和土壤水一維流動的數(shù)值計(jì)算［J］.土壤學(xué)報(bào)，1982，19（2）：141-153.（LEI Zhidong，YANG Shixiu.Numericalmethod of one dimensional flow through unsaturated soils［J］.Acta Pedologica Sinica，1982，19（2）：141-153.（in Chinese））

［8］范嚴(yán)偉，黃寧，馬孝義.層狀土垂直一維入滲土壤水分運(yùn)動數(shù)值模擬與驗(yàn)證［J］.水土保持通報(bào)，2015，35（1）：215-219.（FAN Yanwei，HUANG Ning，MA Xiaoyi. Numerical simulation and verification of layered soil water movement in vertical one-dimensional infiltration［J］. Bulletin of Soil and Water Conservation，2015，35（1）：215-219.（in Chinese））

［9］張旭明，孫建國.光滑有限單元法及其應(yīng)用［J］.河海大學(xué)學(xué)報(bào)（自然科學(xué)版），2012，34（5）：582-585.（ZHANG Xuming，SUN Jianguo.Smoothed finite element method and its application［J］.Journal of Hohai University（Natural Sciences），2012，34（5）：582-585.（in Chinese））

［10］介玉新，秦曉艷，金鑫，等.加筋土高邊坡的有限元分析［J］.水利水電科技進(jìn)展，2012，32（1）：83-86.（JIE Yuxin，QIN Xiaoyan，JIN Xin，et al.FEM analysis of high reinforced soil slopes［J］.Advances in Science and Technology of Water Resources，2012，32（1）：83-86.（in Chinese））

［11］張曉悅，王棟，張曉樂，等.基于水-氣二相流模型的土坡穩(wěn)定性分析［J］.水利水電科技進(jìn)展，2012，32（2）：23-27.（ZHANG Xiaoyue，WANG Dong，ZHANG Xiaole，et al.Stability analysis of soil slope based on water-air two-phase flow model［J］.Advances in Science and Technology of Water Resources，2012，32（2）：23-27.（in Chinese））

［12］解河海，馮杰，馮青，等.考慮大孔隙的土壤入滲模型［J］.水利水電科技進(jìn)展，2011，31（5）：35-38.（XIE Hehai，F(xiàn)ENG Jie，F(xiàn)ENG Qing，et al.Infiltration model for soils considering macrospores［J］.Advances in Science and Technology ofWater Resources，2011，31（5）：35-38.（in Chinese））

［13］彭勝，陳家軍，王金生，等.包氣帶水氣二相流國外研究綜述［J］.水科學(xué)進(jìn)展，2000，11（3）：333-338.（PENG Sheng，CHEN Jiajun，WANG Jinsheng，et al.Review of foreign study on two-phase flow in vadose zone［J］. Advances in Water Science，2000，11（3）：333-338.（in Chinese））

［14］孫冬梅，朱岳明，張明進(jìn).降雨入滲過程的水-氣二相流模型研究［J］.水利學(xué)報(bào)，2007，38（2）：150-156.（SUN Dongmei，ZHU Yueming，ZHANG Mingjin.Water-air twophase flow model for numerical analysis of rainfall infiltration［J］.Journal of Hydraulic Engineering，2007，38（2）：150-156.（in Chinese））

［15］SUN D，ZANG Y，SEMPRICH S.Effects of airflow induced by rainfall infiltration on unsaturated soil slope stability［J］.Transport in Porous Media，2015，107（3）：821-841.

［16］SZYMKIEWICZA，NEUWEILER I，HELMIG R.Influence of heterogeneous air entry pressure on large scale unsaturated flow in porous media［J］.Acta Geophys，2014，62（5）：1179-1191.

［17］KUANG X，JIAO J J，LI H.Review on airflow in unsaturated zones induced by natural forcings［J］.Water Resources Research，2013，49（10）：6137-6165.

［18］TONG F，NIEMIA，YANG Z，et al.A numericalmodel of tracer transport in a non-isothermal two-phase flow system for CO2geological storage characterization［J］.Transport in Porous Media，2013，98（1）：173-192.

［19］TONG F，JING L，TIAN B.A water retention curvemodel for the simulation of coupled thermo-hydro-mechanical processes in geological porous media［J］.Transport in Porous Media，2012，91（2）：509-530.

［20］van GENUCHTEN M T.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils［J］.Soil Science Society of America Journal，1980，44（5）：892-898.

［21］MUALEM Y.A new model for predicting the hydraulic conductivity of unsaturated porous media［J］.Water Resources Research，1976，12（3）：513-522.

［22］COREY A T.The interrelation between gas and oil relative permeabilities［J］.Prod Mon，1954，19：38-41.

［23］查元源，周發(fā)超，楊金忠.一種由土壤剖面含水率估算土壤水力參數(shù)的方法［J］.水利學(xué)報(bào)，2011，42（8）：883-891.（ZHA Yuanyuan，ZHOU Fachao，YANG Jinzhong. An approach to estimate soil hydraulic parameters using soil profile moisture data［J］.Journal of Hydraulic Engineering，2011，42（8）：883-891.（in Chinese））

［24］范嚴(yán)偉，趙文舉，冀宏.垂直一維入滲土壤水分分布與入滲特性數(shù)值模擬［J］.蘭州理工大學(xué)學(xué)報(bào)，2012，38（2）：51-55.（FAN Yanwei，ZHAO Wenju，JI Hong. Numerical simulation of soil moisture distribution and infiltration characteristics of one-dimensional vertical infiltration［J］.Journal of Lanzhou University of Technology，2012，38（2）：51-55.（in Chinese））

Effects of initial degree of saturation of soil on rainfall infiltration

XUE Song，TONG Fuguo，HAO Shuang，WANG Jun
（College ofHydraulic&Environmental Engineering，China Three Gorges University，Yichang 443002，China）

Based on the theory of gas-liquid two-phase flow，the finite element method was used to simulate rainfall infiltration in loam soil with different degrees of saturation to investigate the effects of the initial degree of saturation on rainfall infiltration.Numerical results show that the pore air pressure increases gradually in the initial stage of rainfall infiltration；then，the infiltration rate decreases to a stable value because of the jacking force caused by pore air.The stable infiltration rate ismainly determined by thematrix suction and permeability characteristics of soil，varying with the initial degree of saturation.With the variation of the initial degree of saturation，the infiltration rate reaches its maximum or minimum values.The maximum infiltration rate usually appears in the stage of residual or full saturation，while the minimum infiltration rate occurs in the unsaturated stage，depending on soil unsaturated characteristics.Affected by the matrix suction of soil and the relative permeability ofwater and air，the infiltration rate gradually decreaseswith the increase of the degree of saturation at first.When the infiltration rate drops to theminimum，it starts to rise with the increase of the degree of saturation.According to the results of numerical simulation，a formula is proposed for calculation of the relatively stable infiltration rate for different initial degrees of saturation.

loam soil；two-phase flow；rainfall infiltration；initial degree of saturation；finite elementmethod

TV139.11

A

1006-7647（2016）04-0031-05

10.3880/j.issn.1006-7647.2016.04.006

2015- 11 21 編輯：熊水斌）

國家自然科學(xué)基金（51279090）

薛松（1992—），男，碩士研究生，主要從事基于水氣兩相流的降雨入滲研究。E-mail：125930515@qq.com

童富果（1972—），男，教授，博士，主要從事水工結(jié)構(gòu)和巖土工程多場耦合問題及其數(shù)值方法研究。E-mail：tfg@ctgu.edu.cn