孫 哲， 趙 林， 胡國(guó)杰， 喬永平， 杜二計(jì)， 鄒德富， 謝昌衛(wèi)

（1.中國(guó)科學(xué)院西北生態(tài)環(huán)境資源研究院冰凍圈科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室/藏北高原冰凍圈特殊環(huán)境與災(zāi)害國(guó)家野外科學(xué)觀測(cè)研究站，甘肅 蘭州 730000； 2.中國(guó)科學(xué)院大學(xué)，北京 100049； 3.南京信息工程大學(xué)地理科學(xué)學(xué)院，江蘇 南京 210044）

0 引言

多年凍土約占北半球陸地面積的25%［1］，其對(duì)氣候變化尤為敏感［2］。青藏高原是全球海拔最高的地區(qū)，分布面積高達(dá)約1.06×106km2的多年凍土［3］，并以高溫高含冰量為特征［4-5］。全球氣候變暖的背景下，多年凍土呈現(xiàn)退化趨勢(shì)，表現(xiàn)為多年凍土升溫、埋藏深度增大和地下冰融化等，對(duì)寒區(qū)的工程建筑的穩(wěn)定性、氣候以及水文生態(tài)過程都有著重大影響。多年凍土退化導(dǎo)致地基的力學(xué)穩(wěn)定性降低以及引發(fā)的地表沉降是影響寒區(qū)工程穩(wěn)定性和安全性的主要因素［6］；多年凍土與寒區(qū)生態(tài)系統(tǒng)之間也存在著相輔相生的共生關(guān)系，青藏高原多年凍土退化已經(jīng)引起高寒草甸生態(tài)系統(tǒng)發(fā)生變化，甚至出現(xiàn)荒漠化、沙漠化［7-8］；隨著多年凍土消融，凍土中的碳儲(chǔ)量可能被釋放，從而放大高原氣候變暖效應(yīng)［9］；大量賦存于多年凍土中的地下冰轉(zhuǎn)化為液態(tài)水，部分參與到活動(dòng)層的凍融過程中，而其余部分將參與到區(qū)域乃至全球的水循環(huán)過程中［10］；此外，多年凍土退化導(dǎo)致高原下墊面的能量平衡發(fā)生變化，最終可能觸發(fā)整個(gè)北半球氣候系統(tǒng)的變化［11-12］。因此，準(zhǔn)確評(píng)估和預(yù)估全球氣候變暖背景下大時(shí)間尺度多年凍土熱狀況的變化已經(jīng)成為冰凍圈科學(xué)研究的核心環(huán)節(jié)之一［13］。

瞬態(tài)數(shù)值凍土模式是研究多年凍土溫度場(chǎng)變化的一個(gè)有效手段，因?yàn)槠淠莒`活描述多年凍土剖面的異質(zhì)性、幾何形狀和邊界條件［14］。已有多種凍土模式應(yīng)用于模擬多年凍土的分布及變化［15-19］。在現(xiàn)有的凍土模式中，考慮相變的熱傳導(dǎo)方程是計(jì)算多年凍土溫度場(chǎng)的核心。一般的模擬步驟為首先采用初始恒定的上邊界條件進(jìn)行反復(fù)迭代數(shù)百次后，直至多年凍土垂直剖面溫度場(chǎng)達(dá)到穩(wěn)態(tài)作為初始條件，然后賦予相應(yīng)的上邊界條件變化模擬不同情景下凍土溫度場(chǎng)的變化［20］?，F(xiàn)有用于模擬多年凍土的研究主要關(guān)注近地表溫度場(chǎng)變化的影響，主要集中于對(duì)氣溫和淺表層物理過程和參數(shù)化方案等改進(jìn)和優(yōu)化［18，20-21］。而對(duì)于下邊界條件設(shè)置對(duì)模擬多年凍土熱狀況的影響尚缺乏深入評(píng)估。例如大氣環(huán)流模式中的陸面模式一般采用較淺的下邊界，通常為3～10 m，采用零地?zé)嵬炕蚝愣囟茸鳛橄逻吔鐥l件［21-22］，而一些專門針對(duì)多年凍土計(jì)算的數(shù)值模式則采用較深的下邊界，通常大于50 m，采用恒定熱通量或恒定溫度作為下邊界條件［15，19，23-24］。對(duì)深層多年凍土物理過程不恰當(dāng)?shù)谋憩F(xiàn)，是當(dāng)前凍土模式模擬結(jié)果主要的不確定性來源之一［22，25-28］。多年凍土是歷史氣候及地質(zhì)過程的綜合產(chǎn)物，由于土壤是不良的熱導(dǎo)體，熱傳導(dǎo)沿巖土層深度具有滯后性及衰減性。淺層凍土的溫度變化主要受控于年際氣候變化，而隨著巖土層加深，年際氣候波動(dòng)信號(hào)會(huì)被衰減過濾，深層多年凍土的溫度主要取決于長(zhǎng)期歷史氣候以及巖土層下部地?zé)崮芰抗餐饔茫?9］。因此，如果研究只關(guān)注活動(dòng)層或較淺的多年凍土層年際地溫場(chǎng)的變化，也許可以忽視巖土層下部地?zé)崮芰康挠绊?，相反，如果關(guān)注點(diǎn)在于在氣候變暖的背景下，模擬和預(yù)測(cè)百年或千年尺度下多年凍土退化的時(shí)機(jī)和程度，其取決于深層多年凍土地溫場(chǎng)變化，則應(yīng)必須重視下邊界條件設(shè)置的影響［26］。定量界定“淺層”和“深層”多年凍土，同樣需要評(píng)估在所研究的時(shí)間尺度內(nèi)，不同的下邊界條件設(shè)置對(duì)不同深度多年凍土溫度計(jì)算的影響程度。

因此，本文擬利用藏北高原冰凍圈站在青藏高原五道梁地區(qū)架設(shè)的綜合觀測(cè)站的長(zhǎng)期地溫監(jiān)測(cè)資料，基于一維熱傳導(dǎo)凍土模型，以重建過去50 年的地表溫度序列以及未來兩種升溫情景作為上邊界驅(qū)動(dòng)條件，在控制其他變量一致的基礎(chǔ)上，通過設(shè)計(jì)不同的下邊界方案進(jìn)行模擬實(shí)驗(yàn)，以定量評(píng)估下邊界條件對(duì)多年凍土地溫場(chǎng)數(shù)值模擬的影響。分析結(jié)果可為凍土模擬研究中下邊界條件的選取提供一定參考。

1 研究區(qū)與數(shù)據(jù)

青藏高原五道梁地區(qū)（35°13′N、93°05′E，海拔4 783 m）位于昆侖山和唐古拉山脈之間的可可西里地區(qū)的東邊（圖1）。該地區(qū)位于高原北部連續(xù)多年凍土區(qū)，年均氣溫-5.1 ℃，年降水量為284.8 mm，下墊面為蓋度稀疏的高寒荒漠草原。其氣候特征與多年凍土特征有較好的代表性［30-31］。

圖1 研究區(qū)地理位置及五道梁綜合觀測(cè)站照片F(xiàn)ig.1 Location of the study area and the photo showing the comprehensive observation station in Wudaoliang

藏北高原冰凍圈站在五道梁架設(shè)了綜合觀測(cè)站，對(duì)五道梁氣候變化、地氣物質(zhì)能量交換以及多年凍土的動(dòng)態(tài)變化進(jìn)行了長(zhǎng)期的監(jiān)測(cè)及研究［32-41］。凍土監(jiān)測(cè)包括活動(dòng)層水熱監(jiān)測(cè)和鉆孔地溫監(jiān)測(cè)。活動(dòng)層水熱監(jiān)測(cè)采用美國(guó)Campbell 公司的CR1000數(shù)采儀監(jiān)測(cè)5～60 cm 共12 層的土壤溫度及水分，監(jiān)測(cè)時(shí)段2009—2015年。鉆孔深40 m，并未打穿多年凍土，地溫探頭由凍土工程國(guó)家重點(diǎn)實(shí)驗(yàn)室加工，精度為±0.05 ℃，測(cè)量共26 層（0～40 m）的土壤溫度，監(jiān)測(cè)時(shí)段為2005—2015 年。根據(jù)2015 年的監(jiān)測(cè)數(shù)據(jù)，該年的地表平均溫度為-1.77 ℃，活動(dòng)層為2.4 m，地溫年變化深度（ZAA）大約為7 m，年平均地溫（MAGT）為-1.5 ℃。

2 多年凍土數(shù)值模型

2.1 考慮相變的一維熱傳導(dǎo)控制方程

假設(shè)土壤中的熱量傳輸完全由傅里葉定律所決定，忽視熱對(duì)流及熱輻射的影響。模型的熱傳導(dǎo)控制方程表示為

式中：ceff為有效體積熱容量（kJ·m-3·℃-1）；λ為導(dǎo)熱系數(shù)（W·m-1·℃-1）；T為溫度（℃）；t為時(shí)間（s）；z為深度（m）。

通過顯熱容法考慮考慮土壤中的相變過程以及凍土融土熱參數(shù)的差異。ceff與λ為溫度T的分段函數(shù)為

式中：下標(biāo)f與u分別表示土壤的凍結(jié)與融化狀態(tài)；L為單位體積冰水相變潛熱（334.54 MJ·m-3）；Δw為單位體積相變水量；［T1，T2］為相變溫度區(qū)間。

從五道梁觀測(cè)站現(xiàn)場(chǎng)采集的巖芯土壤樣品進(jìn)行實(shí)驗(yàn)分析，主要包括土壤質(zhì)地、干容重和含水量等。土壤樣品通過直徑5 cm、高5 cm 的采樣器從探坑或鉆孔巖芯中采集。土壤質(zhì)地由德國(guó)UGT 科學(xué)儀器公司生產(chǎn)的SEDIMAT4-12 粒度分析儀進(jìn)行分析，其測(cè)量原理為濕篩法結(jié)合移液管法。干容重和含水/冰量采用烘干法（60 ℃，24 h）測(cè)定。模型中所有的土壤熱參數(shù)基于以上測(cè)定的土壤性質(zhì)資料，參考《凍土地區(qū)建筑地基基礎(chǔ)設(shè)計(jì)規(guī)范》選取。模型中土壤參數(shù)取值見表1。

表1 研究區(qū)各巖土層的物理參數(shù)Table 1 Physical parameters of soil layers in the study area

2.2 上邊界條件

為充分考慮歷史氣候?qū)Χ嗄陜鐾聊M的影響，本文提取位于五道梁的國(guó)家氣象站（觀測(cè)數(shù)據(jù)（http：//data.cma.cn/）中相應(yīng)時(shí)間段的日氣溫?cái)?shù)據(jù)，與綜合觀測(cè)場(chǎng)的日地表溫度觀測(cè)數(shù)據(jù)進(jìn)行線性擬合，關(guān)系式為

式中：GST 為五道梁綜合觀測(cè)場(chǎng)的實(shí)測(cè)地表溫度；X為五道梁國(guó)家氣象臺(tái)站的歷史氣溫序列。然后，根據(jù)式（4）與五道梁的國(guó)家氣象站自1965年的歷史日氣溫?cái)?shù)據(jù)序列重建過去50 年（1965—2015 年）的日地表溫度序列（圖2）。

圖2 重建1965—2015年五道梁日地表溫度序列（a）以及地表溫度實(shí)測(cè)值與計(jì)算值的對(duì)比圖（b）Fig.2 Reconstructed ground surface temperature series from 1965 to 2015（a）and scatter plots showing measured vs.calculated ground surface temperature during the observation period（b）

上邊界的地表溫度變化T用下面正弦方程進(jìn)行擬合。

式中：T0為初始地表年均溫度（℃）；k為升溫速率（℃·a-1）；t為時(shí)間（d）；A為地表溫度年較差的一半（℃）；φ為初始相位。根據(jù)重建1965—2015 年的日地表溫度序列的擬合結(jié)果，T0取-2.7 ℃，k取0.025 ℃·a-1，A取11.7 ℃，φ取-1.94。

此外，根據(jù)第5 次耦合模式比較計(jì)劃（CMIP5）11 個(gè)氣候模式的預(yù)測(cè)結(jié)果，2011—2100年中國(guó)氣溫升溫率為0.006～0.063 ℃·a-1［42］。在RCP8.5 情景下，青藏高原地表氣溫最大升溫率最大能達(dá)到0.057 ℃·a-1［43］。在本研究中，為探究百年尺度氣候變化對(duì)模擬多年凍土溫度場(chǎng)變化的影響，假設(shè)兩個(gè)未來（2015—2100 年）升溫情景：一是假設(shè)重建的歷史地溫升溫率（k=0.025 ℃·a-1）在未來一致持續(xù)，二是在RCP8.5 極端氣候升溫情景下假設(shè)地表溫度與氣溫以相同速率上升（k=0.057 ℃·a-1），驅(qū)動(dòng)采用不同下邊界方案的模型模擬2015—2100 年的未來多年凍土的變化（圖3）。

圖3 上邊界地表溫度的升溫趨勢(shì)Fig.3 Rise trend of the surface temperature as the upper boundary ground

2.3 數(shù)值求解

空間步長(zhǎng)設(shè)置為多年凍土垂直剖面5 m 以上取0.05 m，5 m 以下取0.5 m；時(shí)間步長(zhǎng)設(shè)置為1 d。熱傳導(dǎo)偏微分方程通過MATLAB中ode23t函數(shù)（自由內(nèi)插法的梯形法則）求解。

3 實(shí)驗(yàn)設(shè)計(jì)

3.1 9種不同下邊界方案

現(xiàn)有凍土數(shù)值模擬計(jì)算中，下邊界條件的設(shè)置差異主要包括設(shè)置深度不同和設(shè)置方式不同。綜合現(xiàn)有主要的凍土模式對(duì)于下邊界條件的設(shè)置，本文設(shè)計(jì)3個(gè)下邊界深度：10 m（大氣環(huán)流模式中的凍土層一般設(shè)置的范圍［22］，且位于研究區(qū)地溫年變化深度以下）、40 m（研究區(qū)鉆孔底部）和100 m（一些專門針對(duì)百年尺度多年凍土變化數(shù)值模擬研究所采用下邊界深度［17，20］）；3種下邊界條件的設(shè)置方式：零地?zé)嵬浚ú豢紤]地?zé)嵊绊?，第二類邊界條件）、恒定地?zé)崃鳎ǖ诙愡吔鐥l件）和相應(yīng)巖土層的恒定溫度（第一類邊界條件）。3個(gè)深度和3種設(shè)置方式進(jìn)行相互組合，共設(shè)計(jì)9種下邊界條件方案（表2）。根據(jù)五道梁地區(qū)的深孔（大于60 m）的相關(guān)文獻(xiàn)［44-45］，深層地?zé)崽荻燃s為0.03 ℃·m-1，換算成地?zé)崃鳎ǚ桨?～6）。

表2 各方案下邊界條件的設(shè)置Table 2 Lower boundary conditions of the nine projects

3.2 實(shí)驗(yàn)步驟

第一步：根據(jù)9 種不同下邊界方案，設(shè)置9 種多年凍土幾何垂直剖面。

第二步：在上邊界用重建歷史地表溫度序列的第一年地表溫度正弦曲線反復(fù)迭代600～1 000 次（年），使9 種不同方案的多年凍土垂直剖面溫度場(chǎng)達(dá)到穩(wěn)態(tài)（各層溫度的年際變化小于0.001 ℃）。

第三步：以第二步的各多年凍土垂直剖面的穩(wěn)態(tài)溫度場(chǎng)作為初始條件，施加重建歷史地表溫度序列作為上邊界驅(qū)動(dòng)，模擬不同方案的多年凍土過去50年的溫度場(chǎng)變化。

第四步：在第三步模擬結(jié)果的基礎(chǔ)上（所為新的初始條件），分別進(jìn)一步施加k=0.025 ℃·a-1和k=0.057 ℃·a-1兩種未來地表升溫情景，模擬不同方案的多年凍土2015—2100年溫度場(chǎng)變化。

第五步：各方案第三、四步的1965—2100年多年凍土溫度場(chǎng)模擬結(jié)果結(jié)合實(shí)測(cè)值進(jìn)行相互比較分析。

由于本文主要關(guān)注較深的多年凍土地溫場(chǎng)的變化，其年內(nèi)溫度振幅較小或?yàn)榱?。故本文使用平均絕對(duì)誤差（MAE）來判斷模擬值與基準(zhǔn)值之間的偏差。

式中：Bi為基準(zhǔn)值；Si為模擬值。

4 結(jié)果與分析

4.1 各方案模擬結(jié)果及與實(shí)測(cè)值的對(duì)比

圖4 展示各方案在實(shí)驗(yàn)第二步預(yù)熱結(jié)束后多年凍土垂直剖面達(dá)到穩(wěn)態(tài)的地溫曲線-深度（1965年）、第三步模擬至2015 年與實(shí)測(cè)的2015 年地溫-深度曲線，以及第四步在兩種未來升溫情景下模擬至2100 年（末年）地溫-深度曲線。圖5 進(jìn)一步詳細(xì)展現(xiàn)各方案模擬值與實(shí)測(cè)值的對(duì)比。

圖4 各方案模擬的1965年（預(yù)熱完成）的地溫-深度曲線（a1，b1，c1），模擬和實(shí)測(cè)的2015年地溫-深度曲線（a2，b2，c2），在k=0.025 ℃·a-1（a3，b3，c3）和k=0.057 ℃·a-1（a4，b4，c4）未來升溫情景下模擬的2100年（末年）地溫-深度曲線Fig.4 Comparison of simulated 1965（initialization complete）temperature-depth curve（a1，b1，c1），simulated and measured 2015 temperature-depth（a2，b2，c2），and simulated 2100（last year）temperature-depth curve as k=0.025 ℃·a-1（a3，b3，c3）and k=0.057 ℃·a-1（a4，b4，c4）

圖5 各方案3 m（a1，b1，c1）、5 m（a2，b2，c2）、10 m（a3，b3，c3）、16 m（a4，b4，c4）和30 m（a5，b5，c5）地溫模擬結(jié)果與實(shí)測(cè)對(duì)比（圖中MEA值計(jì)算以實(shí)測(cè)值為基準(zhǔn)）Fig.5 Comparison of observation and simulation of the ground temperature in 3 m（a1，b1，c1），5 m（a2，b2，c2），10 m（a3，b3，c3），16 m（a4，b4，c4）and 30 m（a5，b5，c5）from each project（the MAE in the figure calculated based on the observation value）

方案1～3不考慮來自巖土層下部地?zé)崮芰坑绊?，采用零地?zé)嵬孔鳛橄逻吔鐥l件。在預(yù)熱結(jié)束后方案1～3計(jì)算的1965年地度-深度曲線基本重合［圖4（a1）］，而隨著地表升溫，至2015年，方案1～3的多年凍土地溫模擬結(jié)果相對(duì)于實(shí)測(cè)值明顯偏小，且下邊界設(shè)置越深，偏小越顯著［圖4（a2），圖5（a1～a5）］。原因?yàn)閷?shí)驗(yàn)第二步的預(yù)熱計(jì)算的初始穩(wěn)態(tài)地溫場(chǎng)的物理意義為在之前百年千年長(zhǎng)時(shí)間尺度下處于恒定的地氣熱交換所達(dá)到穩(wěn)定狀態(tài)的過程，實(shí)驗(yàn)第三步反映了氣候變暖改變了上邊界條件導(dǎo)致多年凍土溫度場(chǎng)的變化。因此方案1～3在實(shí)驗(yàn)第二步和第三步中忽視了來自巖土層下部地?zé)崮芰繉?duì)多年凍土地溫場(chǎng)的加熱作用，模擬的地溫場(chǎng)會(huì)顯著偏小。而下邊界設(shè)得越深，更深層的巖土層將吸收更多來自地表升溫的能量，故模擬深層巖土層的溫度就越偏?。?2］。

從相變熱傳導(dǎo)理論的角度，采用恒定地?zé)崃髯鳛榈诙愊逻吔鐥l件的方案4～6相對(duì)更符合實(shí)際的物理過程［14，26］。在預(yù)熱結(jié)束后方案4～6計(jì)算的1965年地溫-深度曲線也基本重合［圖4（b1）］，而隨著地表升溫，方案4～6 的多年凍土2015 年地溫模擬結(jié)果能較好的符合實(shí)測(cè)值［圖4（b2）］，也能較好的再現(xiàn)各巖土層實(shí)測(cè)近十年的升溫情況［圖5（b1～b5）］。就方案4～6相互比較而言，下邊界深度設(shè)置越淺，多年凍土地溫模擬結(jié)果就會(huì)越偏高。這是因?yàn)樵趯?shí)驗(yàn)第三步地溫自上而下逐步升高時(shí)，上層的地溫梯度會(huì)減小，同時(shí)又受到巖土層下部地?zé)崃髯饔茫沟谜麄€(gè)地溫曲線隨時(shí)間呈現(xiàn)“弓形”［46］。實(shí)測(cè)的上層多年凍土垂直剖面的地溫梯度小于深層的垂直剖面地溫梯度。若下邊界深度設(shè)得淺，又采用深層的地溫梯度換算的恒定地?zé)崃髯鳛橄逻吔鐥l件，則會(huì)高估來自巖土層下部地?zé)崃鬏斎氲哪芰?，進(jìn)而模擬出偏高的多年凍土溫度。

采用2015 年恒定地溫作為第一類下邊界條件的方案7～9 的模擬結(jié)果能很好的符合2015 年實(shí)測(cè)的地溫場(chǎng)［圖4（c1）］。但由于固定了下邊界的溫度，其無法反映下邊界附近巖土層的升溫情況［圖5（c1～c5）］。方案8 把下邊界設(shè)在40 m，方案9 設(shè)為在100 m，與各巖土層實(shí)測(cè)的近十年升溫速率相比，方案8 中模擬8 m 以下巖土層的升溫速率相對(duì)于方案9與實(shí)測(cè)值開始偏小，越往深的巖土層，偏小幅度越大；而將下邊界設(shè)置與10 m 的方案7 則無法模擬5 m 以下實(shí)測(cè)的升溫變化。特別是上邊界有明顯的升溫或降溫的情景時(shí)，定溫邊界會(huì)對(duì)模擬結(jié)果有很大影響。如圖4（c3）和圖4（c4）所示，當(dāng)模擬未來百年尺度多年凍土地溫場(chǎng)變化時(shí)，無論是k=0.025 ℃·a-1還是k=0.057 ℃·a-1升溫情景，方案9 顯示10 m和40 m 的巖土層溫度會(huì)發(fā)生較大的變化，故采用較淺深度的定溫邊界條件的方案7 和8 所預(yù)估的結(jié)果是不具備參考價(jià)值的。

4.2 以方案6為基準(zhǔn)各方案的相互對(duì)比

基于研究目的，本文近地表層的各種陸面過程考慮得比較簡(jiǎn)單，使得各方案的近地表層溫度模擬結(jié)果以實(shí)測(cè)值作為基準(zhǔn)值時(shí)MAE 較大。方案6 在100 m 深度以恒定地?zé)崃髯鳛橄逻吔鐥l件，相比其他方案更接近實(shí)際的巖土層下部能量傳輸物理過程，其模擬結(jié)果也能更好的重現(xiàn)研究區(qū)實(shí)測(cè)的多年凍土溫度場(chǎng)。為從理論上定量評(píng)估不同下邊界設(shè)置對(duì)計(jì)算多年凍土溫度場(chǎng)的影響程度，以方案6 的模擬結(jié)果作為基準(zhǔn)值（理論基準(zhǔn)值），定量比較其他方案不同深度巖土層溫度模擬結(jié)果相對(duì)于方案6 的MAE（表3），同時(shí)也參考各方案模擬結(jié)果以實(shí)測(cè)為基準(zhǔn)值時(shí)的MAE（圖5），從熱傳導(dǎo)的角度對(duì)多年凍土的垂直剖面進(jìn)行分層。初步設(shè)定3 個(gè)MAE 的閾值：0.1 ℃、0.5 ℃以及1 ℃。根據(jù)目前使用的溫度探頭的精度為±0.05 ℃，以及多年凍土的工程分類：極不穩(wěn)定（＞-0.5 ℃），不穩(wěn)定（-1～-0.5 ℃），穩(wěn)定過渡（-3～-1 ℃）和穩(wěn)定（＜-3 ℃）；將0.1 ℃的MAE 作為在儀器誤差范圍內(nèi)判定模擬值與基準(zhǔn)值是否有偏差的閾值，將0.5 ℃的MAE 作為判定模擬值與基準(zhǔn)值是否有顯著偏差的閾值（偏差大于0.5 ℃就可能誤判多年凍土分類），將1 ℃作為模擬值與基準(zhǔn)值是否有極顯著偏差的閾值。

表3 以方案6為基準(zhǔn)值，計(jì)算的各方案歷史模擬結(jié)果的平均絕對(duì)誤差Table 3 The MAE of historic simulation from each project calculated based on the sixth project as the truth value

如表3 所示，以方案6 為基準(zhǔn)值時(shí)，其他方案的3 m 以上巖土層溫度場(chǎng)模擬結(jié)果的MAE 均小于0.1 ℃，由于所有方案的上邊界條件設(shè)置相同，說明3 m 以上巖土層溫度變化基本完全由上邊界年際氣候變化所決定，從熱傳導(dǎo)的角度可把3 m 以上的巖土層劃為“近地表層”，其不受巖土層下部的地?zé)崮芰坑绊?。?duì)于不考慮地?zé)嵊绊懙姆桨?～3，其3～15 m 巖土層溫度模擬結(jié)果的MAE 在0.1～0.5 ℃之間（若以實(shí)測(cè)為基準(zhǔn)值，MAE 在0.27～0.44 ℃之間），說明該深度范圍的巖土層溫度變化雖然還是主要受控于年際氣候變化，但隨著深度增加，下邊界設(shè)置的影響開始逐漸顯現(xiàn)，該深度范圍的巖土層可劃為“淺層”；15～30 m 巖土層溫度模擬結(jié)果的MAE 在0.5～1 ℃之間（若以實(shí)測(cè)為基準(zhǔn)值，MAE 在0.44～0.97 ℃之間），說明該深度范圍的巖土層溫度變化已經(jīng)開始受到巖土層下部地?zé)崮芰康娘@著影響。青藏高原的多年凍土地溫年變化深度通常在12～20 m 之間［47］，計(jì)算15～30 m 的多年凍土溫度時(shí)，如果不考慮地?zé)嵊绊?，可能?huì)導(dǎo)致較大的誤差進(jìn)而誤判多年凍土的類型，該深度范圍可劃為“中層”；30 m 以下巖土層的溫度模擬結(jié)果的MAE 已經(jīng)大于1 ℃（若以實(shí)測(cè)為基準(zhǔn)值，方案3 的MAE 同樣大于1 ℃），說明使用不考慮地?zé)嵊绊懙姆桨赣?jì)算30 m以下巖土層的溫度應(yīng)當(dāng)是不可取的。因此，30 m 以下巖的深度范圍可以劃為“深層”。

對(duì)于同樣以恒定地?zé)崃髯鳛榈诙愊逻吔鐥l件的方案4和5，下邊界設(shè)置10 m深度的方案4的淺層地溫模擬結(jié)果的MAE 在0.1～0.5 ℃之間；下邊界設(shè)置40 m 深度的方案5 的深層地溫模擬結(jié)果的MAE 開始大于0.1 ℃（若以實(shí)測(cè)為基準(zhǔn)值，則20 m以下巖土層的MAE開始大于0.1 ℃）。

對(duì)于以恒定溫度作為第一類下邊界條件的方案7～9，同樣將下邊界設(shè)置在100 m 處的方案9 與作為基準(zhǔn)的方案6 的模擬結(jié)果幾乎完全一致（MAE＜0.01 ℃）；將下邊界設(shè)置在40 m 處的方案8中，20 m以下巖土層溫度模擬結(jié)果開始顯現(xiàn)較小的偏差（0.1 ℃＜MAE＜0.31 ℃）；而將下邊界設(shè)置在10 m處的方案7 中，則3 m 以下（淺層范圍內(nèi)）地溫模擬結(jié)果開始出現(xiàn)偏差（MAE＞0.1 ℃），7 m 以下地溫模擬結(jié)果的MAE大于0.3 ℃。

4.3 百年尺度氣候變化對(duì)多年凍土垂直剖面溫度變化的影響

下邊界深度應(yīng)設(shè)置在地?zé)崃骰驕囟缺容^穩(wěn)定的巖土層，即在所研究的時(shí)間尺度內(nèi)，上邊界氣候變化難以影響到的巖土層。目前對(duì)于氣候變暖的背景下青藏高原多年凍土變化的模擬及預(yù)估的研究，主要還是集中于近百年尺度內(nèi)［48］。因此，下面探討在兩種未來升溫情景下，不同下邊界條件的設(shè)置對(duì)百年尺度多年凍土溫度場(chǎng)模擬的影響。

表4和表5分別展示各方案模擬在k=0.025 ℃·a-1和k=0.057 ℃·a-1兩種未來升溫情景下不同深度的巖土層未來近100年來總升溫量?？梢悦黠@看出恒定溫度（第一類邊界條件）作為下邊界條件的方案7、8的對(duì)靠近下邊界巖土層的升溫變化模擬的局限性。若以方案6為基準(zhǔn)，在兩種未來升溫情景下，70 m以下巖土層的升溫總量小于0.01 ℃，這也解釋了采用100 m 恒定溫度作為下邊界條件的方案9 與方案6 模擬結(jié)果基本一致的原因。但是方案3 的模擬結(jié)果顯示，70 m 以下巖土層的升溫總量均大于1 ℃，百年尺度內(nèi)的氣候變化依舊會(huì)影響這些深層巖土層的溫度。在k=0.057 ℃·a-1升溫情景下，方案6 和9 的模擬結(jié)果顯示升溫的近100 年間多年凍土消融量（多年凍土厚度減小）為20.15 m，而方案3顯示升溫的近100 年間多年凍土消融量為15.20 m。這說明百年尺度上邊界變化對(duì)深層多年凍土溫度的影響不僅與氣候變化的幅度有關(guān)，而且與多年凍土消融產(chǎn)生的相變熱的多少有也直接的關(guān)系，氣候變化信號(hào)導(dǎo)致傳到深層導(dǎo)致相應(yīng)位置的負(fù)熱流（熱流方向以向上為正）發(fā)生了改變，進(jìn)而引起巖土層的溫度發(fā)生變化，若多年凍土相變耗熱多，則傳到此處的負(fù)熱流就較小，所以深層地溫變化小，如果相變熱小，深層地溫可能就會(huì)有顯著的變化。

表4 各方案模擬在k=0.025 ℃·a-1升溫情景下多年凍土垂直剖面2015—2100年總升溫量Table 4 Simulated warming of the geothermal field from the nine projects during 2015—2100 as k=0.025 ℃·a-1

表5 各方案模擬在k=0.057 ℃·a-1升溫情景下多年凍土垂直剖面2015—2100年總升溫量Table 5 Simulated warming of the geothermal field from the nine projects during 2015—2100 as k=0.057 ℃·a-1

各方案預(yù)估結(jié)果的比較顯示，當(dāng)上邊界有長(zhǎng)時(shí)間尺度顯著升溫的情景時(shí)，不同的下邊界設(shè)置會(huì)對(duì)多年凍土消融及其產(chǎn)生的相變熱的計(jì)算會(huì)有顯著影響，進(jìn)而模擬結(jié)果中百年尺度氣候變化對(duì)深層多年凍土地溫的影響程度也會(huì)有很大的差異。如圖4（a4）和圖4（b4）所示的結(jié)果，模擬RCP8.5 升溫情景（k=0.057 ℃·a-1）下多年凍土百年尺度下未來變化時(shí)，若采用第二類邊界條件且下邊界設(shè)在淺層的方案1和方案4的模擬結(jié)果顯示到2100年五道梁地區(qū)的多年凍土以完全消融，變?yōu)榉嵌嗄陜鐾恋貐^(qū)；然而若下邊界位置設(shè)置在深層以下的方案2～5的模擬結(jié)果顯示到2100 年五道梁地區(qū)依舊存在著殘余多年凍土，不能劃為非多年凍土區(qū)。目前氣候模式中計(jì)算的多年凍土層設(shè)置很薄的情況，無疑會(huì)對(duì)氣候變暖的背景下全球多年凍土面積變化的預(yù)估結(jié)果造成影響并引發(fā)較大的爭(zhēng)議［26］。

5 討論與結(jié)論

本文基于一維熱傳導(dǎo)凍土溫度數(shù)值模型，以青藏高原五道梁地區(qū)多年凍土為背景，設(shè)計(jì)9 種不同下邊界條件設(shè)置方案進(jìn)行控制模擬實(shí)驗(yàn)，定量評(píng)估在百年尺度氣候變化下不同下邊界條件設(shè)置對(duì)多年凍土溫度場(chǎng)模擬的影響。結(jié)果顯示，五道梁地區(qū)近表層（＜3 m）的地溫完全受控于年際氣候變化，淺層（3～15 m）和中層（15～30 m）的多年凍土溫度場(chǎng)受下邊界條件的影響逐漸顯著，下邊界條件的設(shè)置對(duì)計(jì)算百年尺度氣候變化下多年凍土消融相變熱以及深層（＞30 m）的溫度場(chǎng)具有重大影響。不同多年凍土地區(qū)由于下墊面、地層性質(zhì)和地表溫度振幅等差異，可能導(dǎo)致具體的分層情況有所差異，限于篇幅原因，本文未探討這些因素對(duì)下邊界條件設(shè)置的影響，這也是下一步研究的關(guān)注點(diǎn)。

最后根據(jù)不同的條件與目的，對(duì)于采用基于一維熱傳導(dǎo)方程為核心的凍土模式對(duì)多年凍土溫度場(chǎng)進(jìn)行模擬的研究，本文給出如下的下邊界條件設(shè)置建議：

（1）在條件允許的情況，建議采用實(shí)測(cè)的深層恒定地?zé)崃髯鳛橄逻吔鐥l件且下邊界設(shè)置在深層，使得多年凍土溫度場(chǎng)模擬結(jié)果既能反映巖土層下部的地?zé)崮芰康挠绊懹帜芊从吃跉夂蛴绊懴碌淖兓内厔?shì)。

（2）近表層的地溫完全受控于年際氣候變化，而近表層和淺層的熱流為負(fù)值，且隨深度遞減、周期波動(dòng)，在這個(gè)深度范圍內(nèi)采用單一熱流邊界并不合適。如果研究關(guān)注當(dāng)前研究區(qū)活動(dòng)層厚度或年際活動(dòng)層溫度場(chǎng)變化，考慮計(jì)算效率，可以忽略巖土層下部地?zé)崮芰坑绊?，建議下邊界采用零地?zé)嵬浚O(shè)置在中層或地溫年變化深度以下。

（3）如果研究關(guān)注十年或十年以上多年凍土年平均溫度（MAGT）的變化，應(yīng)當(dāng)不可忽視來自巖土層下部地?zé)崮芰?，建議下邊界采用恒定地?zé)崃?，設(shè)置在中層或深層。恒定溫度（第一類邊界條件）作為下邊界條件時(shí)，無法反映下邊界位置附近的多年凍土溫度的變化。因此，僅僅對(duì)于地溫年際沒有變化或變化很小的情景，恒溫邊界才是可行的。對(duì)于處于零地?zé)崽荻认嘧冞^渡型多年凍土［49］，則建議下邊界設(shè)置在多年凍土底板以下的融土層。

（4）如果研究關(guān)注模擬及預(yù)估世紀(jì)尺度或更長(zhǎng)時(shí)間尺度，多年凍土隨時(shí)間退化甚至完全消失過程，則應(yīng)慎重考慮下邊界條件設(shè)置對(duì)多年凍土消融相變熱以及深層多年凍土溫度場(chǎng)計(jì)算的影響，建議下邊界設(shè)置在深層或多年凍土底板以下的融土層，采用恒定地?zé)崃髯鳛橄逻吔鐥l件。