CMIP6 BCC等模式对青藏高原土壤冻融模拟性能分析

周鑫原 , 吕世华 , 罗江鑫

（成都信息工程大学大气科学学院, 成都 610225）

引言

地球气候系统中的冰冻圈是仅次于大气圈的第二大圈层[1]。冰冻圈具有相变潜热巨大、反照率高、冷储高的特点[2-3]。冻土是重要的陆面强迫因子，土壤冻融过程可以直接参与能量、水分循环，影响局地天气和气候[4-9]。中低纬度带上，冻土面积最广（约占中国多年冻土面积的70%）、厚度最大的地区主要分布在具有特殊地理位置和高海拔的青藏高原[10-11]。相对而言，青藏高原冬季温度较高，冻土厚度较小且稳定性差，对气候变化的响应更为敏感[12]。因此，青藏高原土壤冻融过程与气候变化息息相关，一直以来受到国内外气象学者的高度关注。

在IPCC第五次评估报告中，全球变暖再次得到证实[13]。气象观测站资料以及钻孔监测资料分别显示：1961～2007年，高原中部和东部70多个海拔高于2000 m站点的平均增温率为0.28°C/10 a[13]；过去几十年，高原北部和南部多年冻土下界也有明显的升高趋势，分别升高了25 m、50～80 m[14]。一直以来，高原地区由于海拔过高、地形复杂、气候环境恶劣，还存在大量的无人区，观测资料匮乏，严重制约了对该地区土壤冻融过程研究的深入开展。随着气象观测能力和计算机技术的飞速发展，数值模拟已被广泛应用于气象研究的各个领域，应用数值模拟方法在青藏高原土壤冻融过程研究方面取得的成果更是日益丰富[15-17]。近年来，国家气候中心致力于研发新的模式版本，包括进一步改进物理过程的地球系统模式和高分辨率气候模式，并积极参与到世界气候研究计划（WCRP）发起的第六次国际耦合模式比较计划（CMIP6）中[16]。

在 CMIP6 计划中，国家气候中心通过改进陆面分量，进一步完善物理过程，发布了最新的BCC陆面过程模式，模式数据得到国内外学者的大量分析和使用。本文利用CMIP6 模式模拟的多层土壤温度数据，结合鄂陵湖草地站土壤观测资料和欧洲中心ERA5再分析资料，评估BCC陆面过程模式对高原地区冻土的模拟能力，以期加深对高原冻土冻融过程及其与气候变化相互作用的认识。

1 资料与方法

1.1 资料概况

研究资料包括CMIP6 模式模拟的多层土壤温度数据、鄂陵湖草地站土壤观测资料和欧洲中心ERA5再分析资料，具体说明如下：

（1） BCC-CSM2-MR 和其他CMIP6模式（见表1）通过历史气候模拟试验获取的多层土壤温度逐月数据（Temperature of Soil，TSL），研究时段为 1984～2014 年。由于BCC等模式资料缺少土壤温度逐日数据，选用2012～2013年逐3 h表面温度近似代替（Surface Temperature Where Land or Sea Ice，TSLSI）。上述数据来源于 CMIP6网站 (https://esgf-node.llnl.gov/projects/cmip6/)。

表1 CMIP6模式资料

（2）鄂陵湖草地站土壤观测资料来自中科院西北所。观测站点位于青海省果洛藏族自治州西北部玛多县鄂陵湖西北侧，地理坐标为 34.913°N、97.553°E，邻近擦泽村，海拔4274 m。该站点下垫面为高寒草甸，土壤粒径较粗，砾石和岩石碎块含量较高，周围地形开阔，属于典型的高寒半干旱大陆性气候，主要分布有季节冻土、岛状多年冻土和大片连续多年冻土。由于观测资料有较长时间的缺测，本研究筛选出2012年10月～2013年10月无缺测的土壤温度数据，测量深度分别是5 cm、10 cm、20 cm和40 cm。

（3）欧洲中心再分析资料（ERA5），水平分辨率为31 km，来自欧洲中心官方网站（https://www.ecmwf.int/）。

1.2 研究方法

冻融日期定义：陈渤黎等[18]指出可以将全年分为4个不同的冻融阶段（不考虑盐分对冰点降低的影响），即完全消融（日最低土壤温度＞0°C）、完全冻结（日最高土壤温度＜0°C）、冻结过程（土壤剖面处于冻结过程中）和消融过程（土壤剖面处于消融过程中）。在一个冻融阶段向下一个冻融阶段转变过程中，不能认为刚好满足下一个冻融阶段判定条件的日期即为下一个冻融阶段的开始日。因此，为了通过合理地描述土壤冻融日期来反映真实的冻融过程，在状态转变时若连续3 d均满足下一冻融阶段的判定条件，即认为这3天中的第一天为下一冻融阶段的起始日。

研究使用Mann-kendall趋势检验[19-21]、Sen斜率[22]估计方法评估青藏高原土壤温度气候倾向特征。Mann-kendall方法原本是一种检验两个变量序列之间相关性的方法，而把时间作为其中一个变量，即检验一个统计量与时间（如年份）之间的关系，就相当于检验该统计量的变化趋势[23]。Mann-kendall趋势检验和Sen 斜率估计均为非参数检验，适用于分析时间序列的单调增加或减小趋势。假设数据样本是独立和随机排列的，统计量计算公式如下：

式中sgn为符号函数，在此处即表示：

式中，xi、xj分别为i、j时刻的样本值，j＞i，n为样本量。此时提出原假设H0：时间序列x没有显著的单调变化趋势。当原假设不成立，且n≥10时，统计量S近似服从正态分布。

此时S的平均值为：

S的方差为：

式中，m为并列组（一组具有相同值的样本数据）的数目，ti为第i组的样本数。

此时再做转换：

式中，Zc服从标准正态分布。当Zc＞0时，序列有上升趋势；Zc＜0时，序列有下降趋势。考虑到正态分布的对称性，给出显著性水平α，则当Zc＞Z1-α/2或Zc＜Zα/2时，趋势通过显著性检验，拒绝原假设，即变化趋势显著。

Sen斜率估计统计量Qmed为时间序列中任意两个时刻样本间斜率的中间值，计算方法如下：

式中，xj、xk分别为j、k时刻的样本值，j＞k。

2 冻融日期模拟效果检验

从冻融过程来看，由图1可知在选取的模式中，总体上对冻融过程模拟效果最接近观测的是CMCCCM2-SR5和CMCC-ESM2，其次是NESM3，在所有模式中BCC- CSM2-MR的模拟效果好于大部分模式。从整个冻融过程的时间来看，观测值显示冻融过程从冻结初日（2012年10月5日）至完全消融初日（2013年4月15日）共持续192 d，BCC-CSM2-MR模拟的冻融过程时间比观测增加了76 d。在冻结过程阶段，观测值从10月5日开始至10月17日结束共12 d，而BCC-CSM2-MR模拟结果比观测明显增长46 d，冻结初日时间提前到9月21日，完全冻结初日日期推迟到了11月18日。在完全冻结阶段，观测值从10月17日开始至2月22日结束共计128 d，CNRM-CM6-1与观测值最为接近，而BCC-CSM2-MR模拟值与观测相差23 d，完全冻结初日日期和消融初日日期都有一定程度的推迟。在消融过程阶段，观测值从消融初日2月22日开始至完全消融初日4月15日结束共52 d，BCC-CSM2-MR模拟的消融过程天数与观测较接近，但消融初日推迟到4月13日，完全消融初日推迟到6月16日，比观测值落后两个月左右。通过以上分析，BCC-CSM2-MR模拟的完全冻结和消融过程阶段天数与观测值接近，但BCC-CSM2-MR对冻结过程阶段模拟效果不佳，持续时间增长46 d，导致完全冻结和消融过程阶段日期推迟近两个月。这可能是由于BCCCSM2-MR陆面模式物理参数化过程不完善，能量传输与实际有差别，导致温度下降更慢造成的。

图1 BCC等CMIP6 模式和观测资料反映的冻融过程

通过分析总天数指标（图1）可知：在冻结过程阶段，与观测值最接近的模式是CMCC-CM2-SR5、CMCC-ESM2、MIROC-ES2L 分别相差 11 d、16 d、22 d，而BCC-CSM2-MR 、CNRM-ESM2-1、CNRM-CM6-1、CESM2与观测值相差较大；在完全冻结阶段，与观测值最接近的模式是CNRM-CM6-1，仅相差4 d，其余模式皆有一定程度的偏差；与国际先进同类模式对比，CMCC-CM2-SR5与 MIROC-ES2L对消融过程阶段的模拟效果最佳，相差0 d，BCC-CSM2-MR次之，相差7 d。

通过分析日期指标（表2）可知：与观测冻结初日最接近的是CMCC-ESM2和CMCC-CM2-SR5，分别落后了10 d、12 d，其余模式模拟的冻结初日皆有所提前；除CMCC-CM2-SR5与 MIROC-ES2L外，所选模式对于完全冻结初日的模拟均有落后的现象；对于消融初日，CMCC-ESM2模拟结果与观测值最接近，相差6 d，其他模式都与观测值有一定差距；同国际各模式相比，CMCC-CM2-SR5对完全消融初日的模拟同观测值最为接近，相差14 d。

表2 BCC等CMIP6 模式和观测资料的冻融日期

3 冻土相关物理量模拟效果检验

3.1 冻土深度

冻土深度是青藏高原土壤冻融过程模拟检验的重要指标之一。将土壤温度为0°C所在的深度定义为土壤冻结深度，计算BCC-CSM2-MR等CMIP6模式数据、ERA5再分析数据与鄂陵湖西岸草地观测站点逐月土壤冻结深度的统计量。将所求得的相关系数（Correlation）、标准差（Standard Deviation）以及均方根误差（Root Mean Square Error）整合到泰勒图上（图2）。从变化趋势来看，E3SM-1-0、GFDL-ESM4和ACCESSCM2相对CMIP6其他模式及ERA5来说，对冻土深度变化趋势的模拟效果最佳，相关系数值分别为0.94、0.92、0.88，而BCC-CSM2-MR的相关系数为0.43。总体上看，各模式对于冻土深度随时间波动的模拟还有待提高，多数模式相关系数介于0.4～0.5（均通过了0.05水平的显著性检验）。这可能与各模式土壤层数有限且各不相同有关，各层土壤温度值的非线性变化也将影响插值函数的运算，导致各模式模拟的冻土深度有一定误差。从样本的离散程度来看，观测值的标准差为4 m，IPSL-CM6A-LR-INCA的标准差相对其他模式资料与观测值最接近，标准差为3.5 m，而BCCCSM2-MR标准差为1.7m。多数模式均方根误差值介于2～4 m，其中GFDL-ESM4相对CMIP6其他模式和ERA5，与观测值冻土深度偏差最小，均方根误差值为1.58 m，BCC-CSM2- MR均方根误差值为3.34 m，处于平均水平。

图2 CMIP6模式资料、ERA5再分析资料相对于实测土壤深度的泰勒分布

图3给出了BCC等CMIP6 模式、ERA5再分析及观测资料反映的冻土深度逐月变化特征。如图所示，多数CMIP6模式、ERA5再分析以及观测资料大致表现为单峰特征。2012年10月～次年4月或5月是土壤冻结时段。对于冻土深度增加阶段，从10月～次年2月，观测值冻土深度逐渐由0 m增加至11.24 m。BCC-CSM2-MR模式在这段冻结深度不断增加的过程中，从0.2 m上升至4.3 m，10月～次年2月与观测曲线最为接近，但冻土深度最大值相对观测值偏小且冻土深度增加，终止时间相对于观测推迟一个月。E3SM-1-0对土壤冻结时段变化趋势的模拟效果较好，冻土深度最大值出现月份与观测值一致，但数值上相对于观测值总体偏高。10月～次年1月，多数模式模拟的冻土深度变化趋势与观测结果较为接近。但在随后过程中，最大冻结深度的模拟结果与观测值均有一定偏差，冻土深度增加的终止时间相对于观测推迟一个月，这可能与各模式土壤层数有限且各不相同，以及各层土壤温度值的非线性单调特征有关。对于冻土深度降低阶段，从2～4月，观测值冻土深度总体呈下降趋势，从11.24 m降至基本消融。GFDL-ESM4模式的模拟结果与观测值最为接近，从12.32 m降至0 m。BCC-CSM2-MR、CMCC-ESM2等大多模式及ERA5均能模拟出冻土深度的降低趋势，但模拟最大深度比观测值偏低。综上所述，BCC-CSM2-MR基本能模拟出冻土深度变化曲线呈单峰的形式，模拟效果整体好于部分模式，但BCC-CSM2-MR模拟的最大冻深比观测值偏低，且冻土深度增加的终止时间相对于观测推迟一个月。

图3 BCC等CMIP6 模式资料、ERA5再分析资料及观测资料反映的冻土深度逐月变化

图4给出了ERA5再分析资料、BCC等CMIP6 模式资料反映的青藏高原冻土深度空间分布和模式集合平均结果。如图所示，多个模式资料和再分析资料反映的冻土深度空间分布与模式资料集合平均结果都较为接近，总体上表现为从高原外围向中部和西北部逐渐加深的特征，这可能与高原海拔西高东低、气温东高西低且南高北低导致高原西北部和中部冻土相对深厚有关[24]。对比CMIP6各模式和ERA5再分析资料的冻土最大深度区域（图4）可知，冻土最大深度差异显著，这是因为永久冻土区域冻土深度一般较大，各模式土壤层数有限、各层土壤温度非单调性均会影响最大冻土深度的插值结果。CMIP6各模式对于青藏高原西南部与柴达木盆地的冻土模拟也有一定差异。集合平均结果显示，青藏高原西南部冻土呈两侧高、中间低的分布态势，与该区域的地形分布特征一致，即藏北高原与冈底斯山脉高海拔地区间夹着相对低值区。CMCC-CM2-SR5、CMCC-ESM2、TaiESM1与集合平均结果比较接近，而BCC-CSM2-MR不能模拟出这种分布。除此之外，BCC-CSM2-MR对柴达木盆地冻土深度分布的模拟效果还有待提高，这可能与BCC-CSM2-MR模式网格分辨率较低，一个面积只有一个土壤列代表，小块和岛状多年冻土无法被模式准确模拟有关[25]。

图4 ERA5再分析资料（a）、BCC等CMIP6 模式资料（b～t.依次对应表3中各模式）和模式集合平均（u）反映的青藏高原冻土深度空间分布（单位：m）

3.2 土壤温度气候变化率

基于BCC等CMIP6 模式资料和ERA5再分析资料，提取青藏高原地区土壤温度，计算1985～2014年逐月土壤温度气候倾向率（Sen 斜率估计统计量 Q），详见表3。总的来看，上述资料在1985～2014年土壤温度均呈上升趋势，逐月气候变率基本为正值。集合平均的最高增温率出现在2月和7月；CNRM-CM6-1最高增温率与集合平均结果最为接近，2月为0.029°C·a-1，7 月为 0.04°C·a-1；BCC-CSM2-MR 最高增温率出现在2月和9月，但数值总体偏大。对于最低增温率，集合平均大值出现在1月、5月和11月，分别为0.022°C·a-1、0.024°C·a-1和 0.021°C·a-1，各模式模拟效果均不佳。

表3 青藏高原逐月土壤温度气候倾向率（单位： 10-2°C·a-1）

图5给出了青藏高原土壤温度气候倾向率空间分布（打点区域通过0.05水平的显著性检验）。集合平均结果显示：1985～2014年青藏高原土壤温度气候倾向率自东南向西北呈阶梯状上升的分布特征；高原北部存在部分相对低值区域，这些低值区域对应着高原山脉（昆仑山脉）之间的低海拔地区（柴达木盆地）。已有研究[26]表明，青藏高原具有强烈的山体效应，深刻影响着高原的水热分配，即由边缘越往内部温度越高，并伴随着变干趋势。除了ERA5和部分模式中存在土壤温度呈下降趋势的区域外（分别在高原西北部和东南部），大部分模式在高原地区的模拟结果与集合平均一致，土壤温度气候倾向率为正。CMCC-CM2-SR5能模拟出上述低值区域，与集合平均结果最为接近，但整体上数值有所偏低。BCC-CSM2-MR模拟结果自东向西呈阶梯状上升的分布形势，但东北部地区气候倾向率相较于集合平均偏低，而西部地区偏高。

图5 同图4，但为土壤温度气候变化趋势（单位：°C·a-1，加点区域表示通过了0.05水平的显著性检验）

4 结论与讨论

本文利用CMIP6 模式的多层土壤温度资料、欧洲中心ERA5再分析资料和鄂陵湖草地站土壤观测资料，评估了BCC陆面过程模式等对青藏高原土壤冻融过程的模拟能力，得到如下主要结论：

（1）BCC-CSM2-MR对土壤冻融日期的模拟效果好于大部分模式，其优点在于对冻融总天数，特别是对于消融过程阶段的模拟接近观测值，但BCC模式模拟冻结过程阶段的结束时间较晚，导致后续两个阶段（完全冻结阶段和消融过程阶段）日期都有所推迟。这可能与BCC模式物理参数化过程对相变过程描述不准确，从而使得土壤温度降低过程减缓有关。

（2）BCC-CSM2-MR基本能模拟出冻土深度变化曲线的单峰特征，对土壤冻结时段的前期模拟效果最佳。但BCC-CSM2-MR模式各统计量结果表现一般，最大冻深比观测值偏低，且冻土深度增加的终止时间相对于观测推迟一个月。BCC-CSM2-MR能模拟出从高原外围向中部和西北部逐渐加深的形势，但对青藏高原西南部两侧高、中间低的冻土深度分布特征不能很好地模拟。

（3）BCC-CSM2-MR 模拟1985～2014年土壤温度总体为上升趋势，从东向西呈阶梯状上升的分布特征，与集合平均较为接近。模拟的东北部地区气候倾向率相较于集合平均偏低而西部地区偏高，不能模拟出高原北部存在的部分相对低值区域。

虽然耦合全球气候系统模式能够定量地研究冻土当前状态、未来变化以及对气候变化的影响，但通过以上讨论发现，气候系统模式有限的分辨率和参数化方案的差异仍给模式结果带来了很大的不确定性[27-29]。随着中国区域高分辨率观测数据的收集、格点化观测数据集的建立及其在高分辨率模式模拟检验中的广泛应用，未来将有望使用格点化观测资料对青藏高原陆面过程进行研究。