韓洪武，穆彥虎，虞 洪，丁澤琨，陳 領(lǐng)

(1. 青海省水利水電工程局有限責(zé)任公司，青海 西寧 810001； 2. 中國科學(xué)院西北生態(tài)環(huán)境資源研究院 凍土工程國家重點實驗室，甘肅 蘭州 730000； 3. 中國科學(xué)院大學(xué)，北京 100049)

多年凍土指地表以下一定深度內(nèi)地溫在0 ℃及0 ℃以下至少連續(xù)存在兩年的巖土層。我國為世界第三凍土大國，多年凍土分布面積約占國土面積的22%，主要分布在青藏高原、大小興安嶺及西部高山地區(qū)[1]。作為一種客觀特殊地質(zhì)體，多年凍土是地質(zhì)歷史和氣候變遷背景下受區(qū)域地理環(huán)境、地質(zhì)構(gòu)造、巖性、水文和植被特征等因素共同影響下，通過地氣間物質(zhì)與能量交換發(fā)育而成的。從物質(zhì)組成來看，凍土是由礦物顆粒、冰、未凍水和氣體組成的一種四相體材料[2]。與常規(guī)巖土不同，凍土的物理、力學(xué)、水力性質(zhì)及伴隨凍結(jié)和融化過程中的熱、力、變形現(xiàn)象均與其溫度變化過程密切相關(guān)[3]。

在多年凍土區(qū)進(jìn)行工程建設(shè)，工程活動對包括局地地形、地貌、植被、淺層巖土體類型、地表水體、徑流等局地因素的改變將不可避免地打破原有的地氣能量平衡，進(jìn)而影響到下伏多年凍土水熱狀況，引發(fā)多年凍土的退化過程，具體可表現(xiàn)為活動層厚度增加、地溫升高、融化夾層發(fā)育[4-7]。對于渠道工程，涉及的開挖工程和過水等將對活動層季節(jié)凍融過程和下部多年凍土水熱狀況產(chǎn)生長期且顯著的影響，而這一過程反過來將大幅降低凍土地基承載力，引發(fā)顯著的沉降變形過程，由此威脅到水工結(jié)構(gòu)物的長期穩(wěn)定性，甚至結(jié)構(gòu)功能失效[8]。目前，受工程實踐需求所限，國內(nèi)外凍土區(qū)渠道相關(guān)研究主要集中在季節(jié)凍土區(qū)，包括渠基土的水熱監(jiān)測與模擬、襯砌結(jié)構(gòu)的凍脹機理與過程、適應(yīng)凍脹的襯砌結(jié)構(gòu)優(yōu)化與措施研發(fā)、循環(huán)凍融作用對渠基土、混凝土等材料工程性能的影響等方面[9-14]。而在多年凍土區(qū)，相關(guān)方面的研究十分有限，主要集中在俄羅斯、加拿大等高緯度多年凍土區(qū)的一些水利大壩的設(shè)計、建設(shè)及運營過程中壩體的熱狀況、凍融界限及包括凍脹、融沉等壩體變形的監(jiān)測和模擬[15-19]。

以高海拔多年凍土區(qū)某渠道工程為研究對象，開展渠道工程和多年凍土地基相互作用研究。針對長期過水及氣候變暖影響，在考慮凍融土體內(nèi)水分遷移、冰水相變及土體未凍水含量與溫度非線性關(guān)系基礎(chǔ)上，構(gòu)建凍融土體水-熱耦合數(shù)學(xué)模型，利用數(shù)值模擬手段開展渠道下伏凍土地基長期熱狀況演化過程模擬研究，并通過已有現(xiàn)場監(jiān)測數(shù)據(jù)對模型進(jìn)行驗證。

1 水熱耦合計算模型

為研究渠道下部多年凍土地基長期熱狀況演化規(guī)律，需要建立凍融土體水熱耦合數(shù)學(xué)模型。模型建立過程中的基本假設(shè)包括：土體均勻、連續(xù)且各向同性；凍結(jié)、融化土體內(nèi)的水分遷移均滿足達(dá)西定律；考慮水勢梯度作為凍融土體中液態(tài)水遷移的驅(qū)動力，忽略溫度梯度的影響。

1.1 質(zhì)量守恒方程

根據(jù)質(zhì)量守恒定律，在不考慮水汽遷移的情況下，凍土中總等效體積含水量θ可以寫成：

式中：θ為總等效體積含水量（m3/m3）；θu為液態(tài)水體積含量（m3/m3）；θi為體積冰含量（m3/m3）；qlh為土壤水勢梯度下液態(tài)水通量密度（m/s）；ρi和ρw分別為冰密度和水密度（kg/m3）；t為時間。

根據(jù)Harlan模型，凍土中的液態(tài)水遷移與非飽和融土中液態(tài)水遷移相似[20]。凍土中的液態(tài)水遷移可以用Richards方程來描述，且主要受水勢梯度和導(dǎo)水系數(shù)的影響[21]。在未凍土和凍土中，水勢梯度和導(dǎo)水系數(shù)有差異顯著，但凍結(jié)土體中水分遷移仍可認(rèn)為遵循達(dá)西定律[22]。因此，只考慮水勢梯度作為凍結(jié)土體和未凍土體中液態(tài)水遷移的驅(qū)動力，則液態(tài)水的通量密度[23]可以寫成:

式中：Klh為水勢梯度下土體液態(tài)水導(dǎo)水系數(shù)（m2/s）；y為縱向坐標(biāo)（m）；h為壓力水頭（m）?？紤]將總等效體積含水量作為變量，凍融土體中液態(tài)水質(zhì)量守恒方程可以寫為：

1.2 水力參數(shù)

凍融土體中液態(tài)水遷移與其壓力水頭相關(guān)，可通過土水特征曲線（SWCC）計算壓力水頭。本文中，凍融土體的水力參數(shù)可以用van Genuchten模型和Mualem模型[24-25]進(jìn)行描述：

式中：Se為有效飽和度；Ks為飽和土體導(dǎo)水系數(shù)（m/s）；θl、θs、θr分別為液態(tài)水含量、飽和液態(tài)水含量和殘余含水量（m3/m3）；α為土體進(jìn)氣值的倒數(shù)（m?1）；m= (1 – 1/n)、n、l為經(jīng)驗參數(shù)，文獻(xiàn) [25]建議l取 0.5。

1.3 能量守恒方程

已有試驗研究表明，與熱傳導(dǎo)和冰水相變潛熱相比，土體中對流傳熱量級非常小，在凍融土體的傳熱分析中一般可以忽略[21]。因此，僅考慮熱傳導(dǎo)和冰水相變潛熱，則凍融土體的能量守恒方程[26]可以寫為：

式中：T為溫度；Cm、λm分別為土體等效體積比熱容（J/m3/K）和等效導(dǎo)熱系數(shù)（W/m/K）；L為冰水相變潛熱（3.34×105J/kg）。根據(jù)顯熱容法，土體等效體積比熱容和等效導(dǎo)熱系數(shù)可以寫為：

式中：Tf± ΔT為相變區(qū)間；Cu和λu分別為凍結(jié)溫度以上土體比熱容和導(dǎo)熱系數(shù)；Cf和λf分別為凍結(jié)溫度以下土體比熱容和導(dǎo)熱系數(shù)；Ls為含水巖土介質(zhì)單位體積相變潛熱。

1.4 聯(lián)系方程

對于凍結(jié)區(qū)，以上質(zhì)量守恒、能量守恒兩個方程涉及溫度、體積含冰量、液態(tài)水含量3個未知量，未知量大于方程數(shù)，需要增加聯(lián)系方程。土體凍結(jié)后，未凍水的含量取決于溫度、壓力、水鹽度、礦物學(xué)、土壤比表面積和土壤表面化學(xué)成分等因素。基于已有理論和試驗研究結(jié)果[21]，采用下式來表示土體凍結(jié)過程中的最大未凍水含量：

式中：a和b為土體試驗參數(shù)。體積液態(tài)水含量和含冰量可通過等效體積含水量和溫度確定：

式中：Tf為土體凍結(jié)溫度(℃)。已有研究表明，土體的凍結(jié)溫度并不是一個定值，并且只有土體中的液態(tài)水大于最大未凍水含量時才會有冰生成。

2 數(shù)值計算模型及邊界條件

2.1 數(shù)值計算模型

以高海拔多年凍土區(qū)某渠道為研究對象，工程區(qū)平均海拔4 500 m以上，區(qū)域多年平均氣溫為?5.1 ℃，年平均降水量313.8 mm，屬于大片連續(xù)多年凍土區(qū)。根據(jù)現(xiàn)場鉆探及相關(guān)資料，區(qū)域內(nèi)活動層厚度（季節(jié)融化層深度）在2~3 m，巖性大部分為含礫砂土及礫砂，呈散體狀。多年凍土含冰量特征以富冰凍土為主，局部為飽冰-含土冰層，多分布于多年凍土上限至5 m深度范圍。

考慮到渠道結(jié)構(gòu)對稱性，以渠道中心線為對稱軸，取渠道結(jié)構(gòu)一半建立物理模型。根據(jù)現(xiàn)場工程地質(zhì)勘查結(jié)果和渠道設(shè)計方案，渠基土自上而下簡化為3層，即砂礫土、粉質(zhì)黏土和強風(fēng)化泥巖，建立數(shù)值計算模型如圖1所示。圖1中，渠道底面寬度為25 m，渠道邊坡按1∶3放坡，坡長15 m，渠內(nèi)水深依據(jù)目前實測水深設(shè)置為1 m?？紤]到渠道工程沿水流方向包括結(jié)構(gòu)尺寸及氣溫和水溫差異性不大，渠道地基長期熱狀況的模擬可以簡化為二維非穩(wěn)態(tài)溫度場計算問題。

圖1 渠道物理模型（單位：m）Fig. 1 Physical model of the canal (unit: m)

結(jié)合現(xiàn)場鉆孔巖芯實測干密度和含水量，以及室內(nèi)凍土熱物理參數(shù)研究結(jié)果[21-22, 26-27]，表1給出了3層土體的熱物理參數(shù)。

表1 土體熱物理參數(shù)Tab. 1 Physical-thermal parameters of soil layers

結(jié)合現(xiàn)場工程地質(zhì)勘查和多年凍土地溫實測數(shù)據(jù)，渠道沿線多年凍土地溫分布差異顯著，10~15 m深度年平均溫度（TMAGT）最高為?0.3 ℃，最低可達(dá)?1.5 ℃。依據(jù)《青藏鐵路多年凍土勘測暫行規(guī)定》中關(guān)于凍土含冰量的劃分，渠道沿線凍土體積含冰量（iv）情況則包括少冰（iv≤10%）、多冰（10%

表2 不同含冰類型土的熱物理參數(shù)Tab. 2 Physical-thermal parameters of frozen soils with different ice contents

表3 不同含冰類型土的視比熱容Tab. 3 Apparent specific heat of frozen soils with different ice contents 單位：J·kg?1·℃?1

2.2 邊界條件

為考慮氣候變暖影響，收集整理了距工程區(qū)最近氣象站數(shù)據(jù)。結(jié)果表明，自1961年至2018年期間工程區(qū)氣候呈現(xiàn)顯著的暖濕化趨勢，氣溫上升速率為0.33 ℃/(10 a)。根據(jù)附面層理論、現(xiàn)場觀測和工程區(qū)氣溫上升速率，將數(shù)值模型中邊界LC和CE設(shè)置為天然地表溫度邊界：

從實際傳熱過程看，渠道過水截面部分為第三類對流換熱邊界。然而，受限于邊界對流換系數(shù)選取困難，目前相關(guān)模擬研究中，常常以實測水體溫度作為邊界條件，即采用第一類邊界，在簡化問題的同時確保計算結(jié)果的可靠性[28]。為此，本文在渠道過水截面部分采用第一類邊界條件。已有研究表明，湖泊水體溫度與環(huán)境氣溫存在良好的線性關(guān)系[29]：

式中：Tw(t)、Ta(t)為同一時間尺度上的水溫和氣溫；A和B為常數(shù)。

根據(jù)林戰(zhàn)舉等[30]對工程附近區(qū)域多個湖泊的實際監(jiān)測水溫和氣象站氣溫數(shù)據(jù)，可確定式（13）中的A、B值。

在不考慮地下水輸入、人工熱輸入和防風(fēng)設(shè)置等因素影響時，同一地區(qū)A、B值差異不大[29]。因此，根據(jù)工程區(qū)氣象站氣溫監(jiān)測數(shù)據(jù)，并考慮區(qū)域氣溫升溫速率，可將邊界AB和BL設(shè)定為渠道水溫邊界：

模型中KG為對稱邊界，EF為絕熱邊界。底邊GF為熱流邊界，根據(jù)現(xiàn)場測溫所得該深度地溫梯度和土層熱物理參數(shù)確定熱流值為0.03 W/m2。

渠道施工過程中，淺層土體的開挖導(dǎo)致下伏多年凍土直接暴露在空氣中，勢必對多年凍土的熱狀況產(chǎn)生一定影響。然而，相較于后期長期運營過程，施工過程持續(xù)時間短，且所產(chǎn)生的熱擾動相對較小，因此在本文模擬預(yù)測中未考慮施工開挖過程對下伏多年凍土的熱影響。

3 模擬結(jié)果及分析

3.1 模型驗證

圖2給出了渠道開挖1年后天然地面不同深度地溫的數(shù)值模擬結(jié)果與現(xiàn)場監(jiān)測結(jié)果的對比情況。由監(jiān)測結(jié)果可以看出，天然地面凍土上限保持在2.1 m左右。數(shù)值模擬結(jié)果顯示，最大季節(jié)融化深度時的多年凍土上限同樣在2.1 m左右，與現(xiàn)場監(jiān)測結(jié)果接近。自天然地面以下15 m深度處，現(xiàn)場監(jiān)測結(jié)果和數(shù)值模擬結(jié)果均顯示多年凍土地溫為?0.3 ℃。從兩者的對比看，本文數(shù)學(xué)模型及所采用的熱物理參數(shù)和邊界條件能夠較好地反映該地區(qū)多年凍土熱狀況，包括活動層厚度和下伏多年凍土地溫均表現(xiàn)出較好的一致性。但從活動層范圍內(nèi)及多年凍土頂板附近的地溫對比看，模擬結(jié)果和實測結(jié)果仍存在一定的差異，這些差異與數(shù)值模型中土層的簡化及實測過程中測點深度的確定有關(guān)。

圖2 開挖1年后天然場地不同深度地溫實測與模擬對比Fig. 2 Field observed and numerical simulated temperature profiles at natural ground after the canal construction

3.2 不同年平均地溫條件下渠道地基溫度場

圖3給出了少冰凍土條件下，年平均地溫分別為?0.5、?1.0和?1.5 ℃時渠道開挖后第50年10月15日地基溫度場分布云圖。從圖3（a）可以看出，當(dāng)年平均地溫為?0.5 ℃時，開挖后第50年10月15日天然地表處的最大季節(jié)融化深度為3.5 m，相較于初始最大季節(jié)融化深度2.1 m，增加了約1.5 m，年均增加速率為3 cm/a。這一量值與目前青藏高原活動層增加速率監(jiān)測結(jié)果接近[31]。當(dāng)年平均地溫為?1.0、?1.5 ℃時，兩者的最大季節(jié)融化深度增速較年平均地溫為?0.5 ℃條件下小，反映了氣候變暖條件下不同年平均地溫多年凍土退化響應(yīng)的差異[32]。

圖3 不同年平均地溫條件下渠道開挖后第50年10月15日地基溫度場Fig. 3 Temperature fields of permafrost subgrades with different MAGTs on October 15, 50 years after the canal construction

從圖3可以看出，受渠道過流影響，50年后渠底處、坡腳處及天然地面靠水渠的部分下部30 m范圍內(nèi)多年凍土已經(jīng)退化，反映了渠道強烈的熱侵蝕作用，且這種熱侵蝕在垂向和水平向均比較顯著。多年凍土的退化勢必引發(fā)渠道底部和岸坡的差異融沉變形，導(dǎo)致渠道結(jié)構(gòu)和形態(tài)的變化。不同年平均地溫條件下，渠道的側(cè)向熱侵蝕存在顯著差異，當(dāng)多年凍土地溫為?0.5 ℃時，自岸坡向外約10 m范圍內(nèi)下部多年凍土已退化，而當(dāng)多年凍土地溫為?1.0和?1.5 ℃，岸坡下部仍有多年凍土分布。這一結(jié)果體現(xiàn)了不同年平均地溫多年凍土熱惰性的顯著差異。靠近渠道坡腳部分多年凍土的退化和岸坡下顯著的溫度梯度，可導(dǎo)致岸坡的滑塌和失穩(wěn)，并由此引發(fā)渠道形態(tài)的改變。

3.3 不同含冰量條件下渠道地基溫度場

圖4給出了年平均地溫為?0.5 ℃時，3種含冰量條件下渠道開挖后第50年10月15日地基溫度場分布云圖?？梢钥闯?，3種含冰量條件下渠道地基溫度場差異顯著。當(dāng)渠道下部多年凍土含冰量為少冰時，渠道下部以及自坡肩向外約10 m范圍內(nèi)多年凍土已經(jīng)退化。隨著含冰量的增加，多年凍土熱惰性顯著增加，因此其退化速率明顯減小。與不同地溫條件下的對比（圖3）可以發(fā)現(xiàn)，相較于地溫，多年凍土的含冰量對氣候變暖和人為擾動背景下多年凍土退化響應(yīng)影響更加顯著。這也反映了含冰量作為凍土工程地質(zhì)條件的關(guān)鍵指標(biāo)，決定著凍土工程地基熱力穩(wěn)定性和基礎(chǔ)與結(jié)構(gòu)物長期穩(wěn)定性。

圖4 不同含冰量下渠道開挖后第50年10月15日地基溫度場Fig. 4 Temperature fields of permafrost subgrades with different ice contents on October 15 of the 50th year after the canal construction

當(dāng)多年凍土含冰量為多冰和富冰時，渠道運營50年后其下部仍有多年凍土存在。但受氣候變暖和過水影響，其季節(jié)融化深度顯著增加，渠道下部形成了一個“鍋底狀”的融化盤。當(dāng)多年凍土為多冰凍土?xí)r，渠道中心鍋底狀融化盤深度可達(dá)渠底以下約20 m，而當(dāng)多年凍土為富冰凍土?xí)r，融化盤深度為渠底以下約15 m。含冰量不同時，岸坡外多年凍土的退化過程也存在明顯不同，包括季節(jié)融化深度和多年凍土上限以下地溫分布，含冰量越高，其季節(jié)融化深度越小，隨時間增加過程越緩慢。

3.4 渠道典型位置地溫狀況時間演化過程

從以上分析可以看出，在氣候變暖和渠道過水影響下，渠道不同位置多年凍土的退化過程和速率明顯不同。為分析渠道不同位置地溫狀況時間演化過程，以年平均地溫為?0.5 ℃，含冰量為多冰凍土為例，圖5為20年內(nèi)渠道中心、坡腳、坡肩和天然地面等4個位置地溫-深度曲線。為節(jié)約篇幅，圖5中僅給出了渠道運營第5、10、15和20年10月15日最大融化深度時的地溫-深度曲線。

圖5 渠道開挖后20年內(nèi)渠道不同位置地溫曲線Fig. 5 Temperature profiles at different locations within 20 years after the canal construction

從圖5（a）可以看出，渠道中心受渠道內(nèi)流動水體影響較大，第5、10、15和20年10月15日最大融化深度分別為5.2、8.1、10.2和12.2 m，年均增加速率為0.4~0.6 m/a。同時可以看出，最大季節(jié)融化深度的增加主要發(fā)生在前5年，這與渠道開挖的直接影響密切相關(guān)。整體來看，在渠道開挖和過流影響下，渠道下部多年凍土主要表現(xiàn)為自上而下的退化過程，即多年凍土上限下移顯著而下部多年凍土地溫升溫不顯著，這是由于在?0.5 ℃地溫條件下，多年凍土已經(jīng)接近劇烈相變區(qū)。相對于渠中，坡腳的多年凍土退化速率較慢，但其退化模式基本相同（圖5（b））。岸坡和天然場地多年凍土的退化主要受氣候變暖影響（圖5（c）和（d）），尤其在前20年內(nèi)，渠道內(nèi)過水側(cè)向熱侵蝕尚未發(fā)展到岸坡位置?？梢钥闯?，僅在氣候變暖影響下，兩個位置的多年凍土退化相對緩慢，且退化模式均表現(xiàn)為多年凍土上限的下移和上限以下一定深度范圍內(nèi)多年凍土的升溫。通過4個位置地溫的對比，可以很好地反映開挖過水因素和氣候變暖因素作用下多年凍土退化過程的差異性。

4 結(jié) 語

在多年凍土區(qū)，與道路、房建、輸油氣管道等工程有所不同，水利工程涉及大量的挖方工程和流動水體，因此對下伏多年凍土的擾動更加顯著。以高海拔多年凍土區(qū)某渠道工程為例，利用數(shù)值模擬手段，研究了過水和氣候變暖對下部多年凍土熱狀況的影響，并考慮多年凍土地溫、含冰量兩個重要因素，模擬分析了氣候變暖背景下渠道下部多年凍土地溫的時間演化過程，主要結(jié)論如下：

（1）在氣候變暖背景下，天然場地多年凍土的退化主要表現(xiàn)為活動層厚度增加和多年凍土地溫升高。以區(qū)域歷史氣象增溫率模擬氣候變暖過程，當(dāng)多年凍土年平均地溫為?0.5 ℃，含冰量為少冰時，天然場地活動層年均增加速率為3 cm/a，這一量值與目前實測青藏高原活動層增加速率接近。

（2）受渠道開挖和過水影響，渠道底部、坡腳及岸坡位置多年凍土退化顯著。渠道運營50年后，當(dāng)多年凍土含冰量為少冰時，?0.5、?1.0和?1.5 ℃三種多年凍土地溫條件下，渠底、坡腳及天然地面靠水渠的部分下部30m范圍內(nèi)多年凍土已經(jīng)退化，反映了渠道在垂向和水平向強烈的熱侵蝕作用。不同年平均地溫條件下，多年凍土熱惰性不同，因此渠道的側(cè)向熱侵蝕存在一定差異。當(dāng)多年凍土地溫為?0.5 ℃時，自坡肩向外約9.8 m范圍內(nèi)下部多年凍土已退化；而當(dāng)多年凍土地溫為?1.0和?1.5 ℃，岸坡下部仍有多年凍土分布。

（3）含冰量是影響多年凍土退化速率的重要因素。渠道下部多年凍土含冰量由少冰增加至多冰、富冰時，其熱惰性顯著增加，氣候變暖和過水雙重作用下多年凍土退化速率明顯減小。即使在年平均地溫為?0.5 ℃時，渠道運營50年后其下部仍有多年凍土存在。但是，受氣候變暖和過水雙重影響，渠道下部季節(jié)融化深度顯著增加，形成了一個“鍋底狀”的融化盤，融化盤大小與含冰量有關(guān)。

（4）渠底、坡腳、岸坡及天然場地多年凍土地溫-深度的時間變化過程對比表明，在過水和氣候變暖雙重作用下，渠底和坡腳多年凍土的退化主要為多年凍土上限的快速下移，即自上而下的退化模式；而在氣候變暖單獨作用下，岸坡和天然場地多年凍土的退化表現(xiàn)為活動層緩慢下移的同時，上限以下一定升幅范圍內(nèi)多年凍土的緩慢升溫。