川藏鐵路季節性粗顆粒凍土邊坡水熱過程模擬

梁樹，謝強，趙文，葉唐進, ，蔣良文

川藏鐵路季節性粗顆粒凍土邊坡水熱過程模擬

梁樹1，謝強1，趙文1，葉唐進1,2，蔣良文3

(1. 西南交通大學 地球科學與環境工程學院，四川 成都 610031；2. 西藏大學 工學院，西藏 拉薩 850000；3. 中鐵二院工程集團有限責任公司，四川 成都 610032)

以川藏鐵路東段的季節性粗顆粒凍土邊坡為研究對象，通過建立一個考慮積雪、滲流和年循環氣溫作用下的帶相變的瞬態水熱耦合的飽和-非飽和的多孔多相介質數值計算模型，結合野外實測，分析季節性粗顆粒凍土邊坡的溫度場、水分場的分布特征、凍結深度及其影響因素。研究結果表明：積雪消融入滲改變季節性粗顆粒凍土的水分場，地表可形成最大0.8 m的暫態飽和區。水分場的動態變化提高了熱傳遞速度，增強了凍結能力，邊坡凍結深度增大60%，凍結速率增大30%，融化速率增大200%。地下水熱對流作用抑制土體凍結，加速土體融化，其中坡腳地下水出露邊坡的凍融深度為地下水深埋邊坡的63%，凍結速率為79%，融化速率增大1倍。川藏鐵路新都橋地區季節性粗顆粒凍土邊坡的凍結深度為1.0 m，最大可達到1.9 m。邊坡不同位置的凍結深度不同，在坡肩處最高，坡腳處最低；進入融化期后，因凍結深度小及地下水熱流作用，坡腳處最先融化。

川藏鐵路；季節性凍土；季節性粗顆粒凍土；凍融深度；水熱耦合；數值模擬

粗顆粒土是典型的多孔介質，其孔隙大，孔隙連通率高，滲透系數大，有利于流體的流動。流體流動過程中產生的熱對流，成為粗顆粒凍土熱傳遞的重要方式。季節性粗顆粒凍土主要研究凍融循環作用下的表層非飽和?飽和粗顆粒土體，其熱傳遞過程受到積雪、積雪融化、地下水滲流、外界溫度等多方面的影響，是一個多孔多相介質帶相變的滲流場、溫度場復雜問題[1?3]。對于滲流場，非飽和粗顆粒土凍結后，孔隙中的水相變，滲流主要受到2個方面的影響。一是液態水變為固態冰以后，冰填充孔隙，降低土體的孔隙率，進而降低土體的滲透性。細顆粒土中的孔隙小，水變成冰體積增大，有可能充滿孔隙或隔斷過水通道，因此細顆粒土的滲透性受到的影響很大，往往可將其視為不透水層進行簡化處理[4]。而粗顆粒土的孔隙大，水變成冰的過程中難以充滿孔隙和隔斷過水通道。因此，季節性粗顆粒凍土的滲透性相對復雜，受到含水率、孔隙率、溫度、水的流速等多方面的綜合影響，與熱傳遞過程協調變化。二是負溫對水的基質吸力產生較大的影響，而基質吸力是非飽和滲流的重要驅動力，因此，非飽和滲流受負溫的影響[5]。對于溫度場，由于粗顆粒凍土滲流的復雜性，以滲流為基礎的熱對流作為熱傳遞的主要方式，熱場(溫度場)受到的影響也很大。此外，溫度變化，土體的導熱系數、比熱容均隨之變化，加上相變潛熱，使得溫度場受到的影響因素也較多。川藏鐵路康定?林芝段途徑川西新都橋、理塘等高寒地區，這些地區的坡麓、溝槽等平緩地帶是線路通過的首選。初步調查表明，線路首選的這些平緩地帶多為坡洪積層，以粗顆粒土為主，淺層土體季節性凍融，為粗顆粒凍土。雖然凍土地區一般遵循“寧填不挖”的基本原則，但在實際工程中，如車站、隧道進出庫等，仍然不可避免開挖，產生季節性粗顆粒凍土邊坡。這些地區的這類邊坡地下水埋深淺、流速快，地表積雪易消融入滲，不同于其他細顆粒土或多年凍土問題，其針對性的水熱耦合研究資料還不夠充分。本文對川藏鐵路新都橋地區某G318國道附近的季節性粗顆粒凍土邊坡，進行為期1 a的野外實測。結合多孔多相介質的帶相變的滲流場理論、溫度場理論及其耦合理論，建立了一個考慮積雪、滲流和年循環氣溫作用下的含相變的瞬態水熱耦合的飽和?非飽和的季節性粗顆粒凍土邊坡數值模型，分析邊坡的溫度場、水分場的分布特征、凍結深度及其影響因素。

1 滲流及熱傳導理論

1.1 非飽和滲流理論

根據滲流連續性方程，有：

式中：k(k)為()方向的滲透系數；為源匯項；為體積含水率；為時間；為總水頭；=μ/γ+, μ為基質吸力；γ為水的重度；為位置水頭。

當考慮水的相變時，基質吸力隨負溫的變化而變化。Black等[6]建立了基質吸力與負溫的線性 關系：

式中：Δμ為基質吸力的變化量；為土的溫度(低于0 ℃)；系數1 110 kPa/℃包含了凍融潛熱、對流換熱及熱傳導。

1.2 含相變的熱傳導方程

基本熱傳導方程為：

式中：為溫度；k(k)為()方向的導熱系數；q為熱源匯項；為儲熱系數，應考慮相變，引入未凍結含水率w，并有關系：

式中：為土的體積比熱容，凍結前后不同；為水的潛熱；w為未凍結含水率。

1.3 含相變的熱對流理論

當考慮滲流過程中的對流傳熱時：

式中：c為水的體積比熱容；q(q)為()方向的水流量。

1.4 計算參數

1.4.1 土水特征曲線

粗顆粒土的土水特征曲線可采用grain size- Modified Kovacs模型[7]。

式中：為孔隙率；∞為毛細水上升高度，與吸力有關；，和μ為與材料有關的系數；C為矯正系數。

1.4.2 水力傳導曲線

水力傳導曲線表示土的滲透系數與吸力之間的關系，采用適用于各類土的經驗公式的Fredlund模型[8]。

式中：k為非飽和滲透系數；k為飽和滲透系數；為最低負吸力；為最高負吸力；為到的間隔；′(e)為土水特征曲線的斜率。

1.4.3 未凍結含水率

未凍結含水率是指未凍結含水率與空隙中全部水的百分比，用于表達溫度在略低于0℃時凍融過程中的熱量變化情況。Anderson等[9]提出了一個經驗計算公式。

式中：1和為與土的性質相關的常數。

1.4.4 導熱系數

飽和狀態下的導熱系數[10]：

式中：sat為飽和導熱系數；K為土顆粒導熱系數；K為冰的導熱系數；K為水的導熱系數；為孔隙率；w為未凍結含水率。

非飽和狀態下的導熱系數為：

dry為干燥土體的導熱系數，可通過經驗公式確定：

K為系數，取值為：

為飽和度。

1.5 邊界條件

1.5.1 積雪消融入滲邊界條件

積雪出水量與環境溫度、雪層厚度、輻射及雪的狀態有關。根據陳凱[11]對季節性積雪在融雪期的雪蓋出流機制的研究結果，積雪出水量的計算方法為：

式中：h為雪層中滯后的剩余水量的深度；W為存儲的剩余水量；C3，C4為經驗系數；ρ為液態水密度；d為雪層厚度；ρ為雪層中冰的密度。

受環境溫度、積雪厚度等影響，不同階段的積雪消融入滲量不同。本文將積雪消融入滲分為3個階段：一是積雪初期，降雪量小，溫度較高，雪不易積墊，易融化，雪融化后就地入滲，積雪消融入滲量等于降雪量。二是積雪期，本階段環境溫度較低，融雪速度相對較慢，雪融水以入滲土體為主，入滲量影響因素多，不易計算。三是融雪期，融雪期的溫度回升，雪融化速度快，出水量大，雪融水往往形成地表徑流，積雪入滲量等于最大土體最大入滲量。為了簡化計算，按照工程保守原則，將積雪期的雪融水入滲量按照土體最大入滲量取值。

1.5.2 年循環地面溫度邊界條件

根據趙永峰等[12]在新都橋地區的實測及資料收集，確定年循環地面溫度邊界條件為：

2 數值模型的建立和驗證

2.1 模型建立

本文采用TEMP/W與SEEP/W有限元軟件進行季節性粗顆粒凍土邊坡的水熱數值計算。以川藏鐵路新都橋段線路附近的G318國道2 905 km處的季節性粗顆粒凍土邊坡為原型建立數值模型。原型邊坡的巖土介質為礫土，邊坡坡高4 m，坡度42°，地下水在坡腳溢出。按照原型破邊建立的數值模型1:1的數值模型，模型尺寸長×寬=32 m×12 m，網格尺寸0.2 m。

首先建立年均地面溫度作用下的穩態溫度場模型和考略地下水的穩態非飽和滲流模型。在此基礎上，建立一個年循環地面溫度作用及地下水作用下的瞬態溫度?非飽和滲流耦合模型，計算1 a后的某一時刻模型的狀態作為初始條件。最后考慮積雪消融入滲，建立一個考慮積雪和滲流含相變的季節性粗顆粒凍土邊坡的水熱耦合模型。模型計算流程及邊界條件見表1。通過實驗、工程經驗及反復試算，模型計算采用的基本參數見表2。

表1 模型建立過程

表2 模型基本參數

2.2 模型正確性驗證

野外監測了邊坡坡肩及坡后20 m處(半無限邊界)的不同深度的溫度。監測元件為PT100，監測2016年至2017年整個凍融循環周期。將2016年至2017年的凍融階段的溫度監測結果與數值模擬的結果進行對比分析，見圖1~2。從圖中可以看出，計算結果與實測結果的變化趨勢基本一致；地下?4 m處數值計算與實測的溫度均總體保持在7.8 ℃附近，表明已進入恒溫帶，季節性溫度變化帶小于4 m。地下?1.0 m處與?0.2 m處數值計算的溫度先降低后升高，實測溫度在計算的溫度曲線附近波動，差值普遍小于1。計算結果與實測結果之間的誤差一方面是由于簡化計算模型的精度所致，另一方面是本文研究的是一個凍融循環作用下的邊坡溫度場問題，將地面溫度邊界簡化為簡單的正弦函數曲線，這與實際存在一定的差別。總體而言，數值計算結果能較好地反映邊坡實際溫度場分布及變化規律，證明了所建模型及參數取值的合理性和正確性。

圖1 坡后20 m處不同深度的溫度曲線

圖2 坡肩處不同深度的溫度曲線

3 川藏鐵路季節性粗顆粒凍土邊坡的水熱數值計算

3.1 含水率分布特征

圖3為不同深度的含水率分布特征。初始狀態(2016?11?01)下，埋深0~2.5 m段處于非飽和狀態，含水率隨著深度的增加而增加；埋深超過3 m，土體位于地下水位以下，處于飽和狀態；埋深2.5~3 m，因毛細力的作用，土體仍處于飽和狀態。凍結期前期(2016?11?01~2016?12?09)，在降雪、融雪入滲的作用下，土體含水率逐漸升高。本階段的降雪量往往較小，且白晝溫度大多位于0 ℃以上，融雪速度快，雪不易積墊，往往就地消融入滲，難以形成地表徑流。故本階段的積雪消融入滲量較小，土體含水率提升有限。在凍結期中期(2016?12?09~ 2017?02?23)，積雪量大(期間最大積雪厚度可達30 cm)，加上太陽輻射和白晝氣溫的影響，融雪速度往往較大，積雪消融入滲量大，導致地表形成暫態飽和區，厚度約為0.8 m。在凍結期后期(2017? 02?23~2017?03?14)，積雪融化速度加快，更高海拔地區的積雪也逐漸加入融化，地表徑流流量大，地表仍然保持著暫態飽和區，厚度變化不大。

圖3 不同深度的含水率變化曲線

3.2 溫度場分布特征

圖4~6為溫度場云圖。初始狀態下，溫度隨深度的增加而升高。凍結前期~凍結中期，淺層地表的溫度降低，并逐漸進入凍結狀態，從圖中可以看出，坡肩處的凍結深度最大，坡腳處最小。坡腳處的凍結深度最小，是因為該處地下水埋深最小，溫度較高的地下水熱對流作用強。坡肩處的凍結深度大是因為受到溫度較低的原始地面和坡面雙重作用。進入凍結后期，土體逐漸融化，從圖可以看出，最先融化的部位為坡腳，這也是因為該處地下水埋深最小，溫度較高的地下水熱對流作用強。

圖4 2016?11?10(凍結前期)溫度場

圖5 2017?01?10(凍結中期)溫度場

圖6 2017?03?10(凍結后期)溫度場

圖7為不同深度的溫度分布特征。初始狀態下，土體溫度隨深度的增大而升高；但當埋深大于3 m時，由于土體進入地下水位以下，土體的溫度變化較小。在凍結前期，表層土體溫度迅速降低，至2016?12?09已形成凍結層，厚度約為0.2 m。出現凍結層以后，在非凍結?凍結分界線處表現出溫度隨深度變化曲線的拐點。凍結層內的土體溫度隨著深度的增大呈近線性增大，非凍結層的溫度隨深度增大呈對數形式升高。2017?02?23，土體表層溫度開始回升，至2017?03?14以后，地面溫度達到0 ℃以上，土體開始解凍。至2017?04?02，土體解凍基本結束，地面溫度上升至平均氣溫以上，至2017?04?21，溫度回升已導致地下3 m以上所有土體上升到平均氣溫以上，土體溫度隨著溫度的升高而降低，本次凍融過程結束。

圖7 不同深度的溫度變化曲線

3.3 凍融深度特征

凍融深度變化特征見圖8，圖中縱坐標為凍結?未凍結分界線的埋深，橫坐標為時間。某一時刻有一個埋深時，這個埋深為凍結深度；某一時刻有2個埋深時，較大的埋深代表凍結深度，較小埋深代表進入融化期后邊坡的融化深度。

從圖中可以看出，2016年11月底，土體開始凍結，隨后凍結深度逐漸增大，到2017?03?11時達到最大，深度約為1.2 m。此時地表已于2017? 03?06開始解凍，即最大凍結深度的時間要稍晚于開始解凍時間約5 d。土體進入解凍期后，地下水及地面溫度均高于凍結土體溫度，凍結層的底部和頂部同時解凍，至2017?03?31土體全部解凍。土體解凍從2017?03?06~2017?03?31，共25 d；土體凍結從2016?11?26~2017?03?11(2017?03?06開始解凍，但底部的凍結線仍在向下移動)，共約100 d，解凍速度約為凍結速度的4倍。

圖8 有無積雪入滲的凍結-未凍結分界線深度曲線

3.4 積雪消融入滲對土體凍融的影響

積雪消融入滲對土體凍融的影響(見圖8)主要考慮：1)在凍結初期，雪融水的溫度低于土體溫度，積雪消融入滲帶走土體熱量，有助于土體凍結。 2) 在凍結后期，土體開始融化，凍結的土體溫度低于雪融水的溫度，積雪消融入滲給土體帶來熱量，有助于土體融化。對比有積雪消融入滲和無積雪消融入滲的季節性粗顆粒凍土的凍融深度，可以看出：1) 有積雪消融入滲的最大凍結深度為1.2 m，無積雪消融入滲為0.75 m，有積雪消融入滲的凍結深度是無積雪消融的1.6倍。2) 有無積雪消融入滲的邊坡土體最大凍結深度的時間不同，無積雪消融入滲時，在外界溫度最低點時，土體凍融深度達到最大，凍結層底部邊界線保持穩定，直至地面溫度回升到0 ℃上以后約25 d才快速上升；有積雪消融入滲時，在外界溫度回升到0 ℃之后5 d，土體凍結深度達到最大，隨后凍結層底部邊界線快速上升。這反映了外界溫度達到最低溫度以后，雖然地面溫度仍低于0 ℃，但土體凍結層的底部變化難以向下移動；但如果此時有積雪消融入滲的作用，熱對流將大大提高凍結能力。3) 凍結深度隨時間的曲線斜率反映了土體的凍結速率和融化速率。有無積雪消融入滲的凍結速率分別為1.2/100=0.012 m/d和0.75/80=0.009 4 m/d，有積雪消融入滲的凍結速率是無積雪消融入滲的1.3倍。有無積雪消融入滲的融化速率分別為1.2/20=0.06 m/d和0.75/35=0.02 m/d，有積雪消融入滲的凍結速率是無積雪消融入滲的3倍。

3.5 地下水水位對凍結深度的影響

地下水作為熱源熱對流邊界，往往抑制土體的凍結，加速土體的融化。對比地下水坡腳出露和地下水埋深大(坡腳處埋深8 m)的季節性粗顆粒凍土邊坡的凍融深度，見圖9，可以看出：1)坡腳地下水出露邊坡的最大凍結深度為1.2 m，地下水深埋時為1.9 m，前者是后者的63%。2) 凍結深度隨時間的曲線斜率反映了土體的凍結速率和融化速率。地下水坡腳出露和地下水深埋時邊坡的凍結速率分別為1.2/100=0.012 m/d和1.9/125=0.015 2 m/d，前者是后者的79%。地下水坡腳出露和地下水深埋時邊坡的融化速率分別為1.2/20=0.06 m/d和1.9/65 =0.03 m/d，前者比后者增大1倍。

圖9 不同地下水深淺的凍結-未凍結分界線深度曲線

4 結論

1) 積雪消融入滲提高了季節性粗顆粒凍土邊坡熱交換的速度，增強了凍結能力。本文計算結果為季節性粗顆粒凍土邊坡的凍結深度增大60%，凍結速率增大30%，融化速率增大200%。

2) 地下水作為熱對流邊界，會抑制土體凍結，加速土體融化。本文計算結果為坡腳地下水出露邊坡的凍融深度為地下水深埋邊坡的63%，凍結速率為79%，融化速率增大1倍。

3) 川藏鐵路新都橋地區季節性粗顆粒凍土邊坡的凍結深度為1.2 m，最大可達到1.9 m。邊坡不同位置的凍結深度不同，在坡肩處最高，坡腳處最低；進入融化期后，因凍結深度小及地下水熱流作用，坡腳處最先融化。

4) 季節性粗顆粒凍土邊坡白晝時積雪出水量大，產生積雪消融入滲，在地表形成暫態飽和區，厚度約為0.8 m。

[1] 張玉芝, 杜彥良, 孫寶臣, 等. 季節性凍土地區高速鐵路路基凍融變形規律研究[J]. 巖石力學與工程學報, 2014, 33(12): 2546?2553. ZHANG Yuzhi, DU Yanliang, SUN Baochen, et al. Roadbed deformation of high-speed railway due to freezing-thawing process in seasonally frozen regions[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(12): 2546?2553.

[2] 張玉芝, 杜彥良, 孫寶臣. 季節性凍土地區高速鐵路路基地溫分布規律研究[J]. 巖石力學與工程學報, 2014, 33(6): 1286?1296. ZHANG Yuzhi, DU Yanliang, SUN Baochen. Temperature distribution in roadbed of high-speed railway in seasonally frozen regions[J]. Chinese Journal of Rock Mechanics and Engineering, 2014(12): 2546? 2553.

[3] 張媛, 董建華, 董旭光, 等. 季節性凍土區土釘邊坡支護結構凍融反應分析[J]. 巖土力學, 2017, 38(2): 574? 582. ZHANG Yuan, DONG Jianhua, DONG Xuguang, et al. Analysis of freezing and thawing of slope improved by soil nailing structure in seasonal frozen soil region[J]. Rock and Soil Mechanics, 2017, 38(2): 574?582.

[4] 習阿幸, 劉志輝, 盧文君. 干旱區季節性凍土凍融狀況及對融雪徑流的影響[J]. 水土保持研究, 2016, 23(2): 333?339. XI Axing, LIU Zhihui, LU Wenjun. Processes of seasonal frozen soil frezing-thawing and impact on snowmelt runoff in arid area[J]. Research of Soil and Water Conservation, 2016, 23(2): 333?339.

[5] 吳青柏, 蔣觀利, 蒲毅彬. 青藏高原天然氣水合物的形成與多年凍土的關系[J]. 地質通報, 2006, 25(2): 29?33. WU Qingbai, JIANG Guanli, PU Yibin. Relationship between permafrost and gas hydrates on Qinghai-Tibet Plateau[J]. Geologial Bulletin of China, 2006, 25(2): 29?33.

[6] Black P, Tice A. Comparison of soil freezing curve and soil water curve data for windsor sandy loam[J]. Water Resources Research, 1989, 25(10): 2205?2210.

[7] Aubertin M, Mbonimpa M, Bussière B, et al. A model to predict the water retention curve from basic geotechnical properties[J]. Canadian Geotechnical Journal, 2003, 40(6): 1104?1122.

[8] Fredlund D G, XING Anqing. Equations for the soil-water characteristic curve[J]. Canadian Geotechnical Journal, 1994(31): 521?532.

[9] Johansen O. Thermal conductivity of soils[D]. Trondheim, Norway. (CRREL Draft Translation, 637, 1977), ADA044002.

[10] Anderson D M, Tice A R. The unfrozen interfacial phase in frozen soil water systems, physical aspects of soil water and salt in ecosystems[J]. In Ecological Studies Springer Verlag, Berlin, 1973(4): 107?124.

[11] 陳凱. 融雪水在雪層中的凍融過程及側滲出流研究[D]. 烏魯木齊: 新疆大學, 2012. CHEN Kai. Research on freeze-thaw process and side flow of snowmelt in the snow[D]. Urumuqi: Xinjiang University, 2012.

[12] 趙永峰, 謝強, 王子江, 等. 川藏鐵路季節性凍土區粗顆粒土邊坡溫度場特征[J]. 鐵道建筑, 2017(5): 95?99. ZHAO Yongfeng, XIE Qiang, WANG Zijiang, et al. Characteristics of temperature field of coarse-grained soil slope in seasonal frozen soil region along Sichuan-Tibet Railway[J]. Railway Engineering, 2017(5): 95?99.

Hydrothermal process simulation of seasonal coarse-grain frozen soil slope on Sichuan-Tibet railway

LIANG Shu1, XIE Qiang1, ZHAO Wen1, YE Tangjin1,2, JIANG Liangwen3

(1. Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, Chengdu 610031, China;2. Tibet University Institute of Technology, Lhasa 850000, China; 3. China Railway Eryuan Engineering Group Co., Ltd, Chengdu 610032, China)

This article based on the seasonal coarse-grained frozen soil slope in the eastern section of the Sichuan-Tibet Railway, a transient hydrothermal coupled saturated-unsaturated porous multiphase media numerical model with phase transition was established. This model considered snow cover, seepage and annual cycle temperatures. The temperature field of seasonal coarse-grained frozen soil slope, moisture distribution characteristics, freezing depth and its influencing factors were analyzed combing with field measurements. The results show that the removal and infiltration of snow may change the moisture field of seasonal coarse-grained frozen soil. The surface can form a transient saturation area of up to 0.8 m. The dynamic changes of the moisture field cause the increase of heat transfer rates, the enhancement of freezing ability. It also causes slope freeze depthincreasing by 60%, the freezing rate by 30% and the melting rate by 200%. The soil freezing speed is suppressed by the boundary of the groundwater convection. The frozen-thawing depth of the outcrop slope at the slope foot is 63% of the deep-buried slope of the groundwater, the freezing rate is 79% and the melting rate is increased by 1 time. The freezing depth of seasonal coarse-grained frozen soil slope in the Xinduqiao area of the Sichuan-Tibet railway is 1.0 m, and the maximum reachable is 1.9 m. The different depths of freezing in different positions of the slope are highest at the shoulders and lowest at the foot. After entering the melting period, due to the small depth of freezing and the effect of groundwater heat flow, the first step of melting is at the foot of the slope.

Sichuan-Tibet railway; seasonal frozen soil; the seasonal coarse-grain frozen soil; freeze-thaw depth; hydrothermal coupling; numerical simulation

U213.1

A

1672 ? 7029(2019)06? 1440 ? 08

10.19713/j.cnki.43?1423/u.2019.06.012

2018?05?22

國家自然科學基金資助項目(41662020)；中國鐵路總公司科技研究開發計劃課題(2017G008-F)

趙文(1975?)，男，重慶人，副教授，從事特種巖土工程、邊坡工程治理的研究；E?mail：civil1301@swjtu.edu.cn

(編輯 涂鵬)