武夷山地區(qū)短時(shí)凍期邊坡淺層溫濕度變化的數(shù)值模擬

袁文豪，王榮，林偉

（1.福建船政交通職業(yè)學(xué)院 土木工程學(xué)院，福建 福州 350007；2.武夷學(xué)院 土木工程與建筑學(xué)院，福建 南平 354300；3.福州大學(xué) 土木工程學(xué)院，福建 福州 350116）

福建省位于我國東南沿海，境內(nèi)短時(shí)凍土（土體凍結(jié)時(shí)間在幾小時(shí)至幾天）主要分布在閩西北，面積約2 萬km2。福建短時(shí)凍土區(qū)既有土坡廣泛在每年降雨量很少的冬季出現(xiàn)淺層失穩(wěn)現(xiàn)象，但破壞原因存在困惑?，F(xiàn)場調(diào)查表明，短時(shí)凍區(qū)所引發(fā)的土坡淺層短時(shí)反復(fù)凍融是影響該區(qū)域內(nèi)土質(zhì)邊坡冬季失穩(wěn)的重要原因。因此，研究短時(shí)凍區(qū)氣候下的土坡淺層溫濕變化具有重要意義。與此同時(shí)，國內(nèi)外學(xué)者針對(duì)多年與季節(jié)凍土工程的水熱力耦合進(jìn)行大量研究，取得許多重要的成果[1-3]。然而上述研究對(duì)短時(shí)凍區(qū)土質(zhì)邊坡水熱變化研究還較少。事實(shí)上，南方短時(shí)凍土與多年、季節(jié)凍土相比，其凍結(jié)深度、凍結(jié)及熱融過程均不盡相同[4-5]。基于此，依據(jù)武夷山地區(qū)冬季典型短時(shí)凍融氣候，采用有限差分原理編制凍土水熱耦合程序，模擬武夷山地區(qū)典型短時(shí)凍融氣候條件下邊坡淺層溫濕度場變化規(guī)律。通過研究，可以掌握武夷山地區(qū)典型氣候下溫度場和水分場的變化規(guī)律，為下一步邊坡力學(xué)場的分析提供重要依據(jù)，從而為武夷山地區(qū)典型氣候條件下邊坡工程的設(shè)計(jì)、施工提供借鑒指導(dǎo)。

1 凍土水熱耦合模型與求解

1.1 水熱耦合方程與程序編制

凍融條件下，水熱耦合運(yùn)移問題必須聯(lián)立求解土壤水分運(yùn)動(dòng)和熱流兩個(gè)基本方程[2]，其中水分運(yùn)動(dòng)方程主要基于質(zhì)量守恒定律和非飽和達(dá)西定律而來，對(duì)于水熱耦合問題，水分運(yùn)動(dòng)方程中多了一項(xiàng) “含冰率”；溫度方程基于能量守恒定律和Fourier 熱傳導(dǎo)理論而來。只考慮土壤中的水熱場，忽略應(yīng)力場和溶質(zhì)場，且短時(shí)凍土地區(qū)氣態(tài)水和溫度梯度的影響較小，基于Harlan 模型可得非飽和凍融條件下土壤水熱一維耦合模型為:

式中:θ、θI分別為含水率和含冰率；D（θ）、K（θ）分別為非飽和土水分?jǐn)U散率和導(dǎo)水系數(shù)；T、t 分別為溫度及時(shí)間；Cν為土壤體積比熱容；Kh為熱導(dǎo)率；LI為冰的融解潛熱，一般近似取335 J/g；ρω、ρI分別為水和冰的密度，兩者依次取值1 g/cm3和0.917 g/cm3。

基于交替隱式差分（ADI）方法對(duì)式（1）、式（2）進(jìn)行求解，采用MATLAB 語言編制程序。

1.2 參數(shù)取值

方程（1）中的導(dǎo)水系數(shù)K(θ)通過VG 模型求解，即通過已知的飽和導(dǎo)水率求得土體非飽和導(dǎo)水率[6]，方程如下:

式中:Ks為飽和導(dǎo)水率；Se為有效飽和度；θr和θs分別為殘余含水率和飽和含水率，n 為土水特征曲線參數(shù)。

選取武夷山地區(qū)典型土質(zhì)，按《公路工程試驗(yàn)規(guī)程》 試驗(yàn)采用變水頭試驗(yàn)測得殘積土飽和導(dǎo)水率，土水特征曲線參數(shù)通過定量濾紙法并利用Matlab 中的Isqcurvefit 函數(shù)擬合得出，具體參數(shù)如表1 所示。

表1 水分運(yùn)動(dòng)方程參數(shù)取值Table 1 The parameters of the water movement equation are taken℃

非飽和土水分?jǐn)U散率由以下經(jīng)驗(yàn)公式給出[5，7]:

方程（2）中土壤體積比熱容Cv按凍融土各成分質(zhì)量加權(quán)平均求得，主要考慮組成成分為土骨架、水和冰，各類土骨架比熱容參考值參考文獻(xiàn)[8]，根據(jù)室內(nèi)試驗(yàn)測得的殘積土土骨架類型可得土骨架比熱容，然后根據(jù)水和冰已知的比熱容可得土壤體積比熱容Cv=1 196 J/(cm3·℃)。

熱導(dǎo)率Kh參考文獻(xiàn)[9]中的公式，并根據(jù)室內(nèi)土工試驗(yàn)測得的干密度求得，其值為5 568 J/(cm·℃·s)，未凍含水率與負(fù)溫關(guān)系按下列經(jīng)驗(yàn)公式給出:

1.3 邊界與初始條件

初始條件:溫度T（t=0，z）=5℃，總含水率ω（t=0，z）=0.20，未凍含水率ωl（t=0，z）=0.20，含冰率ωI（t=0，z）=0。

選取武夷山地區(qū)2009 年1 月10 日至2009 年1月16 日共7 d 的短時(shí)凍融氣象數(shù)據(jù)（如表2 所示），計(jì)算模型假定深度為2 m，且地下水埋深大，毛細(xì)水無法為凍結(jié)區(qū)補(bǔ)給水量。坐標(biāo)系取地表為原點(diǎn)，向下為正。

表2 2009 年1 月10 日—2009 年1 月16 日武夷山地區(qū)氣溫資料Table 2 Temperature data for the Wuyishan region from January 10,2009 to January 16,2009

1.4 模型驗(yàn)證

為驗(yàn)證上述計(jì)算模型的準(zhǔn)確性，選取武夷山地區(qū)境內(nèi)省道303 線K 364+900 處邊坡作為本次試驗(yàn)的監(jiān)測點(diǎn)，該地區(qū)具有典型的南方短時(shí)凍土的氣候特點(diǎn)。監(jiān)測該地區(qū)土壤層內(nèi)固定位置溫度和體積含水率的變化規(guī)律，并與數(shù)值模擬結(jié)果對(duì)比。監(jiān)測儀器采用PC-3 型移動(dòng)式自動(dòng)氣象站，該儀器可監(jiān)測環(huán)境的溫濕度。監(jiān)測時(shí)間為2008 年12 月27 日，測點(diǎn)1 和測點(diǎn)2 深度分別為12 和27 cm。2 個(gè)測點(diǎn)的溫度和體積含水率監(jiān)測值和模擬值分別如圖1 和如圖2 所示。

圖1 監(jiān)測點(diǎn)溫度對(duì)比Fig.1 Temperature contrast of monitoring point

圖2 監(jiān)測點(diǎn)體積含水率對(duì)比Fig.2 Comparison of volumetric moisture content at monitoring points

由圖1 所示，2 個(gè)監(jiān)測點(diǎn)溫度的檢測值和模擬值變化趨勢相同，雖然存在一定誤差，但是誤差總體上均小于1 ℃，產(chǎn)生誤差的可能原因有如下:一是數(shù)值計(jì)算有限差分節(jié)點(diǎn)溫度在0 ℃至凍結(jié)溫度范圍，未考慮土壤參數(shù)的變化；二是MATLAB 有限差分計(jì)算采用方塊幾何模型計(jì)算，而實(shí)際情況則為邊坡模型，邊界條件按坡度做出折算后再代入模型，兩者存在一定偏差；三是由于操作誤差，現(xiàn)場監(jiān)測點(diǎn)的埋設(shè)尺寸存在一定偏差。由圖2 所示，2 個(gè)監(jiān)測點(diǎn)體積含水率的實(shí)測值與模擬值基本相同，其日變化范圍大致在0.001～0.002，主要原因是受到蒸發(fā)的影響。綜上，雖然數(shù)值計(jì)算結(jié)果與實(shí)測值存在一定偏差，但該偏差在可控范圍之內(nèi)，因此，可采用本數(shù)值模型計(jì)算短時(shí)凍土邊坡溫濕度。

2 短時(shí)凍土邊坡水熱耦合模擬結(jié)果與分析

2.1 邊坡溫度變化

基于2009 年1 月10 日至2009 年1 月16 日共7 d的短時(shí)凍融氣候條件，通過數(shù)值模擬得到7 d 的邊坡溫度總變化如圖3 所示。由圖可見，邊坡淺層溫度狀況隨氣溫的變化基本呈周期性的變化規(guī)律，且存在邊坡淺層溫度滯后氣溫的現(xiàn)象。同時(shí)，邊坡淺層（0～15 cm）最低溫度在每天均達(dá)到零下，最大值約-5.1 ℃；而最高溫度可達(dá)13.5 ℃左右。

圖3 2009 年1 月10 日—2009 年1 月16 日不同深度邊坡溫度變化Fig.3 Slope temperature changes at different depths from January 10,2009 to January 16,2009

以2009 年1 月11 日為例，邊坡不同深度和不同時(shí)刻處溫度計(jì)算結(jié)果如圖4 所示。

從圖4 中可以看出:

（1）表層土體從夜間21:00 到次日凌晨5:00 之間，土體溫度低于0 ℃。地表溫度剛開始下降時(shí)，地表淺層降溫迅速，5 cm 處2 h 之內(nèi)溫度大約降低3 ℃，但由于低溫天氣持續(xù)時(shí)間有限，而氣溫對(duì)土體的影響存在一定滯后性，故影響深度十分有限。

（2）未凍土段斜率明顯大于凍土段，說明未凍土段的溫度梯度小于凍土段。隨著持續(xù)凍結(jié)，未凍土段和凍結(jié)段的斜率都在提高，說明隨著時(shí)間的增加，凍土層和未凍土層溫度梯度逐漸減小，溫差變小。原因在于在凍結(jié)之初由于地層表面為負(fù)溫，而土層內(nèi)部為正溫，存在較大的溫差，所以凍土層的溫度梯度較大；然而隨著時(shí)間的增加，地層間的溫差漸漸減小，溫度梯度變小。

（3）從0 點(diǎn)至次日0 點(diǎn)，表面土體溫度約8 h 處于0 ℃以下，最低溫度達(dá)-5 ℃，5 cm 處土體約7 h 的時(shí)間處于0 ℃以下，但最低溫度只有-3 ℃，而10 cm處土體約有5 h 的時(shí)間處于0 ℃以下，最低溫度高于-2 ℃，15 cm 處土體位于0 ℃以下時(shí)間小于3 h，20 cm 深度處土體溫度高于0 ℃。鑒于土體凍結(jié)溫度通常低于0 ℃，因此，判斷凍土影響深度為10 cm 左右。

2.2 邊坡水分變化

圖5～圖10 為在7 d 的凍融過程中，土體在不同時(shí)刻不同深度總含水率、未凍含水率和含冰率的變化曲線圖。

圖5 2009 年1 月10 日—2009 年1 月16 日不同深度邊坡總含水率變化Fig.5 The change of total moisture content of slope in different depth from October 10,2009 to January 16,2009

從圖5 和圖6 可見，經(jīng)過一定的凍結(jié)時(shí)間后，上部土層的總體積含水率增大，中間層位土體總體積含水率略有減少，而下部土體總體積含水率基本不變。如經(jīng)過約3 h 的凍結(jié)，地表的含水率由初始的20%提升到25%，地層在0～5 cm 厚度范圍內(nèi)含水量提高明顯，5 至20 cm 深度范圍內(nèi)總含水量降低，而在20 cm深地層以后總含水量基本沒有變化。這是因?yàn)殡S著冰凍，地層內(nèi)毛細(xì)水持續(xù)不斷向凍結(jié)層移動(dòng)，而冰凍層內(nèi)毛細(xì)水或間隙水則移向地表，導(dǎo)致地表含水量顯著提高，而冰凍層以下土層的間隙水或毛細(xì)水移出后，沒有及時(shí)補(bǔ)給，造成含水率降低。

圖6 2009 年1 月11 日邊坡水分變化Fig.6 Moisture change of slope in January 11,2009

從圖7 和圖8 可見，隨著冰凍的持續(xù)，地表以下15～20 cm 深范圍內(nèi)的未凍地層含水率逐漸降低，其含水率由未凍前的20%降低至16%左右，土層孔隙則由冰填充。但隨著凍結(jié)土體的融化，土體未凍含水率大幅度上升，這是因?yàn)楸硗量偤试趦鼋Y(jié)過程中不斷上升，此時(shí)凍結(jié)的土全部融化，從而導(dǎo)致未凍含水率大大上升，如地表處未凍含水率由凍結(jié)時(shí)的16%大幅度上升至后來的25%，基本趨于飽和狀態(tài)。隨后由于重力勢的作用，水分向下遷移，未凍含水率恢復(fù)至正常水平。

圖7 2009 年1 月10 日—2009 年1 月16 日邊坡未凍含水率變化Fig.7 The change of unfrozen moisture content of slope from October 10,2009 to January 16,2009

圖8 2009 年1 月11 日不同深度不同時(shí)刻邊坡未凍含水率變化Fig.8 Change of unfrozen moisture content of slope at different depth and time in January 11,2009

從圖9 和圖10 中可見，土體含冰率主要出現(xiàn)在淺層土體約10 cm 深度范圍內(nèi)，即武夷山地區(qū)土體凍結(jié)深度約為10 cm，其中，表土含冰率最高可達(dá)8%，隨著深度增加，含冰率顯著減小。

圖9 2009 年1 月10 日—2009 年1 月16 日不同深度邊坡含冰率變化Fig.9 The change of ice content in slopes with different depths from October 10,2009 to January 16,2009

圖10 2009 年1 月11 日不同深度不同時(shí)刻邊坡含冰率變化Fig.10 Changes of ice content in slope at different depth and time in January 11,2009

綜上可以得知: 無論是未凍層含水量，還是總含水率，地表淺層的間隙水或毛細(xì)水含量變化都最顯著，而距地表10 cm 深度以下地層總含水率和未凍水含量基本未發(fā)生較大變化。這說明凍結(jié)層水分遷移現(xiàn)象最明顯，水分聚集最嚴(yán)重，武夷山地區(qū)凍結(jié)層厚度約為10 cm。

3 結(jié)論

運(yùn)用凍融狀態(tài)下的溫濕耦合模型，以武夷山地區(qū)的氣候條件為背景，對(duì)該地區(qū)邊坡淺層溫度和濕度變化進(jìn)行研究，得出福建省北部典型短時(shí)凍土區(qū)邊坡溫度場和濕度場變化規(guī)律，結(jié)論如下:

（1）邊坡淺層溫度狀況隨氣溫的變化基本呈周期性的變化規(guī)律，且存在邊坡淺層溫度滯后氣溫的現(xiàn)象。同時(shí)，邊坡淺層（0～15 cm）最低溫度在每天均達(dá)到零下，最大值約-5.1 ℃；而最高溫度可達(dá)13.5 ℃左右，表層土體凍結(jié)時(shí)間在夜間21:00 到次日凌晨5:00 之間。

（2）冰凍之初，未凍土層的溫度梯度小于凍土層，而隨著冰凍的持續(xù)，凍土層和未凍土層內(nèi)溫度梯度降低，溫差逐漸變小。

（3）不論是總含水率，還是未凍水含量，地表淺層的含水率變化幅度都最大，而10 cm 以下的含水率曲線平緩，總含水率和未凍水含量輕微下降或不變。說明水分遷移現(xiàn)象最明顯且水分聚集最嚴(yán)重的部位在于凍結(jié)層本身，武夷山地區(qū)凍結(jié)層厚度約為10 cm。