李銳超 謝瑾博 謝正輝

1 中國科學(xué)院大氣物理研究所大氣科學(xué)和地球流體力學(xué)數(shù)值模擬國家重點(diǎn)實(shí)驗(yàn)室，北京 100029

2 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院，北京 100049

1 引言

凍土在北半球分布廣泛，約占其陸地面積的25%（Muller et al., 1947; Zhang et al., 1999）。由于氣候變暖，近年來多年凍土退化嚴(yán)重、反映地殼表層每年寒季凍結(jié)、暖季融化巖土層的活動層厚度不斷增加（Delisle, 2007; Koven et al., 2013）。近 20年來青藏高原多年凍土相對面積減少4.3%，約3.1×104km2，相應(yīng)的季節(jié)凍土面積增加 3.1×104km2（ 張 國 勝 等 , 2007; Yang et al., 2010; 徐 曉 明 等 ,2017）。多年凍土活動層厚度受氣候條件的驅(qū)動，并與海拔、經(jīng)緯度、活動層巖性、含水特征、地面溫度以及下墊面性質(zhì)有關(guān)。另一方面，凍土活動層厚度變化影響陸氣之間的水分能量交換以及生態(tài)系統(tǒng)的碳循環(huán)過程，對天氣以及氣候過程產(chǎn)生重要影響。隨著活動層厚度的增加，北極多年凍土中儲存著大量的二氧化碳、甲烷等溫室氣體會釋放到大氣中，進(jìn)一步加劇全球變暖（Schuur et al., 2015）。世界氣候研究計(jì)劃WCRP已經(jīng)將活動層厚度的變化研究列為氣候和冰凍圈CLIC計(jì)劃的主要觀測研究內(nèi)容之一（Brown et al., 2000; Nelson et al., 2004）。這些都表明研究多年凍土區(qū)活動層厚度變化及其影響的重要性。

早期開展的環(huán)北極多年凍土區(qū)活動層厚度以及近地表凍土變化的監(jiān)測表明：在過去幾十年里，北極地區(qū)的活動層厚度不斷增加，其中俄羅斯北極地區(qū)活動層厚度增加較多，增加了約22 cm（Harris et al., 2003; Frauenfeld et al., 2004）。近幾十年來青藏高原多年凍土活動層厚度也表現(xiàn)出明顯增加的趨勢，對水文過程、生態(tài)系統(tǒng)碳循環(huán)過程、氣候變化、寒區(qū)工程以及基礎(chǔ)設(shè)施等影響巨大（Zhao et al.,2004; 李元壽等, 2010; Wu and Zhang, 2010; 張中瓊和吳青柏, 2012）。多年凍土區(qū)主要分布在北半球高緯度以及高山地區(qū)。在大尺度空間和時間上活動層厚度觀測數(shù)據(jù)稀少，利用陸面過程模式結(jié)合各種觀測進(jìn)行模擬是研究多年凍土活動層厚度在大尺度空間和時間上變化的重要工具。通過對凍土參數(shù)化方案的改進(jìn)，陸面過程模式對凍土的模擬已經(jīng)取得了很大進(jìn)展（Lawrence et al., 2008, 2012; Xie et al.,2018; Gao et al., 2019）。Guo et al.（2017）利用陸面過程模式 CLM4.5（Lawrence et al., 2012）模擬的土壤溫度插值得到活動層厚度，由此研究了多年凍土分布以及活動層厚度的變化趨勢。Gao et al.（2019）將雙向Stefan算法耦合到通用陸面過程模式CLM4.5中，發(fā)展考慮土壤凍結(jié)融化界面動態(tài)變化的陸面過程模式CAS-LSM。該模式不僅能夠模擬土壤凍結(jié)融化界面的動態(tài)變化，同時也可以模擬活動層厚度。在陸面模式中引入了求解活動層厚度的新方法。該方法有利于推進(jìn)陸面模式中多年凍土活動層模擬研究，并為活動層模擬提供額外新的有用信息。另一方面，多年凍土以及活動層厚度的模擬也受到大氣強(qiáng)迫數(shù)據(jù)集的影響（Guo et al.,2012; Lawrence et al., 2012）。Guo et al.（2012）采用動力降尺度對大氣強(qiáng)迫數(shù)據(jù)進(jìn)行了校正減小青藏高原多年凍土模擬的偏差。Lawrence et al.（2012）通過改進(jìn)大氣強(qiáng)迫數(shù)據(jù)中的氣溫和降水估計(jì)，模擬的多年凍土退化速率降低了29%。

本研究我們采用陸面過程模式CAS-LSM，利用陸面模式比較計(jì)劃 LS3MIP（Van den Hurk et al.,2016） 中 4 套 強(qiáng) 迫 數(shù) 據(jù) GSWP3（ Kim, 2017）、WFDEI（ Haddeland et al., 2011; Weedon at al.,2014）、 CRU-NCEP（ Viovy and Ciais, 2009）、Princeton（Sheffield et al., 2006），驅(qū)動 CAS-LSM進(jìn)行模擬試驗(yàn)，通過對活動層厚度及其變化趨勢在不同強(qiáng)迫下的模擬與分析，研究不同大氣強(qiáng)迫作用下陸面過程模式CAS-LSM對于活動層厚度的模擬與不確定性。

2 數(shù)據(jù)和方法

2.1 數(shù)據(jù)

本研究采用的大氣強(qiáng)迫數(shù)據(jù)是陸面模式比較計(jì)劃LS3MIP中的4組強(qiáng)迫數(shù)據(jù)，分別是CRU-NCEP（Climatic Research Unit-NCEP forcing data）（Viovy and Ciais, 2009）、 GSWP3（ Global Soil Wetness Project forcing dataset）（ Kim, 2017）、 Princeton（Princeton meteorological forcing dataset）（Sheffield et al., 2006）、 WFDEI（ water and global change forcing data methodology applied to ERA-Interim data）（Haddeland et al., 2011; Weedon at al., 2014）。強(qiáng)迫數(shù)據(jù)包括氣溫、風(fēng)速、濕度、太陽輻射和降水。CRU-NCEP是兩套數(shù)據(jù)的組合：氣候研究中心CRU TS v3.2的氣候數(shù)據(jù)和美國國家環(huán)境預(yù)報中心NCEP的再分析數(shù)據(jù)。其中降水、溫度、太陽輻射、相對濕度是基于CRU的數(shù)據(jù)，長波輻射、氣壓、風(fēng)速是由NCEP數(shù)據(jù)插值得到的，空間分辨率0.5°（緯度）×0.5°（經(jīng)度），時間分辨率 6 h，時間周期1901～2010年。GSWP3是基于20世紀(jì)的再分析數(shù)據(jù)動力降尺度而來的，利用CRU TS v3.21、全球降水氣候中心GPCC v7和地表輻射收支數(shù)據(jù)集分別對溫度、降水、長波輻射和短波輻射進(jìn)行偏差校正?？臻g分辨率0.5°（緯度）×0.5°（經(jīng)度），時間分辨率 3 h，時間周期 1901～2014年。Princeton是由 NCEP的再分析數(shù)據(jù)與CRU的數(shù)據(jù)以及基于衛(wèi)星的降水?dāng)?shù)據(jù)結(jié)合而來，其中降水?dāng)?shù)據(jù)是基于全球降水氣候項(xiàng)目（GPCP）和熱帶降雨測量項(xiàng)目（TRMM）統(tǒng)計(jì)降尺度而來?？臻g分辨率0.5°（緯度）×0.5°（經(jīng)度），時間分辨率 3 h，時間周期 1901～2012 年。WFDEI 包含兩部分?jǐn)?shù)據(jù)，WFD（WATCH force Data）（1901～1957年）基于歐洲中期天氣預(yù)報中心重新整理的ERA-40 數(shù)據(jù)得到的，WFD（WATCH force Data）（1958～2001年）基于歐洲中期天氣預(yù)報中ECMWFERA-40再分析數(shù)據(jù)得到的，其中溫度和降水均來自ERA-40，并且利用全球降水氣候項(xiàng)目（GPCP）和氣候研究中心（CRU）數(shù)據(jù)進(jìn)行偏差校正，空間分辨率0.5°（緯度）×0.5°（經(jīng)度），時間分辨率3 h，時間周期1901～2001年。WFDEI是基于ECMWFERA-Interim再分析數(shù)據(jù)得到的，其中溫度和降水均來自ERA-Interim，并且利用全球降水氣候項(xiàng)目（GPCP）和氣候研究中心（CRU）數(shù)據(jù)進(jìn)行偏差校正，空間分辨率0.5°（緯度）×0.5°（經(jīng)度），時間分辨率3 h，時間周期1979～2014年。這些強(qiáng)迫數(shù)據(jù)已被應(yīng)用與陸面模式比較LS3MIP（Van den Hurk et al., 2016）。

活動層厚度的觀測數(shù)據(jù)來自北極活動層觀測網(wǎng)（CALM）（Brown et al., 2000）（http://www.gwu.edu/～CALM/[2019-09-05])站點(diǎn)觀測數(shù)據(jù)。CALM的觀測數(shù)據(jù)時間周期是1990～2010年。本研究選取了1996～2010年期間CALM的50個觀測站的活動層厚度數(shù)據(jù)，因?yàn)榇蠖鄶?shù)觀測站在這段期間是有連續(xù)的記錄。這些站點(diǎn)主要集中在俄羅斯地區(qū)以及阿拉斯加地區(qū)。

2.2 陸面過程模型CAS-LSM的凍融過程方案

CAS-LSM是以美國國家大氣科學(xué)中心所發(fā)展的陸面過程模式 CLM4.5 為基礎(chǔ)（Oleson et al.,2013）, 耦合了自主發(fā)展的考慮地下水側(cè)向流（Zeng et al., 2016a, 2016b）、人類取用水活動（Zou et al.,2014, 2015; Zeng et al., 2017）、 凍 結(jié) 融 化 界 面（Gao et al., 2016, 2019）以 及 河流 氮 輸送 過 程（Liu et al., 2019）的方案，該模式已應(yīng)用于中國黑河流域研究（Xie et al., 2018）。模式較好地再現(xiàn)了地下水、蒸發(fā)蒸騰和凍土的分布，與地表溫度、熱通量時間變化吻合較好。本研究的工作主要涉及土壤凍結(jié)融化界面的動態(tài)變化這一模塊，該模塊描述如下：

凍結(jié)融化界面的計(jì)算方法（Jumikis, 1978; Gao et al., 2016, 2019）：

對公式（1）兩邊取平方得：

對公式（2）中的D求微分得：

其中，在陸面模式中土壤垂直分為15層，我們定義i為土壤層數(shù)，為 第i層土壤的深度，為第i層土壤的厚度，為第i層土壤的熱傳導(dǎo)率，為第i層土壤的熱傳輸阻力為凍融界面通過第i層土壤所需的凍融指數(shù)，首先計(jì)算凍融界面通過第一層土壤所需的凍融指數(shù)，即，此過程中從0增長至，將平均分成層，可近似表示為

同理可得，

那么從地表算起的凍結(jié)或融化深度可表示為

該算法與陸面過程模式CAS-LSM是實(shí)時耦合。計(jì)算時需要輸入地表溫度以及各層土壤的相關(guān)性質(zhì)，同時將得到土壤凍結(jié)融化界面信息加入到原有土壤分層中。原土壤分層信息不變，增加土壤凍結(jié)融化界面這一層。新增加這一層土壤是按照陸面模式的分層方式（Oleson et al., 2013）加入的。計(jì)算新土壤分層的土壤溫度，之后將新土壤分層的溫度返回原土壤分層。從而更加合理地模擬高寒區(qū)域凍融狀況，具體過程可參見（Oleson et al., 2013; Gao et al.,2016, 2019）。CAS-LSM不僅能夠模擬土壤凍結(jié)融化界面的動態(tài)變化，同時也可以模擬活動層厚度。該方法有利于推進(jìn)陸面模式中多年凍土活動層模擬研究，并為活動層模擬提供額外新的有用信息。

2.3 試驗(yàn)設(shè)計(jì)

本研究采用陸面過程模式CAS-LSM進(jìn)行4組數(shù)值試驗(yàn)，分別采用4組不同的大氣強(qiáng)迫數(shù)據(jù)（GSWP3、WFDEI、CRU-NCEP、Princeton）驅(qū)動CAS-LSM，進(jìn)行全球模擬，陸面模式的模擬時間積分步長為1800 s，模擬開始時首先將大氣強(qiáng)迫數(shù)據(jù)插值到陸面模式的模擬時間積分步長上，然后進(jìn)行模擬，模擬時間為1901～2010，空間分辨率為0.9°（緯度）×1.25°（經(jīng)度）。為使模式達(dá)到平衡，利用CRU-NCEP 20年（1901～1920年）的大氣強(qiáng)迫數(shù)據(jù)循環(huán)強(qiáng)迫，結(jié)果作為4組數(shù)值模擬試驗(yàn)的初始場。最后對1901～2010進(jìn)行正式模擬（月輸出結(jié)果），結(jié)果分析時段為1960～2009年。4組試驗(yàn)除了大氣強(qiáng)迫數(shù)據(jù)不同外，其余設(shè)置（如模擬區(qū)域、步長、初始場等）均相同。

2.4 分析方法

在陸面模式中多年凍土的定義：在模式土壤分層的前10層（相當(dāng)于3.8 m）土壤中，至少有1層的土壤溫度連續(xù)24個月低于0℃（Lawrence et al.,2012）。應(yīng)該指出的是這里多年凍土的定義代表的是“近地表多年凍土”。一般來說，近地表多年凍土比深層多年凍土對氣候變化更敏感（Lawrence et al., 2008）。

活動層厚度的定義是多年凍土在一年中最大融化深度（Lawrence et al., 2012）。本研究中的活動層厚度是通過CAS-LSM計(jì)算得到。CAS-LSM通過耦合Stefan算法來模擬活動層厚度。一般來說，在陸面模式中活動層厚度可以通過土壤溫度插值計(jì)算得到（Lawrence et al., 2012; Guo et al., 2017）。因此我們也給出了通過土壤溫度插值計(jì)算得到活動層厚度，并且將兩種方法所得到活動層厚度與觀測值進(jìn)行了比較。為了說明CAS-LSM模擬的活動層厚度是合理的。對于氣候態(tài)的活動層厚度驗(yàn)證，我們將每個站點(diǎn)都有記錄的觀測數(shù)據(jù)進(jìn)行氣候態(tài)平均，然后與站點(diǎn)相對應(yīng)的網(wǎng)格數(shù)據(jù)平均值進(jìn)行比較驗(yàn)證。使用平均偏差、均方根誤差和相關(guān)系數(shù)來評價模擬結(jié)果與觀測結(jié)果的一致性。采用最小二乘擬合法計(jì)算趨勢作為線性回歸的斜率，采用Mann- Kendall檢驗(yàn)評價線性趨勢的統(tǒng)計(jì)意義。

圖1 1996～2010 年（a）CRU-NCEP、（b）GSWP3、（c）Princeton、（d）WFDEI 4 組大氣強(qiáng)迫數(shù)據(jù)模擬的活動層厚度（藍(lán)點(diǎn)）、土壤溫度插值的活動層厚度（紅點(diǎn)）與CALM觀測數(shù)據(jù)驗(yàn)證Fig. 1 Simulated active layer thickness validation (blue point) using four atmosphere forcing data (a) CRU-NCEP (Climatic Research Unit-NCEP forcing data), (b) GSWP3 (Global Soil Wetness Project forcing dataset), (c) Princeton (Princeton meteorological forcing dataset), and (d) WFDEI(water and global change forcing data methodology applied to ERA-Interim data) and the soil temperatures interpolationed active layer thickness validation (red point) based on observations from CALM (Circumpolar Active Layer Monitoring) sites during 1996–2010

3 模擬多年凍土區(qū)活動層厚度的變化

3.1 站點(diǎn)驗(yàn)證

首先將不同強(qiáng)迫下模擬所得到的活動層厚度數(shù)據(jù)與觀測數(shù)據(jù)進(jìn)行驗(yàn)證，圖1給出了4種強(qiáng)迫下（CRU-NCEP、GSWP3、Princeton、WFDEI）模擬的1996～2010年平均活動層厚度與CALM觀測數(shù)據(jù)驗(yàn)證的散點(diǎn)圖。藍(lán)色的點(diǎn)為CAS-LSM模擬得到的活動層厚度，紅色的點(diǎn)為土壤溫度插值得到的活動層厚度數(shù)據(jù)。由圖1可知不同強(qiáng)迫下CASLSM模擬得到的活動層厚度和觀測值比較接近，土壤溫度插值得到的活動層厚度和觀測值也比較接近。同時與觀測相比都表現(xiàn)出偏大的趨勢?；顒訉雍穸扔^測數(shù)據(jù)的平均值為0.86 m，CAS-LSM模擬得到的活動層厚度分別為 1.47 m、1.28 m、1.48 m、1.26 m。土壤溫度插值得到的活動層厚度分別為1.28 m、1.19 m、1.4 m、1.18 m。CAS-LSM 模擬得到的活動層厚度與觀測的平均偏差分別為0.61 m、0.42 m、0.62 m、0.40 m，均方根誤差分別為 1.04 m、0.85 m、0.92 m、0.78 m，相關(guān)系數(shù)分別為 0.44、0.53、0.61、0.64。土壤溫度插值得到的活動層厚度與觀測值的平均偏差分別為0.42 m、0.33 m、0.54 m、0.32 m，均方根誤差分別為0.85 m、0.86 m、0.81 m、0.79 m，相關(guān)系數(shù)分別為 0.47、0.43、0.62、0.55。從平均偏差、均方根誤差與相關(guān)系數(shù)上可以看出，這兩種該方法模擬得到的活動層厚度差別不大，從而表明CAL-LSM模擬的活動層厚度合理。從平均偏差、均方根誤差與相關(guān)系數(shù)上可以看出，強(qiáng)迫數(shù)據(jù)WFDEI模擬的活動層厚度與觀測數(shù)據(jù)最為接近，整體結(jié)果優(yōu)于其他3種強(qiáng)迫數(shù)據(jù)的模擬。

圖2 （a）CRU-NCEP、（b）GSWP3、（c）Princeton、（d）WFDEI 4 組大氣強(qiáng)迫數(shù)據(jù)模擬的 1960～2009 年活動層厚度的空間分布Fig. 2 Simulated mean active layer thickness using four atmosphere forcing data (a) CRU-NCEP, (b) GSWP3, (c) Princeton, and (d) WFDEI during 1960–2009

3.2 多年凍土活動層厚度空間分布及變化趨勢

我們主要分析多年凍土活動層厚度50年（1960～2009年）平均的空間分布及變化趨勢。從圖2可以看出，4組試驗(yàn)?zāi)M的活動層厚度在氣候態(tài)（1960～2009年）平均上表現(xiàn)出空間相似性。除青藏高原外，模擬的活動層厚度隨緯度的增加而減小，活動層厚度在靠近多年凍土與季節(jié)性凍土分界線的區(qū)域深度較深。表1中區(qū)域平均的活動層厚可以看出CRU-NCEP和Princeton模擬的活動層厚度較深，WFDEI和GSWP3模擬的活動層厚度較淺。圖3中活動層厚度的變化趨勢基本為正值，表現(xiàn)逐年增加的趨勢，區(qū)域平均的變化趨勢分別為0.3 cm/a、0.295 cm/a、0.27 cm/a 和 0.5 cm/a，4 組模擬的平均趨勢為0.34 cm/a。青藏高原地區(qū)（28°N～38°N，73°E～104°E）、北美地區(qū)（56°N～70°N，70°W～165°W）、歐亞大陸北部地區(qū)（54°N～72°N，72°E～178°E）的區(qū)域平均活動層厚度以及變化趨勢見表1。4組模擬在這3個地區(qū)的區(qū)域平均厚度及變化趨勢分別為2.53 m、1.16 m、1.38 m 和0.56 cm/a、0.35 cm/a、0.37 cm/a。圖 4 模擬的是青藏高原地區(qū) （ 28°N～ 38°N， 3°E～ 104°E）、 北 美 地 區(qū)（56°N～70°N，70°W～165°W）、歐亞大陸北部地區(qū)（54°N～72°N，72°E～178°E）1960～2009 年活動層厚度年際變化。從圖4中可以看出，1960～2009年活動層厚度總體上表現(xiàn)出增加的趨勢。

表1 1960～2009年利用4組大氣強(qiáng)迫數(shù)據(jù)模擬的活動層厚度區(qū)域平均值及其變化趨勢Table 1 Simulated regional average active layer thickness and its trends from 1960 to 2009 using four atmosphere forcing data

圖3 （a）CRU-NCEP、（b）GSWP3、（c）Princeton、（d）WFDEI 4 組大氣強(qiáng)迫數(shù)據(jù)模擬的 1960～2009 年活動層厚度變化趨勢Fig. 3 Simulated active layer thickness trends (cm/a) using four atmosphere forcing data (a) CRU-NCEP, (b) GSWP3, (c) Princeton, and (d) WFDEI during 1960–2009

圖4 （a）CRU-NCEP、（b）GSWP3、（c）Princeton、（d）WFDEI 4 組大氣強(qiáng)迫數(shù)據(jù)模擬的 1960～2009 年區(qū)域平均活動層厚度年時間序列。黑線：青藏高原地區(qū)（28°N～38°N，73°E～104°E），紅線：北美地區(qū)（56°N～70°N，70°W～165°W），藍(lán)線：歐亞大陸北部地區(qū)（54°N～72°N，72°E～178°E）Fig. 4 Simulated active layer thickness annual time series using four atmosphere forcing data (a) CRU-NCEP, (b) GSWP3, (c) Princeton, and (d)WFDEI during 1960–2009. Black line: Qinghai–Tibet plateau region (28°N–38°N, 73°E–104°E), red line: North America (56°N–70°N, 70°E–165°E),blue line: north of Eurasia (54°N–72°N, 72°E–178°E)

對于不同強(qiáng)迫造成的活動層厚度模擬差異，我們主要從強(qiáng)迫數(shù)據(jù)本身分析。首先我們采用Stefan方法模擬活動層厚度，由于地表溫度是形成多年凍土的關(guān)鍵因素，Stefan利用地表溫度來計(jì)算單層土柱凍融深度，是一種半經(jīng)驗(yàn)方法。不同強(qiáng)迫模擬地表溫度的不同會造成活動層厚度模擬的差異，從圖5中我們可以看出，地表溫度變化趨勢基本為正值，也就是說地表溫度整體呈現(xiàn)出增加的趨勢，變化趨勢分別為 0.036℃/a、0.036℃/a、0.035℃/a和0.061℃/a，4組模擬的平均趨勢為0.042℃/a。由于地表溫度增加，進(jìn)而會導(dǎo)致活動層厚度的增加，地表溫度的增加趨勢是與活動層厚度的增加是一致的。圖3和圖5中青藏高原的部分區(qū)域出現(xiàn)了活動層厚度變薄、地表溫度降低的趨勢。然而強(qiáng)迫數(shù)據(jù)中的2 m氣溫所在的這些區(qū)域都是增加的。這可能由于陸面模式在該地區(qū)的地表參數(shù)、物理參數(shù)化方案等造成的，具體還需要針對這一地區(qū)結(jié)合觀測數(shù)據(jù)進(jìn)行深入研究?？傮w來看，除青藏高原的部分區(qū)域外，2 m氣溫以及地表溫度的變化趨勢與活動層厚度的變化趨勢是一致的。4種不同強(qiáng)迫模擬的地表溫度差異主要是由于強(qiáng)迫數(shù)據(jù)中2 m氣溫的差異造成的，因此我們分析了2 m的氣溫變化。從圖6可以看出，2 m氣溫變化趨勢基本為正值，表現(xiàn)出增加的趨勢，也就是說2 m氣溫整體呈現(xiàn)出增加的趨勢，變化趨勢分別為0.036℃/a、0.035℃/a、0.035℃/a和0.054℃/a，4組模擬的平均趨勢為0.040℃/a。除青藏高原的部分區(qū)域外，2 m氣溫以及地表溫度及活動層厚度的變化趨勢是一致的。

圖5 （a）CRU-NCEP、（b）GSWP3、（c）Princeton、（d）WFDEI 4 組大氣強(qiáng)迫數(shù)據(jù)模擬的 1960～2009 年地表溫度變化趨勢Fig. 5 Simulated ground temperature trends (℃/a) using four atmosphere forcing data (a) CRU-NCEP, (b) GSWP3, (c) Princeton, and (d) WFDEI during 1960–2009

圖7模擬的是1960～2009年多年凍土區(qū)2 m氣溫與活動層厚度年際變化，從圖中可以看出，2 m氣溫的變化趨勢與活動層厚度變化趨勢顯著相關(guān)，相關(guān)系數(shù)分別為0.82、0.83、0.73、0.87。隨著氣溫的升高，活動層厚度表現(xiàn)出增加的趨勢。表明活動層厚度的模擬與2 m氣溫有很強(qiáng)的相關(guān)性。圖8模擬的是1960～2009年多年凍土區(qū)2 m氣溫變化趨勢與活動層厚度變化趨勢相關(guān)系數(shù)的空間分布。從相關(guān)系數(shù)的空間分布上可以看出，不同強(qiáng)迫下2 m氣溫變化趨勢與活動層厚度變化趨勢整體上顯著相關(guān)，但不同的區(qū)域會有差異。俄羅斯東部和西部、加拿大東部等地區(qū)均呈現(xiàn)出較高的相關(guān)性，說明這些地區(qū)活動層厚度模擬受氣溫影響較大。俄羅斯中部以及阿拉斯加地區(qū)相關(guān)性較低，說明這些地區(qū)除了氣溫影響外，同時海拔、活動層巖性、含水特征以及影響地面溫度變化進(jìn)程的地形特征和下墊面性質(zhì)等也對活動層厚度模擬影響較大。在陸面模式中可表示為土壤性質(zhì)、模式參數(shù)、物理參數(shù)化方案等對活動層厚度模擬的影響，具體區(qū)域還需要結(jié)合觀測數(shù)據(jù)進(jìn)一步進(jìn)行研究?？傮w來看，強(qiáng)迫數(shù)據(jù)中2 m氣溫對于活動層厚度模擬的影響還是占主要因素。

圖6 （a）CRU-NCEP、（b）GSWP3、（c）Princeton、（d）WFDEI 4 組大氣強(qiáng)迫數(shù)據(jù)模擬的 1960～2009 年 2 m 氣溫變化趨勢Fig. 6 Simulated 2-m air temperature trends (℃/a) using four atmosphere forcing data (a) CRU-NCEP, (b) GSWP3, (c) Princeton, and (d) WFDEI during 1960–2009

4 討論與結(jié)論

本研究主要是在不同大氣強(qiáng)迫作用下，對陸面模式CAS-LSM模擬的多年凍土活動層厚度進(jìn)行評估。分別采用4種不同的大氣強(qiáng)迫數(shù)據(jù)（SWP3、CRU-NCEP、Princeton、WFDEI）驅(qū)動CAS-LSM，并分析不同強(qiáng)迫下的活動層厚度空間分布及變化趨勢。通過與CALM活動層的觀測數(shù)據(jù)比較，CASLSM模擬得到的活動層厚度與觀測值接近較為一致，但表現(xiàn)出偏大的趨勢。主要是由于Stenfan公式在計(jì)算活動層厚度時，假設(shè)吸收的熱量全部用于多年凍土的融化，忽略了土壤本身溫度變化的感熱影響（Jumikis, 1978; Woo et al., 2004），因此造成了模擬結(jié)果偏大。由于地表溫度是形成多年凍土的關(guān)鍵因素，Stefan利用地表溫度來計(jì)算單層土柱凍融深度。不同強(qiáng)迫模擬地表溫度的不同會造成活動層厚度模擬的差異，而地表溫度模擬的差異主要是由于強(qiáng)迫數(shù)據(jù)中的2 m氣溫的差異造成的。通過對圖3、5、6的分析表明：總體來看，除青藏高原的部分區(qū)域外，2 m氣溫以及地表溫度的變化趨勢與活動層厚度的變化趨勢是一致的。青藏高原的部分區(qū)域出現(xiàn)了活動層厚度變薄、地表溫度降低的趨勢而2 m氣溫是增加的。這可能由于陸面模式在該地區(qū)的地表參數(shù)、物理參數(shù)化方案等造成的，具體還需要針對這一地區(qū)結(jié)合觀測數(shù)據(jù)進(jìn)行深入研究。

對于不同強(qiáng)迫造成的活動層厚度相對變化的不確定性，我們首先計(jì)算活動層厚度的相對變化，定義為t年活動層厚度，年到年的活動層厚度變化趨勢為，那么年相對于年活動層厚度的相對變化：

那么2009年相對于1960年的活動層厚度相對變化：

不同強(qiáng)迫造成的活動層厚度相對變化的不確定性范圍是基于不同強(qiáng)迫數(shù)據(jù)模擬相對變化的標(biāo)準(zhǔn)差求得的。從活動層厚度區(qū)域平均上，4種強(qiáng)迫模擬的活動層厚度相對變化分別為11.2%、13.5%、10.3%、23.5%，不確定性范圍為14.6%±6.1%。圖9模擬的是多年凍土區(qū)活動層厚度在不同強(qiáng)迫下相對變化標(biāo)準(zhǔn)差的空間分布，標(biāo)準(zhǔn)差越大說明不同強(qiáng)迫對于本地區(qū)的模擬差異越大，不確定性越大。這些地區(qū)的強(qiáng)迫數(shù)據(jù)對于活動層厚度相對變化的模擬影響越大。從圖9中可以看出在青藏高原、俄羅斯的西部地區(qū)、阿拉斯加地區(qū)等地區(qū)差異較大。因此這些地區(qū)的不確定性范圍也較大，在這些地區(qū)不同強(qiáng)迫對于活動層厚度相對變化模擬影響也較大。

圖7 （a）CRU-NCEP、（b）GSWP3、（c）Princeton、（d）WFDEI 4 組大氣強(qiáng)迫數(shù)據(jù)模擬的 1960～2009 年多年凍土區(qū) 2 m 氣溫（紅線）與活動層厚度（藍(lán)線）年際變化趨勢Fig. 7 Simulated interannual trends of 2-m air temperature (red line) and active layer thickness (blue line) using four atmosphere forcing data (a)CRU-NCEP, (b) GSWP3, (c) Princeton, and (d) WFDEI during 1960–2009

本研究的主要結(jié)論如下：

（1）通過與多年凍土區(qū)活動層厚度的觀測數(shù)據(jù)以及土壤溫度插值得到的結(jié)果比較，CAS-LSM模擬的結(jié)果與觀測值比較接近?；赪FDEI模擬的活動層厚度要優(yōu)于其他3種強(qiáng)迫數(shù)據(jù)模擬的。

（2）對于不同強(qiáng)迫造成的活動層厚度模擬的差異，主要原因是由于強(qiáng)迫數(shù)據(jù)中2 m氣溫的差異造成的。2 m氣溫的差異造成模擬的地表溫度不同，從而導(dǎo)致活動層厚度模擬的差異。通過分析表明：除青藏高原的部分區(qū)域外，活動層厚度的變化趨勢與地表溫度以及2 m氣溫變化趨勢是一致的。總體來看，強(qiáng)迫數(shù)據(jù)中2 m氣溫對于活動層厚度模擬的影響還是占主要因素。不同大氣強(qiáng)迫模擬的活動層厚度區(qū)域平均和變化趨勢范圍分別為1.1～1.25 m和 0.27～0.51 cm/a。

圖8 （ a） CRU-NCEP、（ b） GSWP3、（ c） Princeton、（ d）WFDEI 4組大氣強(qiáng)迫數(shù)據(jù)模擬的1960～2009年多年凍土區(qū)2 m氣溫與活動層厚度相關(guān)系數(shù)的空間分布Fig. 8 Spatial distribution of correlation coefficients between simulated 2-m air temperature and simulated active layer thickness using four atmosphere forcing data (a) CRU-NCEP, (b) GSWP3, (c)Princeton, and (d) WFDEI during 1960–2009

圖9 多年凍土區(qū)活動層厚度在不同強(qiáng)迫下相對變化標(biāo)準(zhǔn)差的空間分布Fig. 9 Spatial distribution of relative variation standard deviation of active layer thickness in permafrost under different forces data

（3）不同強(qiáng)迫數(shù)據(jù)模擬的活動層厚度變化趨勢與2 m氣溫變化趨勢整體上相關(guān)，但不同的區(qū)域會有差異。同時不同強(qiáng)迫數(shù)據(jù)對于模擬活動層厚度相對變化的影響較大。模擬的活動層厚度相對變化分別為 11.2%、13.5%、10.3%和23.5%，不確定性范圍為14.6%±6.1%。在青藏高原、俄羅斯的西部地區(qū)、阿拉斯加地區(qū)等地區(qū)的不確定性范圍較大，在這些地區(qū)不同強(qiáng)迫對于活動層厚度相對變化模擬影較大。