多年冻土约占北半球陆地面积的1/4。全球多年冻土中储存着约18300×10⁸吨的土壤有机碳，相当于大气中二氧化碳含量(8280×10⁸吨)的2.2倍^[1~2]。北半球多年冻土中还存储着约11370×10⁹~36550×10⁹m³的地下冰，相当于全球海平面升高3~10cm的水含量^[3]。此外，表层冻土的季节性冻融过程与地表生态、水文及能量收支有着紧密的联系^[4~9]。冻土消融会导致土壤有机碳释放、地下冰消融以及冻土季节性冻融规律的改变，进而对气候、水文与水资源、冻土工程稳定性以及地表生态环境产生重要影响^{[4, 10]}。因此，在当前气候变暖以及将来气候可能进一步变暖的背景下^[11~14]，冻土消融及其影响成为气候变化研究领域的热点问题。

青藏高原上广泛分布着多年冻土和季节冻土。多年冻土约1.06×10⁶km²，季节冻土约1.46×10⁶km^2[15]。由于青藏高原位于北半球中低纬度，因此其多年冻土属性较暖，对气候变化尤为敏感；研究表明气候变暖引起的青藏高原多年冻土消融速度是高纬度多年冻土的1.4倍^[16]。此外，作为“亚洲水塔”，青藏高原是东南亚一些主要河流的源头。被誉为“天路”的青藏铁路的约50 %路段位于多年冻土之上。多年冻土退化会对亚洲水资源供应、青藏铁路的稳定性等产生重要影响^[17]。因此，开展青藏高原冻土变化及其影响研究尤为重要。

作为冻土研究中最为关键的参量，土壤温度对冻土变化具有重要指示作用，已经引起了越来越多的关注。例如，国际上启动的“国际极地年计划”就旨在调查多年冻土的土壤温度状况^[18]。基于野外钻孔观测资料，青藏高原土壤温度的日、月和年际变化(1996~2006年)特征已经得以研究^[19~21]；基于气象站观测资料，土壤温度研究的时间尺度被拓展到了年代际尺度(1960~2005年)^{[22, 23]}。但由于青藏高原观测站点稀少，基于观测的研究在区域代表性上受到了限制。于是学者们开始尝试评估和利用卫星遥感和再分析资料研究青藏高原土壤温度，突破区域限制^[24~26]。然而，再分析资料主要开始于1979年，并且土壤深度基本在3.5m以内，难以研究更长时间和更深层的土壤温度变化。值得庆幸的是，随着模式和计算能力的不断改进和发展，数值模拟方法已经能够弥补这些不足。

本研究利用最新版本的通用陆面模式(CLM4.5) 对过去百年青藏高原0~35m土壤温度变化进行了研究，旨在从长时间和大空间尺度上揭示青藏高原土壤温度的演化规律，为青藏高原冻土及陆地生态系统变化研究提供依据。

本文所用的大气强迫场为模式自带的CRUNCEP资料^[27]，包括气温、降水、气压、风速、比湿和辐射变量。空间分辨率为0.5°，时间分辨率为6小时。资料时期为1901~2010年，覆盖了全球范围。CRUNCEP资料融合了CRU(Climate Research Unit，0.5°×0.5°)观测资料^[28]和NCEP(National Centers for Environmental Prediction，2.5°×2.5°)再分析资料^[29]。对与冻土模拟相关较紧密的气温、降水和辐射变量，CRUNCEP资料在发展过程中利用月平均CRU观测资料对NCEP资料进行了校正，保证了这些变量在月尺度上的准确性。当前CRUNCEP资料有5个版本，本文用的是第4个版本，这套资料已经被用于驱动通用陆面模式进行植被生长和蒸散发等研究^{[30, 31]}。

土壤温度再分析资料ERA-Interim(European Centre for Medium-Range Weather Forecasts Re-Analysis Interim)^[32]和CFSR(National Oceanic and Atmospheric Administration Climate Forecast System Reanalysis)^[33]用于与模拟结果进行对比。ERA-Interim资料的空间分辨率为0.75°，CFSR的空间分辨率为0.5°，它们的时期均为1979~2016年。这两套再分析资料没有直接提供对比深度(0.2m)的土壤温度，本文中通过线性插值方法获得。我们之所以用这两套再分析资料是因为它们已被证实在青藏高原上相比其他再分析资料具有更高的精度^{[25, 26, 34]}。我们将模拟结果与它们进行对比来说明模拟结果相比再分析资料的优劣。

气象站观测的土壤温度资料用于验证模拟结果。这套资料观测了8层土壤的土壤温度，深度分别为0.05m、0.1m、0.15m、0.2m、0.4m、0.8m、1.6m和3.2m，其中0.4m及以上的浅层观测资料具有较长观测时期，从20世纪60、70年代开始观测，数据较为连续；0.4m以下土壤资料基本从2004年开始才有观测数据。考虑到利用尽可能长序列的观测资料进行验证模式，我们仅采用了位于冻土活动层之内的0.2m深处土壤观测资料。综合观测时期和数据缺测等因素，我们在青藏高原选取了20个气象站点，其分布如图 1所示，坐标信息如表 1所示。验证时，先将所有站点的观测资料和与其对应模拟格点的模拟资料做平均，然后将二者进行对比。

图 1(Fig. 1)

图 1 研究区域(25°~40.5°N，75°~105°E)和观测站点分布图 圆形代表气象站站点，五角星代表野外钻孔监测站点 Fig. 1 Distribution of study area(25°~40.5°N, 75°~105°E) and observation sites. Circles represent meteorological stations and five-pointed stars represent field borehole monitoring sites

表 1(Table 1)

表 1 气象站站点和野外钻孔监测站点的信息 Table 1 Information of meteorological stations and field borehole monitoring sites 站名 纬度/°N 经度/°E 海拔/m 会泽 26.42 103.28 2110 贵州 36.03 101.43 2237 威宁 26.87 104.28 2238 西宁 36.72 101.75 2295 维西 27.17 99.28 2326 丽江 26.87 100.22 2392 波密 29.87 95.77 2736 格尔木 36.42 94.90 2808 恰卜恰 36.27 100.62 2835 林芝 29.67 94.33 2992 大柴旦 37.85 95.37 3173 中甸 27.83 99.70 3277 昌都 31.15 97.17 3306 兴海 35.58 99.98 3323 泽当 29.25 91.77 3552 拉萨 29.07 91.13 3649 日喀则 29.25 88.88 3836 隆子 28.42 92.47 3860 江孜 28.92 89.60 4040 狮泉河 32.05 80.08 4279 FH1 34.69 92.89 4938 WD2 35.13 93.04 4707 CM2 35.40 93.54 4482 Wuli 34.47 92.73 4625

表 1 气象站站点和野外钻孔监测站点的信息 Table 1 Information of meteorological stations and field borehole monitoring sites

野外钻孔观测资料用于验证深层(1m和6m)土壤温度模拟结果。考虑选取具有较长观测时期(10年以上)的站点和数据的可获得性，收集了4个站(FH1、WD2、CM2和Wuli)的观测资料，其分布如图 1所示，坐标信息如表 1所示。数据的观测时期为1996~2010年，部分站的资料要短于该时期。这些数据主要来自文献[20, 35]。有关观测资料的观测方法、后处理等详细信息可见相关文献[20, 35]。

本研究利用了最新版本的通用陆面模式CLM4.5(Community Land Model，version 4.5)。模式明确考虑了冻土冻融过程，包含了一个冰点降低方程^[36]，允许液态水和冰在温度低于0℃时可以共存，增强了水的下渗能力。模式还考虑了土壤有机物的水热属性，减小了模拟结果的暖偏差^[37]。此外，土壤层被拓展到了15层、约50m深，克服了深层多年冻土的热绝缘问题^[37]。这些改进使得CLM4.5成为当前较为成熟的模拟冻土的模式之一^[38]。需要输入的大气强迫资料包括气温、降水、气压、风速、比湿和向下短波辐射，时间分辨率为6小时。CLM模式目前已经广泛用于冻土及陆面模拟研究^[38~41]。

利用CRUNCEP资料驱动CLM4.5模式进行了模拟试验。植被和土壤质地使用了模式自带的数据集。植被由中等分辨率成像分光辐射计(MODIS)产品反演获得^[42]。土壤质地基于国际岩石圈—生物圈计划(IGBP)土壤数据集发展而成^[43]。模拟输出的时间分辨率为每天，模拟的空间分辨率为0.5°。模拟区域为全球范围，本文选取青藏高原区域的土壤温度模拟数据进行分析。模式的初始化过程为利用1979年的大气数据驱动CLM4.5循环模拟了100年，再用1901年的大气数据驱动CLM4.5继续循环模拟了5年。在初始化过程结束期，所有层区域平均土壤温度的年际变化都低于0.003℃。将此时的模式状态做初始场，开展瞬时模拟，模拟时期为1901~2010年。

纳什效率系数(Nash-Sutcliffe Efficiency, 简称NSE)和相关系数用于评估模拟结果与观测值的吻合程度。NSE由Nash和Sutcliffe^[44]提出，计算如下：

(1)

公式(1) 中Q_obs代表观测值，Q_sim代表模拟值，Q_obs代表观测值的平均，T为总样本数。NSE表示模拟和观测值接近1︰1线的程度，其对极值敏感。纳什效率系数的范围为-∞到1，1表示观测和模拟值完全相同。系数越接近1，表明观测和模拟值越接近；系数接近0，表示模拟值接近观测值的平均水平，总体结果可信；系数远小于0，则模拟结果不可信。

模式验证时，首先利用最邻近插值方法将模拟结果插值到观测站点上，然后再与站点观测值进行对比。土壤温度增温趋势采用普通的最小二乘法和线性拟合方法计算获得。变化趋势的显著性利用t检验方法。

模式模拟出了观测的0.2m土壤温度年际变化和趋势，模拟效果比再分析资料合理(图 2)。模拟和观测结果的相关系数为0.92，NSE为0.82。它们的线性趋势十分接近，分别为0.27℃/10年和0.28℃/10年。与再分析资料的相比，不论是相关系数、NSE、还是趋势，模拟结果都要明显优于两套再分析资料(表 2)。这一结果表明这套模拟资料具有相对再分析资料更高的精度，意味着它们在将来可以用于与土壤温度相关科学问题的研究。

图 2(Fig. 2)

图 2 1972~2005年0.2m土壤温度变化(相对1981~2000年)模拟结果和再分析资料与观测值的对比 CC代表相关系数，“*”代表相关系数通过95 %显著性水平，NSE代表纳什效率系数，趋势的单位为℃/10年 Fig. 2 Comparison of changes in simulated and reanalysis soil temperature at a depth of 0.2m to observations during the period 1972~2005 relative to 1981~2000. CC represents correlation coefficient, "*" represents the correlation coefficient exceeding 95 % significance level, and NSE represents the Nash-Sutcliffe efficiency

表 2(Table 2)

表 2 1979~2005年0.2m土壤温度本文模拟结果、ERA-Interim和CFSR资料与观测值对比的统计 Table 2 Statistics of comparison of simulations, ERA-Interim and CFSR data to observed soil temperature at a depth of 0.2m during the period 1979~2005 观测 模拟 ERA-Interim CFSR 相关系数 0.91 0.74 0.73 NSE 0.78 0.50 0.07 趋势(℃/10年) 0.31 0.34 0.09 0.39

表 2 1979~2005年0.2m土壤温度本文模拟结果、ERA-Interim和CFSR资料与观测值对比的统计 Table 2 Statistics of comparison of simulations, ERA-Interim and CFSR data to observed soil temperature at a depth of 0.2m during the period 1979~2005

模拟的深层(1m和6m)土壤温度变化也表现合理(图 3)。基于钻孔观测资料的验证显示，模拟结果较好地再现了观测的变化，深层较浅层更好。相关系数在0.52~0.94之间，NSE在-0.57~0.88之间。相对FH1、WD2和Wuli站，CM2站的模拟效果较差。由于这种对比是模拟格点平均值与单个站点观测值之间的对比，因此存在尺度不匹配问题。加之青藏高原地形复杂，使得尺度不匹配问题更为突出，这可能放大了某些地区模拟结果的偏差。这一原因和图 3显示的模拟性能证明了模拟结果的合理性。

图 3(Fig. 3)

图 3 1981~2010年1m(a，c，e，g)和6m(b，d，f)土壤温度变化(相对各站点的观测时期)模拟结果与各站点观测时期观测值的对比 CC代表相关系数，“*”代表相关系数通过95 %显著性水平，NSE代表纳什效率系数，趋势的单位为℃/10年 Fig. 3 Comparison of simulated change in soil temperature at a depth of 1m(a, c, e, g)and 6m(b, d, f)to observations during the period 1981~2010 relative to observation period of each site. CC represents correlation coefficient, "*" represents the correlation coefficient exceeding 95 % significance level, and NSE represents the Nash-Sutcliffe efficiency

整个青藏高原1m深度土壤温度都显著增加，西北部增温高于东南部，冬季增温高于夏季(图 4)。从空间分布上看，整个青藏高原土壤温度都显示为增加趋势，增温最高的区域在高原西北角，增温最低区域在高原东南部(图 4a)。从区域平均线性变化来看，总体上1901~2010年青藏高原土壤温度呈增加趋势，增加率为0.08℃/10年。1979~2010年土壤温度快速增加，增加率为0.38℃/10年。可以看出，土壤温度在20世纪40年代出现过一次明显增加，这与观测的全球气温变化相一致^[11]。从季节来看，冬季增温高于夏季，1901~2010(1979~2010) 年冬季和夏季的增温率分别为0.10(0.41)℃/10年和0.08(0.36)℃/10年。

图 4(Fig. 4)

图 4 1901~2010年年平均(a，b)、夏季(c，d)和冬季(e，f)1m土壤温度的趋势(℃/10年)分布(a，c，e)和区域平均线性变化(相对1981~2000年)(b，d，f) “+”代表趋势通过95 %的显著性检验，“*”代表相关系数通过95 %显著性水平 Fig. 4 Distribution(a, c, e)of trend(℃/decade) in soil temperature at a depth of 1m during the period 1901~2010 and the related area-average linear change(b, d, f)relative to 1981~2000 for annual(a, b), summer(c, d)and winter(e, f). "+" represents the trend exceeding 95 % significance level, and "*" represents the correlation coefficient exceeding 95 % significance level

15m深度土壤温度增加的时空特征与1m深度相似，但增加率相对较低，季节变化不明显(图 5)。年平均地温指地下土壤温度季节变化幅度为0的深度的土壤温度，它能很好的指示多年冻土对长期气候变化的响应。这一温度在青藏高原上一般指10~15m深度的温度^[20]。因此，我们分析了15m土壤温度变化。可以看到，15m土壤温度的增加模态与1m的相似，都为西北部高于东南部，1979~2010年增温更快。不同的是15m土壤温度的增加率低于1m深度。1901~2010(1979~2010) 年15m土壤温度的增加率为0.05(0.27)℃/10年，无季节差别。

图 5(Fig. 5)

图 5 1901~2010年年平均(a，b)、夏季(c，d)和冬季(e，f)15m土壤温度的变化趋势(℃/10年)分布(a，c，e)和区域平均线性变化(相对1981~2000年)(b，d，f) “+”代表通过趋势95 %的显著性检验，“*”代表相关系数通过95 %显著性水平 Fig. 5 Same as Fig. 4, but for soil temperature at a depth of 15m

图 6(Fig. 6)

图 6 青藏高原经向(a)和纬向(b)平均的1901~2010年土壤温度变化趋势(℃/10年)分布 深度中的“-”号代表向下，所有趋势都通过了95 %的显著性检验 Fig. 6 Distribution of trend in soil temperature(℃/decade)during the period 1901~2010 for meridional (a) and latitudinal (b) mean over the Tibetan Plateau. "-" in the value of depth represents underground. Trends for all grids exceeding 95 % significance level

青藏高原土壤温度变化具有明显的深度、经度和纬度依赖性(图 6)。随着深度从0增加到35m，土壤增温逐渐减弱。这可能与深层多年冻土消融吸热，抑制增温有关。随着经度增加，土壤增温逐渐减弱(图 6a)。随着纬度增加，土壤增温的规律是：先增加，而后略减小(图 6b)。

图 7(Fig. 7)

图 7 1901~2010年整个青藏高原地区(a，d)、多年冻土区(b，e)和季节冻土区(c，f)1m(a，b，c)和15m(d，e，f)土壤温度与气温变化(相对1981~2000年)的关系 敏感性指土壤温度对气温变化的敏感性，以子图(a)为例，气温增加1℃，对应土壤温度增加0.79℃ Fig. 7 Relationship between soil temperature and air temperature changes at depths of 1m(a, b, c) and 15m(d, e, f)over entire Tibetan Plateau(a, d), permafrost(b, e), and seasonally frozen ground(c, f)regions during the period 1901~2010 relative to 1981~2000

1m深度土壤温度与气温的相关关系明显好于15m深度(图 7)。1m深度土壤温度对气温变化的敏感性为0.79℃/℃(即气温增加1℃，1m深度土壤温度增加0.79℃)，明显高于15m深度的0.37℃/℃(图 7)，这表明浅层土壤温度受气候变化影响比深层更大。就多年冻土和季节冻土区分开来看，多年冻土区1m(15m)土壤温度增加率为0.08(0.07)℃/10年，高于季节冻土区的0.06(0.05)℃/10年(图略)。这里的多年冻土和季节冻土区的划分根据青藏高原冻土图^[45]。对气温变化的敏感性，多年冻土区也高于季节冻土区(图 7)，1m深度多年冻土温度对气温变化的敏感性为0.80℃/℃，季节冻土温度为0.78℃/℃，这表明多年冻土热状况对气温变化较季节冻土更为敏感。

本研究基于最新版本的通用陆面模式CLM4.5对过去百年青藏高原土壤温度演化进行了模拟研究，得到如下结果：

(1) 利用CRUNCEP驱动CLM4.5能合理地模拟出观测的青藏高原土壤温度的年际变化(相关系数为0.92，纳什效率系数为0.82) 和线性趋势(模拟和观测的趋势分别为0.27℃/10年和0.28℃/10年)，其结果要明显优于再分析资料。

(2)1901~2010年青藏高原1m和15m土壤温度增加率分别为0.08℃/10年和0.05℃/10年，1979~2010年土壤温度增加更快，增加率为0.38(1m)℃/10年和0.27(15m)℃/10年。

(3) 土壤增温特征为青藏高原西北部增温高于东南部、冬季增温高于夏季、浅层增温高于深层、增温随经度增加而减小、增温随纬度增加先增加而后略减小。

(4) 土壤温度增加主要受气温控制，浅层(1m)土壤温度对气温变化较深层(15m)更敏感，其敏感度分别为0.79℃/℃和0.37℃/℃；多年冻土温度对气候变化较季节冻土更敏感，其1m(15m)深度敏感性分别为0.80℃/℃和0.78℃/℃(0.38℃/℃和0.36℃/℃)。

上述结果对认识过去百年青藏高原冻土及陆地生态系统变化具有一定的指示意义。基于气象站和野外钻孔观测资料的验证证明了这些模拟结果的合理性。同时，对比发现本模拟结果要比再分析资料更接近观测的土壤温度变化，并且这套资料具有较长时期(1901~2010年)、较高时空分辨率(每天和0.5°)和较深深度(0~35m，15层)的优势，这些使得在将来可以广泛地运用这套模拟结果去研究与青藏高原土壤温度相关的各种变化(比如冻土)、联系(比如遥相关)及物理机制。该套模拟数据可以联系论文作者获得。

由于本研究利用了在冻土模拟方面较为成熟的最新版本的通用陆面模式，同时采用了模式自带的CRUNCEP大气资料。CRUNCEP资料融合了CRU观测资料，且分辨率较高，在区域陆面模拟方面具有一定的优势。这些优势在一定程度上提高了本文模拟结果的可信度。由于青藏高原地形十分复杂，对模拟分辨率要求较高。本研究中的分辨率为0.5°，这一分辨率对青藏高原而言可能仍较粗，给模拟结果带来了不确定，将来可以采取更高分辨率的模拟试验研究。此外，研究表明CLM模式中的植被在青藏高原上存在一定的误差^[46]，这可能在一定程度上影响了模拟的准确度。将来可以采用国内的《中国1︰1000000植被图》^[47]数据集对模式地面和土壤数据进行校正，减少模拟结果的不确定。

鉴于土壤温度对冻土及陆面过程、生态等的重要性，之前学者已经对它们的日变化、月变化、季节变化、年际变化、年代际等进行了研究^[19~23]。本研究在更长时间尺度(110年)上对青藏高原土壤温度变化进行了数值模拟研究。揭示了过去百年青藏高原土壤温度演化的时空特征，阐明了土壤温度演化与土壤深度、纬度、经度和气温变化的联系。将来的研究可以集中在过去百年青藏高原土壤湿度和地表水热交换过程变化及其机理的模拟研究方面。

致谢: 感谢同行评审专家和编辑部老师提出的宝贵修改意见!本研究受国家自然科学基金项目(批准号：41405087和41401086) 和公益性行业(气象)科研专项项目(批准号：GYHY201406001) 共同资助。