黃帥 丁琳 丁黔 胡爭 王鶴 陳克政

(黑龍江大學(xué)，哈爾濱，150000) (黑龍江省生態(tài)地質(zhì)調(diào)查研究院) (東北林業(yè)大學(xué))

中國多年凍土主要分布在東北大小興安嶺及西部青藏高原地區(qū)。據(jù)統(tǒng)計，現(xiàn)存中國多年凍土總面積為1.59×106km2，其中東北高緯多年凍土面積為0.24×106km2(約占15%)[1]。受氣候變化、人類活動的影響，東北多年凍土正在發(fā)生快速的變化[2-4]。而隨著基礎(chǔ)設(shè)施建設(shè)的發(fā)展，諸如高鐵、管道、隧道等工程項(xiàng)目在中國東北地區(qū)逐步開展，建立在凍土區(qū)上的寒區(qū)工程也備受快速變化的凍土所影響，因此開展區(qū)域凍土熱狀態(tài)的模擬，對寒區(qū)工程建設(shè)、生態(tài)環(huán)境變化評估具有重要意義[5-7]。

對中國東北部永久凍土的研究，主要依賴于地質(zhì)調(diào)查、公路和鐵路以及中俄輸油管道等工程走廊沿線的研究[8-13]，不能完全反映中國東北部多年凍土的空間變化和異質(zhì)性。隨著凍土學(xué)科與遙感、GIS等技術(shù)快速發(fā)展，凍土模型、凍土熱狀態(tài)模擬方法發(fā)展迅速，計算方法總體上分為經(jīng)驗(yàn)?zāi)Ｐ?、平衡模型、?shù)值模型(物理模型)[14-16]。

在東北大興安嶺地區(qū)，由于其地理特性，目前開展該區(qū)域的凍土分布模擬研究較少。由于該區(qū)域地形地貌復(fù)雜，難以開展大范圍的地質(zhì)調(diào)查，因此該區(qū)域?qū)τ趦鐾恋乃疅帷⑽锢韰?shù)累計較少。凍結(jié)數(shù)模型與等效緯度模型，曾被應(yīng)用與模擬東北地區(qū)的凍土分布[7，17]，研究主要依據(jù)氣象站提供的地表、大氣溫度數(shù)據(jù)對該區(qū)域的凍土參數(shù)進(jìn)行統(tǒng)計回歸，得到經(jīng)驗(yàn)?zāi)Ｐ停Ｐ途群艽蟪潭壬鲜艿搅藲庀笳緜€數(shù)的影響。本研究在大興安嶺地區(qū)內(nèi)沿北向南選取4個試驗(yàn)點(diǎn)(二十四站、塔爾根、臨海、古源)，應(yīng)用庫德里亞夫采夫(Kudryavtsev)模型(平衡模型)、一維瞬態(tài)傳熱的數(shù)值模型(物理模型)對研究區(qū)多年凍土活動層厚度與熱狀態(tài)進(jìn)行模擬；對比實(shí)測地溫與模型模擬地溫，檢驗(yàn)?zāi)Ｐ湍M精度；依據(jù)模擬精度、計算成本，比較2種模型的實(shí)際應(yīng)用價值。旨在為未來區(qū)域凍土分布模擬提供參考。

1 研究區(qū)概況

大興安嶺位于黑龍江省北部地區(qū)，該地區(qū)具有高緯度、低海拔、低溫時間漫長、植被覆蓋率高、積雪覆蓋持久等地理特征[18]，廣泛分布著興安-貝加爾型多年凍土(生態(tài)敏感型凍土)，凍土的分布特征具有顯著的緯度地帶性(見圖1)，局地因子也一定程度地影響著該地區(qū)多年凍土的發(fā)育。受氣候變化的影響，近50 a來大興安嶺地區(qū)多年凍土正在發(fā)生著快速變化[19-20]。

圖1 研究區(qū)凍土分布

受凍土退化的影響，大興安嶺地區(qū)凍土災(zāi)害頻發(fā)，諸如滑坡、崩塌、熱融湖塘擴(kuò)展、熱融喀斯特等[21-22]。在退化過程中，土壤中有機(jī)碳逐漸分解并釋放，地表水熱平衡遭到破壞，影響了地表植被演替，也加速了氣候變暖的進(jìn)程[2]。

本研究在大興安嶺地區(qū)內(nèi)設(shè)置4個試驗(yàn)點(diǎn)，沿北向南分別為二十四站、塔爾根、臨海、古源(見表1)，從緯度地帶性看，從北向南年均氣溫逐漸降低，但年均地溫卻受到了局地因子影響，如二十四站年均地溫高于南部站點(diǎn)的地溫。

表1 試驗(yàn)點(diǎn)年均氣溫與年均地溫(2015—2017年)

試驗(yàn)點(diǎn)附近，僅有新林氣象站(122°31′E，52°58′N)與試驗(yàn)點(diǎn)相鄰。該氣象站記錄的2014—2015年的氣溫與地表溫度(見圖2)，2014年年均氣溫為-2.3 ℃、年均地表溫度為4.9 ℃。由圖2可見：冬季地表溫度顯著高于氣溫，說明植被與積雪具有較好的保溫效果，局地因子影響效果顯著。

圖2 新林氣象站記錄的2014—2015年氣溫與地表溫度

2 研究方法

2.1 一維傳熱數(shù)值模型

一維傳熱模型是依據(jù)熱擴(kuò)散偏微分方程建立的物理模型[16]，模型構(gòu)建見式(1)～(3)。

ρ·C(?T/?t)=(?/?x)[λ(?T/?x)]+

(?/?y)[λ(?T/?y)]；

(1)

(2)

(3)

式中：C為考慮凍土相變的等效比熱(單位為J·kg-1·℃-1)；ρ為土的天然密度(單位為kg·m-3)；λ為土的導(dǎo)熱系數(shù)(單位為W·m-1·℃-1)；Cu、Cf分別為融土及凍土的比熱容(單位為J·kg-1·℃-1)；L為水的相變潛熱(單位為J·kg-1)；λu、λf分別為融土及凍土的導(dǎo)熱系數(shù)(單位為W·m-1·℃-1)；W、Wt分別為凍土的總含水量及含冰量(用百分比表示)；Tp、Tb分別為凍土劇烈相變區(qū)的上、下界溫度(單位為℃)；T為溫度變量(單位為℃)；t為時間變量(單位為s)；x為沿某方向的空間變量(單位為m)。

未凍水含水量根據(jù)式(4)確定：

Wt=a|T|-b。

(4)

式中：a、b為量綱為1的常數(shù)，根據(jù)不同土質(zhì)類型由擬合經(jīng)驗(yàn)公式確定。

模型邊界條件以地表上2 m處氣溫作為輸入，見式(5)：

T(t)=Ta+A0sin[(2π/8 760)t+π/2]。

(5)

式中：Ta為年平均地表溫度(單位為℃)；A0為氣溫年較差(單位為℃)。

2.2 庫德里亞夫采夫模型

庫德里亞夫采夫(Kudryavtsev)模型是一種常用的平衡方程模型，用于計算凍土的活動層厚度與多年凍土頂板溫度[16]。計算方法見式(6)～式(8)。

(6)

(Ts+QL/2Cu)]}-QL/2Cu；

(7)

(2AZCu+QL)。

(8)

式中：Z為最大季節(jié)凍結(jié)深度或最大季節(jié)融化深度(單位為m)；Ts為年均地表溫度(單位為℃)；As為地表溫度年振幅(單位為℃)；τ為周期(8 640 h)；QL為土體中水的相變潛熱(單位為J·m-3)；Az為季節(jié)凍結(jié)(融化)層內(nèi)溫度的平均較差(單位為℃)；ZC單位為m，無實(shí)際物理意義。

2.3 模型輸入?yún)?shù)

2.3.1 土層參數(shù)

根據(jù)現(xiàn)場鉆孔取樣，得到的地層土的工程分類見圖3。

圖3 試驗(yàn)點(diǎn)地層信息

土體熱參數(shù)根據(jù)試驗(yàn)數(shù)據(jù)與文獻(xiàn)進(jìn)行賦值[23]，未凍水參數(shù)(a、b)根據(jù)文獻(xiàn)給出的參數(shù)進(jìn)行賦值[24]，土體參數(shù)見表2。

表2 研究區(qū)土體的物理參數(shù)

2.3.2 一維傳熱模型的邊界條件

試驗(yàn)點(diǎn)地表溫度由現(xiàn)場埋設(shè)傳感器監(jiān)測所得(見圖4，2015年)；氣溫數(shù)據(jù)采用美國國家航空航天局(NASA)提供的氣溫再分析數(shù)據(jù)產(chǎn)品，時間分辨率為8 d、空間分辨率為1 km(見圖5)；應(yīng)用三角函數(shù)擬合獲得年均氣溫與氣溫振幅。

圖4 試驗(yàn)點(diǎn)地表溫度

圖5 試驗(yàn)點(diǎn)氣溫

3 結(jié)果與分析

3.1 模型計算值與鉆孔監(jiān)測值對比

根據(jù)試驗(yàn)點(diǎn)的鉆孔實(shí)測地溫(見圖6)，通過0 ℃等溫線法獲得試驗(yàn)點(diǎn)的活動層厚度與季節(jié)凍結(jié)深度，并計算年均地溫(見表3)。

圖6 試驗(yàn)點(diǎn)監(jiān)測地溫

表3 試驗(yàn)點(diǎn)活動層厚度與年均地溫

利用一維熱傳導(dǎo)模型計算試驗(yàn)點(diǎn)凍土地溫，對比實(shí)測與模擬值的最大季節(jié)凍結(jié)深度月份與最大季節(jié)融化深度月份的地溫(見圖7)，發(fā)現(xiàn)一維熱傳導(dǎo)數(shù)值模型對試驗(yàn)點(diǎn)的地溫模型效果較好，預(yù)測值與實(shí)測值趨勢基本一致。

圖7 試驗(yàn)點(diǎn)凍土地溫的熱傳導(dǎo)模型模擬值與實(shí)測值

根據(jù)一維傳熱模型模擬的地溫曲線，可繪制地溫等值線圖(見圖8)；根據(jù)0 ℃最低位置可得到最大季節(jié)凍深與最大季節(jié)融深，對比庫德里亞夫采夫模型與一維傳熱模型活動層計算值(見表4)。

表4 試驗(yàn)點(diǎn)最大季節(jié)凍深和最大季節(jié)融深的模擬值與實(shí)測值

圖8 塔爾根地溫等值線

3.2 模型誤差分析

對比庫德里亞夫采夫模型、一維傳熱模型對試驗(yàn)點(diǎn)的活動層厚度與最大季節(jié)凍結(jié)深度模擬結(jié)果，一維傳熱模型的計算誤差在0.14～0.18 m、庫德里亞夫采夫模型的計算誤差在0.19～0.30 m，一維傳熱模型的模擬精度優(yōu)于庫德里亞夫采夫模型。

模型的計算誤差來源于多個方面。對于一維傳熱模型，一方面模型本身只考慮了傳熱問題，而活動層是水熱遷移最為敏感的一層土層，由于水熱遷移作用會影響熱土體中的熱傳導(dǎo)與熱對流，從而對模擬結(jié)果造成誤差；另一方面，在地氣熱對流過程中，冬季的植被與積雪具有保溫層的作用，使地表溫度高于氣溫，夏季的植被與積雪阻礙了地表向大氣的熱傳遞，因此對于試驗(yàn)點(diǎn)中地表為裸地的古源南試驗(yàn)點(diǎn)，其模擬誤差最小。

對于庫德里亞夫采夫模型，一方面，在模型計算時，需假設(shè)地層為勻質(zhì)土，而實(shí)際場地的地層往往較為復(fù)雜，這會造成一定的誤差；另一方面，模型缺乏對局地因子影響參數(shù)，大興安嶺地區(qū)多年凍土受生態(tài)環(huán)境影響顯著，諸如坡向、坡度、積雪、植被、海拔等因素均會對凍土賦存狀態(tài)造成影響，若不考慮局地因子效應(yīng)，則會造成一定的模型計算誤差。

4 討論與結(jié)論

對比一維傳熱模型、庫德里亞夫采夫模型，在計算成本方面，庫德里亞夫采夫模型顯著優(yōu)于一維傳熱模型，這是由于模型本身構(gòu)成決定的；傳熱模型是依據(jù)偏微分方程構(gòu)建的，而庫德里亞夫采夫模型則是依據(jù)平衡等式建立的，其求解效率高于偏微分方程，但在計算精度方面，物理模型顯著高于平衡方程模型。另一方面，物理模型能夠?qū)鐾恋牡販氐膭討B(tài)發(fā)展過程進(jìn)行全程模擬，能夠反映出更多的凍土熱狀態(tài)參數(shù)，諸如年均地溫、凍土頂板溫度、多年凍土層厚度等。

在模型深入研究過程中，一方面可以考慮將水分場加入傳熱模型，構(gòu)建凍土的水熱耦合模型；另一方面，需要增加對于地表水熱參數(shù)(如：地表溫度、土壤含水量等)的時空變化反演的問題研究。

本研究應(yīng)用一維傳熱模型、庫德里亞夫采夫模型計算了大興安嶺地區(qū)4個試驗(yàn)點(diǎn)的凍土活動層厚度、最大季節(jié)凍結(jié)深度與地溫，結(jié)果表明：一維傳熱模型對4個試驗(yàn)點(diǎn)的凍土地溫模擬值與實(shí)測值具有一致的趨勢，模擬結(jié)果較好，對于活動層厚度與最大季節(jié)凍結(jié)深度的模擬，具有良好的精度，誤差范圍在0.14～0.18 m之間，能夠較好地滿足工程計算的需求。庫德里亞夫采夫模型在計算活動層厚度與季節(jié)凍結(jié)深度時，模型計算誤差在0.19～0.30 m，基本能夠滿足工程計算需求。在計算成本方面，庫德里亞夫采夫模型優(yōu)于一維傳熱模型，該模型能夠較好地應(yīng)用于工程計算，提高工程效率。