劉剛+童富果+習(xí)念念+郝霜

摘要：基于多相流及多場耦合理論，采用有限單元法計算了三峽庫區(qū)木魚包滑坡體溫度隨氣溫變化情況，進而依據(jù)溫度與水體黏滯性的相關(guān)關(guān)系分析氣溫變化對坡體滲流的影響。計算表明，坡內(nèi)溫度受氣溫影響僅限于坡體淺層區(qū)域，距地表深度越深，影響程度越小，在時間上的滯后性越大。溫度變化可改變水的黏滯性，進而引起土體滲透參數(shù)的改變，計算表明因氣溫變化引起的坡體滲透性改變較為顯著。坡體滲透系數(shù)受氣溫變化的影響隨深度呈非線性分布，距地表越深受氣溫變化的影響越弱，通常大于雨水最大入滲深度以后，坡體滲透系數(shù)不再發(fā)生顯著變化。另外，坡體滲透系數(shù)隨氣溫呈現(xiàn)季節(jié)性變化規(guī)律，夏季氣溫高于坡內(nèi)土體溫度，滲透系數(shù)隨深度的增加而減??；冬季坡體內(nèi)溫度高于氣溫，滲透系數(shù)隨深度的增加而增加。

關(guān)鍵詞：氣溫；滲透系數(shù)；坡體滲流；多相流；季節(jié)性

中圖分類號：TV223.4 文獻標識碼：A 文章編號：1672-1683（2017）03-0158-06

Abstract：Based on the multi-phase flow & multi-field coupling theory，this paper adopted the finite element method to calculate the temperature distribution of Muyubao landslide as the atmospheric temperature changes.Further，it analyzed the influence of atmospheric temperature changes on the slope seepage according to the correlation between temperature and viscosity of water.The results of calculation showed that the effect of atmospheric temperature on temperature of the slope body is limited to the shallow surface；as the depth increases，the influence decreases and the hysteresis increases concurrently.The mobility of the water is affected by the change of the slope body′s temperature and it can change the water permeability in the inner slope.The effect of atmospheric temperature on permeability coefficient of the slope has a nonlinear relationship with depth.The effect of atmospheric temperature on the permeability coefficient declines as the depth increases.Generally，the permeability coefficient of the slope no longer change significantly when the depth is deeper than the maximum depth of infiltration.The permeability coefficient of the slope shows seasonal variation with the atmospheric temperature.The permeability coefficient decreases with the increase of depth when the atmospheric temperature is higher than the soil in summer.In contrast，the permeability coefficient increases with the increase of depth in winter.

Key words：air temperature；permeability coefficient；slope seepage；multiphase flow；seasonal

坡體滲流通常屬于非飽和滲流的范疇，涉及到土（石）、水、氣組成的多相物質(zhì)的耦合作用[1]。坡體滲流通常與坡內(nèi)水體滲透性、坡體初始含水率、滲透邊界等諸多因素有關(guān)[2-3]，其中水體滲透性是影響坡體滲流的關(guān)鍵因素之一。水在土體中的滲透性主要受水的黏滯性、相對飽和度以及土體本征滲透性的影響[4]，而氣溫變化對土體本征滲透性以及相對飽和度的影響較小，因此氣溫變化主要通過改變水的黏滯性來影響坡內(nèi)水體滲透性。在坡體溫度場方面，許多學(xué)者做了大量研究并取得了豐碩的成果[2，5-10]。然而，較多研究者在計算坡體溫度場時僅考慮了熱量在土體中的傳導(dǎo)，未考慮水、熱耦合的熱傳輸作用，導(dǎo)致坡體溫度場計算結(jié)果不準確[11-12]。此外，部分研究者未能基于氣溫引起的坡體溫度改變結(jié)果進一步計算其對坡內(nèi)水體滲透性的影響[13-14]，常因溫度變化不顯著而得出氣溫變化對坡體滲流無影響的結(jié)論。本文基于多相流及多場耦合的理論和方法，完整地考慮了熱傳導(dǎo)、熱傳輸及水熱耦合效應(yīng)，準確模擬坡體溫度場分布，定量分析氣溫變化對坡內(nèi)水體滲透性的影響。

1 計算模型

1.1 水、氣、熱傳輸控制微分方程

依據(jù)連續(xù)介質(zhì)力學(xué)的基本理論及方法，將土體視為是由氣、液、固三相介質(zhì)組成的連續(xù)體，其內(nèi)部水、氣、熱的傳輸在數(shù)學(xué)上可表述為基于時間和空間的偏微分方程描述。水（液相）在土體中的流動主要源于水壓力梯度驅(qū)動，并受土體孔隙率、土體含水量以及水的黏滯性等因素影響，基于液態(tài)水的質(zhì)量守恒可導(dǎo)出土體內(nèi)水流動的控制微分方程[14]：

上述三個方程是基于時間和空間的非線性偏微分方程組，求解時需對其進行時間和空間離散。本文空間離散采用Garlerkin有限單元法，時間離散采用差分法。以孔隙氣壓Pg、水飽和度Sl以及溫度T為未知量，通過對水、氣、熱三個方程的循環(huán)迭代計算實現(xiàn)每個時間步內(nèi)非線性方程組的求解。

1.2 本構(gòu)模型及參數(shù)

計算所涉及的本構(gòu)關(guān)系主要包括土水特征曲線、水的相對滲透系數(shù)、氣的相對滲透系數(shù)以及土體熱傳導(dǎo)系數(shù)等，具體模型及其參數(shù)表述如下。

（1）土水特征曲線。土水特征曲線是用來描述非飽和土體強度與滲透特性的關(guān)系曲線，本文采用常用的Van Genuchten模型[16-17]，可表述為基質(zhì)吸力與飽和度之間的函數(shù)關(guān)系：

（4）水、氣粘黏滯系數(shù)。水、氣黏滯性的大小與溫度、壓強等因素有關(guān)，在邊坡降雨過程中，黏滯系數(shù)隨壓強變化不大，溫度是影響?zhàn)禂?shù)大小的主要因素，由Guvanasen and Chan模型[21-22]可得水、氣黏滯系數(shù)為

式中：λs、λw、λg分別為固、液、氣三相各自的熱傳導(dǎo)系數(shù)，其值分別為λs=3.20，λw，=0.58，λg=0.024，（單位：W/（m·℃）。η1、η2為模型參數(shù)，其中η1可表示為孔隙率的函數(shù)，其值可表示為η1=0.0692n-0.7831，η2可由孔隙率及飽和度表示出來，η2=0.59S1.487n-0.0404l。

此外，計算涉及到的主要參數(shù)為孔隙率n=0.35，土體本征滲透性k=4.0×10-12m2，重力加速度g=9.8 m/s2，氣體密度ρg=1.29 kg/m3，液相密度ρl=1×103 kg/m3，固體、液體、氣體的比熱容分別為cs=2 080，cl=4 183，cg=1 004（單位：J/（kg·℃））。

1.3 幾何模型

木魚包滑坡位于長江右岸，距三峽大壩壩址56 km，地屬湖北省秭歸縣沙鎮(zhèn)溪鎮(zhèn)范家坪村。木魚包滑坡前緣高程135 m，寬約1 450 m，剪出口位于滑坡后緣北段，最高標高520 m，東側(cè)大樂溝大體近南北向延伸，長約620 m，西側(cè)邊界基本沿鵝卵石溝延伸，約1 200 m，后緣滑壁平直光滑，長約數(shù)百米。木魚包滑坡體表層為松散堆積層，下層為擾動破壞的層狀石英砂巖巖體。選取木魚包滑坡體主滑動剖面為研究對象，采用二維平面滲流模型，沿坡體縱向取1 500 m，豎直向取600 m。對坡面進行有限元網(wǎng)格剖分時，為更好的模擬坡表降雨入滲，坡體表層區(qū)域采用較小網(wǎng)格，最小網(wǎng)格厚度約為0.25 m。最終計算網(wǎng)格共計2 033個單元和1 971個節(jié)點，見圖1。

1.4 邊界及初始條件

邊界條件可分為滲流邊界條件和溫度邊界條件。對于水氣流動，滑坡體后緣側(cè)面及底部均為不透水邊界，其邊界水壓力及其變化可由計算得到。坡體前緣側(cè)面位于庫水位以下，為已知水頭邊界。水面以上為降雨入滲邊界，降雨強度據(jù)三峽庫區(qū)多年平均降雨資料得出（圖2）。對傳熱過程而言，滑坡體后緣側(cè)面設(shè)為不透熱邊界，底部邊界受氣溫影響較小，溫度相對穩(wěn)定，計算中設(shè)定為20 ℃。坡表為已知溫度邊界，其中庫水位以下部分的溫度等于水溫，水位以上部分的溫度等于氣溫，并隨氣溫變化（圖3）。

滲流及溫度初始條件涉及到飽和度和溫度的數(shù)值大小及其空間分布，通常難以準確給出。為獲取較為合理的滲流及溫度初始條件，本分析在等溫場、等飽和度條件下，結(jié)合三峽庫區(qū)多年實測降雨及氣溫情況進行長時間模擬計算，以計算所得滲流場及溫度場作為本研究的初始條件。2 氣溫與坡體滲流的相關(guān)性分析

坡內(nèi)水體滲透性主要由土體孔隙大小及其特性決定，同時也受土體變形、土壤含水量、水的黏滯性等因素的影響。溫度變化對坡內(nèi)土體孔隙特性影響較小，而對坡內(nèi)水體的黏滯性影響較為顯著，因此，本節(jié)先分析氣溫變化對坡體溫度的影響，然后研究由坡體溫度改變引起的黏滯性變化，進而研究氣溫變化對坡內(nèi)水體滲透性的影響。

2.1 氣溫對坡體溫度的影響

在考慮水、熱耦合作用的水氣二相流計算中，基本變量主要有孔隙氣壓力、飽和度以及溫度，在每一個時間步計算結(jié)束后均能得到對應(yīng)時刻坡體內(nèi)各基本變量的值，計算域內(nèi)某一時刻坡體溫度場分布情況見圖4。由圖4可知，坡體溫度變化劇烈的區(qū)域主要集中在坡體表層，坡體內(nèi)大部分區(qū)域的溫度維持在20 ℃。

坡體溫度變化與氣溫季節(jié)性變化特性密切相關(guān)。通過坡體溫度全年的變化情況（圖5）可知，坡體溫度變化與氣溫變化具有相同的趨勢，但在時間上具有一定的滯后性。例如全年最低氣溫出現(xiàn)在12月18日，而坡體不同深度處的最低溫度均要晚于該時間。通過比較可以發(fā)現(xiàn)冬季的滯后性較夏季明顯，這主要與冬季降雨量較少，坡體的溫度變化主要由土體的熱傳導(dǎo)作用引起有關(guān)。由于夏季降雨充沛，雨水入滲引起的熱傳輸作用占主導(dǎo)地位，因此坡體溫度對氣溫變化較敏感，滯后性較冬季小。

由熱傳遞的規(guī)律可知，熱量總是由溫度高的地方向溫度低的地方傳遞。當(dāng)外界氣溫高于坡體溫度時，熱量向坡體內(nèi)傳輸，夏季坡體溫度均低于氣溫，而冬季則呈現(xiàn)出坡體內(nèi)溫度高于外界氣溫的現(xiàn)象。距地表深度越深，氣溫變化引起的溫度改變量越小。坡體內(nèi)溫度變化量通常小于氣溫改變量，在地下0.5 m處的變化幅度可達20 ℃，達到一定深度后，坡體溫度不再隨氣溫變化而發(fā)生改變，且波動小于1 ℃?？梢姎鉁刈兓瘜ζ麦w的影響范圍通常限于坡體淺層區(qū)域（約5～10 m左右），具體的深度與土體的孔隙率、含水量、降雨入滲邊界等因素有關(guān)。

2.2 坡體溫度對坡體滲流的影響

坡體滲流主要受坡內(nèi)孔隙特性以及坡內(nèi)水體流動性的影響，而往往溫度變化對均質(zhì)坡體內(nèi)的孔隙影響較小，因此溫度變化主要通過改變坡內(nèi)水體流動性從而影響坡體滲流。水的黏滯性與溫度的關(guān)系可根據(jù)公式（7）得出[21]，通過計算可知水的黏滯性受其溫度影響較為顯著，如水溫5 ℃與35 ℃相比，黏滯性相差接近一倍（圖6）。

式中：k為土體本征滲透性；μl為水的粘滯系數(shù)；klr為水的相對滲透系數(shù)。

為了便于比較，以20 ℃時的滲透系數(shù)為基準，計算各點滲透系數(shù)的相對變化量ω，圖7為坡體中段某一豎向剖面上控制點的滲透系數(shù)相對變化量ω的全年變化情況。坡體滲透系數(shù)與坡體溫度呈現(xiàn)相同的變化趨勢，具有明顯的季節(jié)性變化規(guī)律，且由于時間上滯后性的存在，坡體內(nèi)不同深度處的滲透系數(shù)并不在同一時刻取得極值。根據(jù)式（10）可知，在坡體飽和度一定的情況下，坡體溫度越高，水的黏滯性越小，坡體的滲透系數(shù)越大。距地表深度大于5 m以后，滲透系數(shù)不再隨時間發(fā)生較大變化，這跟雨水入滲的最大深度相一致[25-26]，最大入滲深度以上既有土體對溫度的傳導(dǎo)作用又有水、熱耦合的傳輸作用，而最大入滲深度以下僅有土體的熱傳導(dǎo)作用。

選取夏季和冬季具有代表性的幾天，對坡體內(nèi)不同深度處滲透系數(shù)的分布情況進行分析。由圖8可知滲透系數(shù)在沿深度向上的變化并非是線性的，隨深度的增加，滲透系數(shù)逐漸增加（或減小）并最終趨于穩(wěn)定。當(dāng)氣溫高于坡體內(nèi)溫度時（5月-6月），隨距地表深度的增加滲透系數(shù)不斷減?。环粗?，當(dāng)氣溫低于坡體內(nèi)溫度時（10月-11月），滲透系數(shù)隨距地表深度的增加而增加。

坡體不同深度處的滲透系數(shù)受氣溫影響的程度不同，將滲透系數(shù)相對變化量ω在全年的波動幅度定義為影響率。圖9為氣溫變化對滲透系數(shù)的影響率隨深度分布圖，在地下1.25 m深度處，由氣溫變化引起的滲透系數(shù)變化量可達43%。隨著深度的增加，氣溫變化對滲透系數(shù)的影響率越來越小，當(dāng)深度為8.75 m時不同季節(jié)滲透系數(shù)的變化幅度不足3%，可認為氣溫變化對木魚包滑坡體的影響深度在這一范圍內(nèi)。

3 主要結(jié)論

本文基于水氣二相流基本原理，考慮熱量在土體中的傳導(dǎo)以及水、熱耦合的傳輸作用，準確計算坡體內(nèi)溫度場分布，定量分析氣溫變化對坡體滲透系數(shù)的影響，從而研究氣溫變化與坡體滲流的相關(guān)關(guān)系。結(jié)果表明，氣溫變化對木魚包滑坡體滲流的影響深度在地表以下10 m左右，坡體滲透系數(shù)隨深度呈非線性分布，其變化率隨深度的增加而減小，并最終趨于穩(wěn)定。坡體滲透系數(shù)隨氣溫呈現(xiàn)季節(jié)性變化規(guī)律，但在時間上具有一定的滯后性，深度越深滯后性越明顯。夏季坡體滲透系數(shù)通常大于冬季，故在相同降雨條件下，夏季降雨入滲的深度更深、強度更大，由此引起土體抗剪強度減小，不利于邊坡穩(wěn)定，因此，夏季降雨更易導(dǎo)致滑坡。

參考文獻（References）：

[1] 胡冉，陳益峰，周創(chuàng)兵.考慮變形效應(yīng)的非飽和土相對滲透系數(shù)模型[J].巖石力學(xué)與工程學(xué)報，2013，32（6）：1279-1287.（HU Ran，CHEN Yi-feng，ZHOU Chuang-bing.A relative hydraulic conductivity model for unsaturated deformable soils[J].Chinese Journal of Rock Mechanics and Engineering，2013，32 （6）：1279-1287.（in Chinese））

[2] 孫冬梅，朱岳明，張明進.降雨入滲過程的水-氣二相流模型研究[J].水利學(xué)報，2007（2）：150-156.（SUN Dong-mei，ZHU Yue-ming，ZHANG Ming-jin.Water-air two-phase flow model for numerical analysis of rainfall infiltration[J].Hydraulic Engineering，2007 （2）：150-156.（in Chinese））

[3] 劉東燕，辜文杰，侯龍.降雨及地下水位對三峽庫區(qū)非飽和土邊坡穩(wěn)定性的影響[J].水利水電技術(shù)，2013，44（7）：111-115.（LIU Dong-yan，GU Wen-jie，HOU Long.Rainfall and underground water level of Three Gorges Reservoir Area of unsaturated soil slope stability[J].Water Resources and Hydropower Engineering，2013，44 （7）：111-115.（in Chinese））

[4] Tong F，Niemi A，Yang Z，et al.A numerical model of tracer transport in a non-isothermal two-phase flow system for CO2 geological storage characterization[J].Transport in Porous Media，2013，98（1）：173-192.

[5] 袁俊平，褚飛飛，季李通.飽和—非飽和非穩(wěn)定滲流數(shù)值分析中初始狀態(tài)的研究[J].水利水電技術(shù)，2006（09）：5-7.（YUAN Jun-ping，ZHU Fei-fei，JI Li-tong.Study on initial condition in numerical analysis of saturated-unsaturated seepage flow[J].Technology of Water Conservancy and Hydroelectric Power，2006 （09）：5-7.（in Chinese））

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

[7] 肖蕾，駱世威.三峽庫區(qū)某庫岸邊坡滲流及穩(wěn)定性分析[J].水利水電技術(shù)，2014（2）：108-111.（XIAO Lei，LUO Shi-wei.Analysis on seepage and stability of a bank slope in Three Gorges Reservoir Area[J].Water Resources and Hydropower Engineering，2014 （2）：108-111.（in Chinese））

[8] Ali A，Huang J，Lyamin A V，et al.Boundary effects of rainfall-induced landslides[J].Computers and Geotechnics，2014，61：341-354.

[8] 高紅貝，邵明安.溫度對土壤水分運動基本參數(shù)的影響[J].水科學(xué)進展，2011，22（4）：484-494.（GAO Hong-bei，SHAO Ming-an.Effect of temperature on soil moisture parameters[J].Advances in Water Science，2011，22 （4）：484-494.（in Chinese））

[9] 陳益峰，周創(chuàng)兵，童富果，等.多相流傳輸THM全耦合數(shù)值模型及程序驗證[J].巖石力學(xué)與工程學(xué)報，2009，28（4）：649-665.（CHEN Yi-feng，ZHOU Chuang-bing，TONG Fu-guo，et.A numerical model for fully coupled THM processes with multiphase flow and code validation[J].Chinese Journal of Rock Mechanics and Engineering，2009，28 （4）：649-665.（in Chinese））

[10] 陳玉超.凍融環(huán)境下巖土邊坡穩(wěn)定性研究初探[D].西安科技大學(xué)，2006.（CHEN Yu-chao.Preliminary study on rock and soil slope stability under the freezing-thawing condition[D].Xi′an University of Science and Technology，2006.（in Chinese））

[11] 張彬，劉寧，王紅云.基于溫度場對邊坡降雨入滲的模擬分析[J].水資源與水工程學(xué)報，2011，22（6）：141-144.（ZHANG Bin，LIU Ning，WANG Hong-yun.Simulation analysis of landslide rainfall infiltration based on temperature field[J].Journal of Water Resources and Water Engineering，2011，22 （6）：141-144.（in Chinese））

[12] 張有天，劉中.降雨過程裂隙網(wǎng)絡(luò)飽和-非飽和、非恒定滲流分析[J].巖石力學(xué)與工程學(xué)報，1997（2）：9-16.（ZHANG You-tian，LIU Zhong.Saturated/unsaturated，unsteady seepage analysis of rock fractured networks due to the percolation of rainfall[J].Chinese Journal of Rock Mechanics and Engineering，1997 （2）：9-16.（in Chinese））

[13] 吳志偉，宋漢周.基于Lu模型的淺部地溫場與滲流場耦合研究[J].水利學(xué)報，2014.（WU Zhi-wei，SONG Han-zhou.Study on shallow geothermal field and seepage field coupling based on Lu Model[J].Journal of Hydraulic Engineering，2014.（in Chinese））

[14] 劉建強，黃金廷.飽和—非飽和帶溫度動態(tài)與氣象因素關(guān)系[J].地下水，2013（02）：21-23.（LIU Jian-qiang，HUANG Jin-ting.Relation between temperature dynamic and meteorological factors in saturated-unsaturated zone[J].Groundwater，2013 （02）：21-23.（in Chinese））

[15] Tong F G，Jing L，Zimmerman R W.A fully coupled thermo-hydro-mechanical model for simulating multiphase flow，deformation and heat transfer in buffer material and rock masses[J].International Journal of Rock Mechanics and Mining Sciences，2010，47：205–217.

[16] 林鴻州，于玉貞，李廣信，等.土水特征曲線在滑坡預(yù)測中的應(yīng)用性探討[J].巖石力學(xué)與工程學(xué)報，2009（12）：2569-2576.（LIN Hong-zhou，YU Yu-zhen，LI Guang-xin，et al.On application of soil-water characteristic curves to landslide forecast[J].Chinese Journal of Rock Mechanics and Engineering，2009 （12）：2569-2576.（in Chinese））

[17] 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：892-898.

[18] Mualem Y.A new model for predicting the hydraulic conductivity of unsaturated porous media[J].Water Resources Research，1976，12：513-522.

[19] 蔡國慶，趙成剛，劉艷.一種預(yù)測不同溫度下非飽和土相對滲透系數(shù)的間接方法[J].巖土力學(xué)，2011，32（5）：1405-1410.（CAI Guo-qing，ZHAO Cheng-gang，LIU Yan.An indirect method for predicting permeability coefficients of unsaturated soils at different temperatures[J].Rock and Soil Mechanics，2011，32 （5）：1405-1410.（in Chinese））

[20] Corey A T.The interrelation between gas and oil relative permeabilities[J].Producers Monthly，1954，19：38-41.

[21] Guvanasen V，Chan T.A three-dimensional numerical model for thermos hydro mechanical deformation with hysteresis in a fractured rock mass[J].International Journal of Rock Mechanics and Mining Sciences，2000，37：89-106.

[22] Pritchett J W，Rice M H，Riney T D.Equation of state for water-carbon dioxide mixtures：implications for Baca reservoir[J].Report DOE/ET/27163-8，UC-66a，La Jolla，CA （1981）

[23] Tong F G，Jing L，Zimmerman R W.A fully-coupled finite element code for modeling thermo-hydro-mechanical processes in porous geological media.In：Proceedings 43rd US Rock Mechanics Symposium，28-30 June 2009 Asheville，NC，USA，Smeallie P，ed.Document ID：09-120.

[24] 孔祥言編著.高等滲流力學(xué)[M].合肥：中國科學(xué)技術(shù)大學(xué)出版社，1999.（KONG Xiang-yan.Seepage mechanics[M].Hefei：University of Science & Technology China Press，1999.（in Chinese））

[25] 陳偉，莫海鴻，陳樂求.非飽和土邊坡降雨入滲過程及最大入滲深度研究[J].礦冶工程，2009，29（6）：13-16，21.（CHEN Wei，MO Hai-hong，CHEN Yue-qiu.Study on rainfall infiltration process and the biggest infiltration depth for unsaturated soil slope[J].Mining and Metallurgical Engineering，2009，29 （6）：13-16，21.（in Chinese））

[26] 張士宇，王瑞鋼.降雨對高填土路堤的入滲深度的確定及有限元穩(wěn)定分析[J].路基工程，2004（5）：17-21.（ZHANG Shi-yu，WANG Rui-gang.Calculation of depth for rainfall infiltration into subgrades and stability analysis with finite methods[J].Subgrade Engineering，2004 （5）：17-21.（in Chinese））