本研究选用了18个CMIP6 模式的历史气候试验和情景模式比较计划模拟结果（https： //esgf-node.llnl.gov/projects/cmip6/）， 如表1所示。历史气候模拟试验是基于实际观测的外强迫驱动下对1850 -2014年的历史气候进行的演变。该试验用于评估模式对气候变化的模拟能力， 还可以用来分析气候模式与观测记录的一致性等（周天军等， 2019）。情景模式比较计划模拟结果将选用更新后的高排放情景SSP5-8.5， 新辐射强迫情景SSP3-7.0， 更新后的中等排放情景SSP2-4.5， 更新后的低排放情景SSP1-2.6， 到2100年人为辐射强迫分布分别达到8.5， 7.0， 4.5和2.6 W·m^-2。其中SSP 描述了在没有气候变化或者气候政策影响下， 未来社会的可能发展情形， SSP1、 SSP2、 SSP3和 SSP5 分别代表了可持续发展、 中度发展、 局部发展和常规发展 4 种路径（张丽霞等， 2019）。研究选用的模式变量包括近地面气温（tas）、 积雪深度（snd）、 雪水当量（snw）。青藏高原区域范围选取为25°N -40°N， 75°E -105°E， 包括其主体与周边地区（李振朝， 2013）。进行分析研究所选用的当前和未来时段分别为1985 -2014年、 2015 -2100年。

为了对CMIP6模式的模拟结果进行评估， 本研究将使用中国气象数据网（http： //data.cma.cn/）所提供的中国地面气温月值0.5°×0.5°格点数据集（V2.0）， 中国地面降水月值0.5°×0.5°格点数据集（V2.0）作为观测资料， 主要用于评估青藏高原的气温、 积雪深度等与冻土变化密切相关的变量。上述格点数据集利用最新整编的中国地面高密度台站数据和相关的空间插值方法得到， 有相对较高的空间分辨率， 并进行了相对严格的验证， 具有较高的可信度（赵煜飞等， 2014； 张强等， 2009）。

数据分析过程中， 由于各模式间、 模式与观测间数据分辨率差异较大， 为了便于定性定量的分析比较， 统一采用双线性插值的方法将模式结果插值到0.5°×0.5°的格点上。本文将对高原冻土相关的气候变量的空间分布、 趋势、 误差等以及对冻土面积的模拟能力进行评估。由于多模式的集合平均通常将会产生更好的模拟效果（Sillmann et al， 2013）， 在此基础上根据评估结果， 将CMIP6模式结果采用算术平均的方式进行多模式集合平均（The method of the multi-model ensemble， MME）， 对未来冻土面积变化进行预估。

冻土面积的估算使用地面冻结数模型（SFI）， 由于青藏高原观测资料的匮乏， 该模型可以利用较少的气候变量对地表层多年冻土进行诊断。该模型由Nelson and Outcalt（1987）在空气冻结数的基础上修正提出， 计算公式如下：

式中： D D F^{+}是考虑了积雪绝对修正性后冻结天数的总和， 即气温在零度以下的天数和； D D T是根据正弦气候计算出的融化天数的总和， 即气温在零度以上的天数和。该模型指数用于确定多年冻土与季节冻土、 多年冻土区中的连续与不连续多年冻土区间的平面分布界线。指数范围从0.0~1.0， 0.5以下表示无多年冻土， 0.5~0.6表示存在零星的多年冻土， 0.60~0.67表示存在广泛但不连续的多年冻土， 0.67以上表示存在连续的多年冻土（南卓铜， 2012； Salter and Lawrence， 2013； 常燕等， 2016）。Burn and Nelson（2006）指出模型只能捕捉不连续和连续性多年冻土。

该模型在使用时需要输入年最高气温、 最低气温以及平均冬季积雪深度， 用以产生正弦气候态。这里， 将前一年的10月到次年的9月作为一年时间。在输入年最高、 最低气温时采用研究时段的前20年平均（例如1985年的数据时使用1966 -1985年的平均数据）。平均冬季积雪深度则使用气温低于零度月份的积雪深度 Z_i与气温 T_i加权平均。计算公式如下：

另外雪水当量（snw）将与积雪深度（snd）用于计算积雪密度， 积雪密度用于估计积雪的导热率， 而平均冬季积雪密度与积雪深度的加权方式相同。对于未提供积雪深度或雪水当量数据的模式， 将采用250 kg·m^-3当作积雪密度用于计算（Salter and Lawrence， 2013）。

观测资料的积雪深度则需要根据月降水量数据计算出， 计算公式如下（Anisimov and Nelson， 1997； Nelson and Outcalt， 1987）：

式中： P_i为当气温低于零度时， k（ i=1， 2， 3， …， 12）月份时的降水量； {\rho }_i代表积雪密度取250 kg·m^-3； \phi表示所考虑位置的纬度。

在SFI方法中， 近地面气温、 积雪深度这两个气候变量与冻土的变化密切相关。温度作为影响冻土形成和变化最主要的因子， 温度变化将影响着青藏高原多年冻土的发育与分布， 随着气候持续变暖， 青藏高原多年冻土呈区域性退化的状态（石亚亚等， 2017）； 积雪对冻土的存在与消退也有着很大的关系， 雪盖对地面既有保温也有增湿的作用， 积雪消融吸收热量又可降低地面温度， 对土壤的水热属性有着极大的影响（王婷等， 2019； 高荣等， 2003）。为了有效地计算出青藏高原地表层多年冻土的面积， 将对各模式气温和积雪深度的模拟结果与观测资料进行分析比较， 验证模式模拟高原气温和积雪深度的能力。通常多年冻土的界限与年均温度的0 ℃线较为吻合， 因此将CMIP6各模式在高原地区的模拟结果与观测资料进行比较， 从图1可以看出， 各模式模拟的年均气温0 ℃线与观测值除了在高原东部以及南部差异较为显著外， 其余地方在高原的位置大致相同。其在高原东部与南部差异大， 可能是与此地区海拔以及地形变化较大的原因， 而各模式不能很好地捕捉地形以及海拔变化对温度变化带来的影响。

图2显示了1985 -2014年间18个模式、 多模式集合平均和观测资料在青藏高原上的年均温度的空间分布。从图2中能够看出， 各模式在青藏高原区域模拟的年均温度的空间分布与观测资料较为一致， 气温由东向西逐渐降低。在青藏高原东部与东南部都存在着气温偏高的区域， 中西部地区存在气温偏低的大值区， 模拟的最低温度从-30~-12 ℃不等， 大部分模式模拟的温度空间分布与观测资料相似， 但在数值上有较大的差异， 特别是冷中心气温差值可达10 ℃以上。观测资料和模式模拟结果的空间相关系数为0.555~0.727（表2）， 其中空间相关系数最好模式的是EC-Earth3。虽然空间相关系数可以表示模式对主要气温中心位置的描述能力， 但不能表明数值的差异大小。而均方根误差（RMSE）则显示了模拟温度与观测资料之间的平均误差， 这18个模式的RMSE从0.429~9.460 ℃不等， 均方根误差最小的模式是CESM-WACCM， 这表明该模式所模拟的结果与观测资料数值最为接近。同时表2也给出了各模式模拟的年平均温度与观测气温的时间相关系数， 利用T检验对各模式相关性进行检验能够看出18个模式中有10个模式通过了0.01的水平显著性检验， 其中CESM-WACCM和CanESM5模式模拟结果与观测资料的相关性最好， 分别是0.787和0.744。多模式的集合平均（MME）与观测资料的年均气温时间和空间相关系数较多个模式有了很大的提升， 特别是时间相关系数MME优于17个模式， 可以看出集合平均后得到的模拟效果更好。

图3给出了1985 -2014年各模式和观测资料在青藏高原上平均冬季积雪深度的空间分布。从图3中能够看到， 各模式模拟的积雪范围与观测资料所得出来的范围除青藏高原东南部存在差异以外其他部分是较为一致的， 但各模式模拟的积雪深度值基本都大于观测资料， 这可能是由于观测资料的积雪深度是利用降水转化而造成的。多个模式在青藏高原南部或东南部有一个明显的积雪区， 同时CAS-ESM2、 CESM2-WACCM、 CMCC-CM2-SR5以及CMCC-ESM2所模拟的主要积雪区域范围较大包含了三江源以及祁连山脉地区。大多数模式模拟的冬季平均积雪深度由东向西明显减小， 但其中CanESM5和KIOST-ESM以及IPSL-CM6A-LR分别在高原西部和高原中部存在大面积的积雪区域， 积雪深度呈现由东向西减小-增大的趋势。多模式集合平均后的冬季平均积雪深度在高原南部山脉以及高原中东部的“三江源”地区呈现出主要的积雪区域， 且在高原西部有一个积雪深度变大的区域， 总体上积雪深度呈现出由南向北逐渐减少， 由东向西逐渐减少-增大的趋势， 这样的分布特征与韦志刚和吕世华（1995）以及罗江鑫等（2020）的研究成果较为一致， 这说明多模式集合平均后能产生更好的模拟结果。

在利用SFI方法计算地表层多年冻土面积的过程中， 关于气温需要输入年均最高、 最低气温用于计算冻结和融化的持续时长； 关于积雪需要输入冬季平均积雪深度和密度， 图4和图5分别展示了各模式和多模式集合平均在1985 -2014年的年均气温、 年均最高和最低气温、 平均冬季积雪深度相对于观测资料的距平。从图4中可以看出， 对于年均气温， 18个模式模拟的结果中只有CESM2-WACCM相对于观测资料出现暖偏差； 在年均最高气温上， ACCESS-ESM1-5、 CESM2-WACCM、 CIESM、 CMCC-CM-SR5、 CMCC-ESM2、 MIROC6以及MPI-ESM1-2-LR出现暖偏差以外其余模式均出现冷偏差； 而在年均最低气温上， 只有ACCESS-CM2和ACCESS-ESM1-5这两个模式出现较为明显的暖偏差， 其余17个模式都呈现冷偏差， 其中CanESM5和IPSL-CM6A-LR出现了高于4 ℃甚至是10 ℃以上的冷偏差。各模式之间年均气温的变化幅度超过9.5 ℃以上， 年均最高气温变化幅度最大超过8.4 ℃以上， 相较于年均气温和年均最高， 年均最低气温变化幅度更大超过11.3 ℃以上， 因此从各模式温度距平的变化幅度可以看出， 各模式之间在青藏高原冬季气温也就是年均最低气温上的模拟出现较大的差异。相对于观测资料， 绝大多数模式在一定程度低估了高原气温， 特别是冬季气温。从图5能够看出， CMIP6各模式对于青藏高原积雪深度的模拟距平差异也是较大， 18个模式中有2个模式模拟的冬季平均积雪深度的距平的绝对值>1 m； 最大距平值为1.66 m（IPSL-CM6A-LR）， 最小距平值为0.17 m（ACCESS-CM2）。各模式对于积雪深度的模拟与观测资料相比， 均高估了积雪深度， 这可能是由于大多数模式在气温的模拟上都低估了高原的气温所造成的， 也可能是观测资料的积雪深度是由降水资料转换而带来的偏差。而多模式集合平均（MME）的年均气温、 年均最高和最低气温相对于观测资料的平均值均偏低， 分别为-2.44 ℃、 -0.77 ℃、 -2.94 ℃（图4）， 而平均冬季积雪深度的距平值为0.56 m（图5）， 这能够看出多模式集合平均后的结果是优于大多数模式的， 但此结果相对于观测资料仍然存在一定的偏差， 这也可能导致后续的研究与实际情况出现偏差。

IPCC自1990年发布第一次评估报告以来， 全球的平均气温升幅一直攀升， 相关学者观测并预估到更高的气候变暖将发生在高纬度、 高海拔地区。随着全球气候变暖的发生， 青藏高原作为气候的“敏感器”， 高原上的气温持续升高（王绍令等， 1996； Guo and Wang， 2012； 杨建平等， 2013）。从图6中能够看出， 观测资料关于青藏高原的近地面年均气温从1986 -2014年气温升高约0.49 ℃·（10a）^-1， 而多模式集合平均则是0.27 ℃·（10a）^-1。在1985 -2014年间， 在1985 -1998年间观测资料距平明显低于多模式集合平均距平， 而在1998年后观测资料的距平变化明显高于多模式集合平均的变化。由于CIESM和KIOST-ESM模式缺少SSP3-7.0情景下的模拟数据， 因此SSP3-7.0情景下只有16个模式， 这也使得在历史模拟阶段SSP3-7.0的模拟结果与其余情景下的数据出现了0.037 ℃左右的偏差。SSP1-2.6、 SSP2-4.5、 SSP3-7.0和SSP5-8.5四种排放情景下， 青藏高原年平均气温整体基本呈线性上升的趋势， 其中SSP1-2.6在2077年温度到达峰值后逐渐减小， SSP2-4.5在2077年后升温幅度逐渐放缓， SSP3-7.0和SSP5-8.5在2045年后升温幅度加剧。2015-2100年间各情景下气温增长率分别为0.19， 0.36， 0.53和0.68 ℃·（10a）^-1； 在2015 -2040年间， 各情景下的气温变化幅度较为一致， 在此之后各情景下的气温变化差异逐渐增大， 2041 -2100年间气温升温率分别为0.09， 0.33， 0.59和0.80 ℃·（10a）^-1。到2100年， SSP1-2.6、 SSP2-4.5、 SSP3-7.0和SSP5-8.5的年均气温分别增长了2.20， 3.51， 5.01和6.33 ℃。

随着全球气候变暖， 青藏高原气温持续升高， 这对位于青藏高原的冻土来说有着巨大影响， 持续升高的气温必将引起冻土面积的减少。相关研究均表明气候变暖导致高原上的多年冻土呈现区域性退化， 表现为冻土的下界的升高、 冻土年平均温度升高、 活动层厚度增大、 冻结深度减小等现象（王绍令等， 1996； 高荣等， 2003； Cheng and Wu， 2007； 蒋元春等， 2020）。随着冻土的区域性退化， 多年冻土的分布和面积也随之变化。研究采用多种方法对青藏高原多年冻土面积进行了估算， 由于所采用研究数据与分析方法不同多年冻土面积差异也较大， 分别为150×10⁴ km²（周幼吾和郭东信， 1982； 金会军等， 2000）、 140.1×10⁴ km²（Li and Cheng， 1999）、 129.4×10⁴ km²（李新和程国栋， 1999）、 119.5×10⁴ km²（南卓铜等， 2002）、 106×10⁴ km²（Zou et al， 2017）。

为保持多模式模拟结果与青藏高原冻土分布的一致性， 研究中多年冻土的冻结指数取为0.6。利用SFI模型， 使用观测资料计算青藏高原1985 -2014年多年冻土的分布结果［图7（b）］与青藏高原冻土分布［图7（a）］较为一致。使用CMIP6模式资料计算， 多模式集合平均模拟的1985 -2014年多年冻土主体分布与高原多年冻土图较为一致， 但明显能够看出其并不能完全表现出多年冻土的分布， 在高原的藏南冈底斯山脉以北、 念青唐古拉山脉附近以及喜马拉雅山脉地区的模拟结果范围偏大， 不能很好地捕捉零星、 岛状的多年冻土分布情况。造成此类差异的原因可能是因为选用历史模拟时段和实际冻土图的时段存在差异的关系； 也可能是因为CMIP6模式的分辨率相对较粗、 对青藏高原的气温和积雪的模拟有着较大的冷偏差、 对青藏高原地形地势变化捕捉能力以及各模式本身物理机制和参数化存在不足而导致的。当SFI指数取为0.6时， 多年冻土面积约为134.52×10⁴ km²（包含冰川和湖泊面积）， 该数值在之前研究中的106×10⁴~150×10⁴ km²的范围之内， 说明本数值能够用于计算高原多年冻土的面积。

对未来高原多年冻土空间分布的变化情况依旧采用多模式的集合平均， 利用CMIP6情景模式比较计划所提供的数据进行预估。我们将21世纪分为近期（2021 -2040年）、 中期（2041 -2060年）、 末期（2081 -2100年）三个阶段， 并从三个阶段来预估高原多年冻土的变化（图8）。从图8中能够明显地看出， 不同的情景下由于外强迫的不同， 气温升高的速率也不一致， 多年冻土的面积在各情景下均呈现减小的趋势， 且随着外强迫的增大， 减小的趋势渐为明显； 多年冻土从其边界开始出现明显退化， 其退化区域主要出现在多年冻土区的东部与南部边缘地区， 且退化趋势呈现东南至西北一线， 这种变化趋势与气温从东至西升高一致。这与Li and Cheng（1999）的研究， 高原的南部与东部的多年冻土可能会经历大范围的退化甚至是可能消失的结果一致。至21世纪中期， 各情境下多年冻土的面积分别减小10.3%， 14.5%， 13.4%和18.5%； 至2100年时， 四种情景下高原多年冻土主要分布在高原的西北区域， 多年冻土面积最小的是在SSP5-8.5情景下仅存35.39×10⁴ km²。但值得注意的是， 21世纪末期时多区域多年冻土已经完全退化， 但本研究数据网格较大并不能完全捕捉到零星、 面积较小的多年冻土区域， 且该数据代表的是地表层多年冻土面积消失， 必然存在当表面的多年冻土退化时， 地表深层仍然存在着多年冻土的情况， 同时也意味着维持多年冻土的这些气候条件不再可持续（Slater and Lawrence， 2013）。

通过对不同情景下CMIP6多模式集合平均后不同阶段的面积计算发现， 高原的多年冻土面积呈现减少的趋势。21世纪近期 SSP1-2.6、 SSP2-4.5、 SSP3-7.0和SSP5-8.5四种不同情景下模拟的多年冻土面积分别为136.81×10⁴ km²、 136.63×10⁴ km²、 135.30×10⁴ km²和132.49×10⁴ km²。21世纪末期， 随着人为强迫的差距加大， 模式模拟的冻土面积减少趋势逐渐出现较大差距， 各情景下多年冻土面积分别仅存93.19×10⁴ km²、 70.59×10⁴ km²、 46.33×10⁴ km²和35.3910⁴ km²。SSP1-2.6情景下， 多模式集合平均后的气温刚好达到2.2 ℃， 到2100年时多年冻土面积较历史模拟减少了41.33×10⁴ km²约占30.7%； SSP2-4.5情景下， 多年冻土面积较历史模拟减少了63.93×10⁴ km²约47.52%； SSP3-7.0情景下， 至2100年时多年冻土面积减少约88.20×10⁴ km²占65.6%； SSP5-8.5情景下， 至21世纪末期时多年冻土基本退化， 面积减少约99.14×10⁴ km²占73.7%。