1. 引言
1.1. 研究背景
青藏高原，位于亚洲大陆中南部、中国西南部，是全球平均海拔最高的高原，是世界上最大的高海拔中低纬度多年冻土区，也是我国最大冻土分布区之一[1]。江河源区位于青藏高原腹地，是东亚气候变化的敏感区之一，是我国重要的水资源补给地之一。江河源区是黄河、长江和湄公河的源头，对国家水资源安全、生态环境保护与区域可持续发展有重大而深远的影响[2]。
冻融是由于季节或昼夜热量变化在表土及以下一定深度形成的反复冻结–解冻的土壤物理过程，这种现象在高纬度和高海拔地带低温生态系统普遍存在[3]。冻土中土壤水分随温度变化而产生的由融化到冻结、由冻结到融化的过程被称为冻融循环[1]。土壤冻融过程是青藏高原陆面过程中最突出的特征之一，研究土壤冻融过程的时空分布特征是理解其变化机理的基础[4]。
1.2. 研究现状
青藏高原的动力作用和热力作用，在不同程度上影响着高原及周围地区的气候变化[5]。青藏高原近地表土壤冻融循环极大地影响着土壤和大气之间的物质交换[1]。土壤的冻结和融化过程，会导致水热特性及变化规律发生改变，会导致地气能量的交换，会导致寒区水循环及生态环境的变化[6]。高原地表非绝热加热受土壤水热传输的影响，尤其是在地表过程快速变化的春季，地表融雪、融冻过程对地表非绝热加热以及地气之间的能量和水分交换产生显著影响[7]。冻融交替变化的过程中影响融雪水下渗，进而改变土壤水热平衡状态，且对陆地生态系统、水文过程及碳循环、气候系统等也具有深刻影响[8]。
Wang [9]等(2015)利用636个气象站点的近地表温度，分析了1956~2006年中国地区近地表土壤冻结首日、终日、冻结期和冻结天数对气候变化的响应，可比较准确反应青藏高原的土壤冻融过程，但由于青藏高原上站点稀缺，数据不连续，因此研究区域有限[10]。
对于整个青藏高原地区长时间尺度土壤冻融过程的时空分布研究还较少，使用在青藏高原地区适用性好的再分析资料可以研究大范围地区土壤冻融过程的、长期变化。为了找到适用性好的再分析资料[11]，陈渤黎等(2014)在若尔盖站冻融期对陆面模式CLM进行了模拟性能的检验，结果表明CLM4.0较CLM3.0在模拟结果上有了一定提高。Yang [12]等(2018)改进了陆面模式CLM4.5的冻融参数化方案，使得模式对于冻融过程中土壤温湿度的模拟结果有了改进。Deng [13]等(2020)基于改进的冻融参数化方案，对比分析了陆面模式CLM4.5和CLM5.0对青藏地区水热传输的模拟效果，进一步验证了该方案的普适性。
在国外，Menzel [14]等(2003)利用41个气象站数据，分析德国近地表土壤冻结时间随气温升高而缩短。Henry (2008)利用加拿大31个站点的数据分析近地表土壤冻融循环机制，结果得到近地表土壤冻结天数逐渐降低。Anandhi [15]等(2013)利用23个站点分析了美国堪萨斯州长时间序列近地表冻融变化，结果表明，近地表土壤冻结首日和冻结终日分别推迟和提前的趋势[10]。
1.3. 研究意义
全球变暖背景下，高原冻土出现不同程度的退化，其退化与生态环境、水分循环等息息相关，和气候系统间存在相互关系[16]。由于气候变暖，青藏高原也有明显变暖的趋势。冻土环境变化对陆地表面与大气之间的物质与能量交换、陆地景观格局演化的影响十分剧烈[17]。研究土壤冻结深度具有非常重要的意义，然而在过去的研究中，很少涉及土壤冻结深度的研究[18] [19]，主要是因为长时间序列的观测数据较少，只能够获取有限区域的土壤温度，阻碍土壤冻结深度的评估，进而影响分析其气候动力机制[20] [21]。目前，有些研究已经对江河源区的冻土特征进行了初步探索，但是对于江河源区土壤冻融循环的时空变化特征，尚缺乏系统性和全面性的研究。通过研究，加深理解江河源区冻土演变循环过程，有助于提高对其陆面过程及气候变化的认识。
江河源区作为我国重要的水源涵养区和生态屏障[22]，研究青藏高原江河源区土壤冻融循环时空变化特征，是理解其分布和变化机理的重要基础，对于理解生态环境的气候变化响应具有重要意义，对于理解全球和局部气候、水汽循环、陆地能量平衡、高原区域水文循环、植被生长和土壤呼吸等过程的生态系统循环以及人类活动具有重要影响[23]，对农业生产、土壤资源的有效利用以及生态气候与水文环境的预测具有重要的指示作用[24]。摸清冻土的分布格局及其时空变化对开展自然科学研究、保障生态环境和人类生产活动安全具有重要意义。
2. 资料和方法
2.1. 研究区概况
青藏高原是亚洲内陆高原，也是世界海拔最高、中国最大的高原，地理位置介于25˚N~40˚N、68˚E~105˚E之间，是全球最大的中低纬度多年冻土区[25]。土壤冻融是青藏高原地区独特的水文现象，本质是大气环境波动下引起的季节冻土活动层中的水热动态变化，是土壤水分在固、液两相之间转换的结果[26]。
青藏高原江河源区，是黄河、长江和湄公河的源头，参考文献，以30˚N~37˚N和88˚E~104˚E区域为研究对象[27]，研究这个范围内江河源区土壤冻融循环，针对该指定江河源区范围开展之后的讨论研究，江河源区范围如图1红框范围所示。江河源区位于青藏高原中心腹部地区其中包括以那曲为中心的南羌塘大盆区、以昌都为中心的藏东高山峡谷区、以玉树为中心的三江源区。那曲亚寒带地区，干旱且昼夜温差大，大风地区，昌都气候以寒冷为主，夏季温和湿润冬季寒冷干燥，有充足的阳光，降水集中分布不均匀，主要集中在夏季，玉树为高寒气候，冷季长达7~8月，年均温3.9℃ [28]。
Figure 1. Distribution map of rivers and observation stations in the source region of the Yangtze River on the Qinghai-Tibet Plateau
图1. 青藏高原江河源区及观测站点分布图
2.2. 资料概况
本文使用的数据包括1979年~2018年空间分辨率为0.1˚ × 0.1˚ CLM5.0模拟江河源区不同深度土壤温度资料和气温资料，美国地球物理数据中心的海拔数据。计算CLM5.0在江河源区的适用性涉及的站点如表1，并在图1中标注个站点具体位置。
Table 1. Latitude and longitude information of stations on the Qinghai-Tibet Plateau
表1. 青藏高原站点经纬度信息
| 站点 |
纬度(˚N) |
经度(˚N) |
| CST03 |
33.54 |
101.58 |
| CST04 |
33.46 |
102.43 |
| NST01 |
33.53 |
102.08 |
| NST21 |
33.54 |
102.08 |
| NST22 |
33.53 |
102.1 |
| NST25 |
34 |
101.59 |
| 称多 |
33.40 |
97.30 |
| 曲麻莱 |
34.58 |
95.35 |
| 雁石坪 |
33.58 |
92.07 |
| 玛沁 |
34.80 |
99.70 |
| 天峻 |
38.16 |
99.16 |
| 那曲 |
31.15 |
91.74 |
| 纳木错 |
30.77 |
90.96 |
| 亚东 |
27.90 |
89.26 |
2.3. 研究方法
本文验证CLM5.0土壤温度数据的适用性，其次使用通过验证的CLM5.0数据，对近四十年土壤冻融循环进行时空变化分析。最后，由处理数据得到的土壤冻融关键参数研究与气候数据(年平均气温，海拔高度)的响应。
2.3.1. 模拟资料的验证方法
适用性的研究过程中，计算相关系数、均方根误差、平均偏差。最常用的相关系数是皮尔逊相关系数(Pearson Correlation Coefficient)，其计算公式如下：
(1)
其中，r是皮尔逊相关系数，xi是观测值的第i个值，yi是模拟值的第i个值。
均方根误差(RMSE)是衡量预测值与实际值之间差异的标准之一。均方根误差表示预测误差的标准偏差，数值越小，表示预测模型的精度越高。偏差(Bias)是衡量估计值或预测值与实际值之间差异的指标。偏差表示估计量的期望值与真实参数值之间的差异。RMSE和Bias的计算公式如下：
(2)
(3)
其中，
是第i个数据点的实际值，
是第i个数据点的预测值，n是数据点的数量。
2.3.2. 线性趋势分析和相关性分析
使用通过验证的CLM5.0数据，对近四十年土壤冻融循环进行时空变化分析。研究近40年冻结开始时间、融化开始时间和冻结持续时间的变化趋势，用线性趋势描述，线性倾向表示不同参量的正负相关[29]。
(4)
其中，x是年份，y是土壤冻融循环关键参量。
相关性分析是统计学中用来衡量两个或多个变量之间关系的一种方法。用于探索土壤冻融过程关键参量与气温、海拔之间是否存在关联，以及关联的方向和强度。但相关性不等于因果性，即使两个变量之间存在强相关关系，也不意味着一个变量的变化导致另一个变量的变化。
2.3.3. 研究方法
为了将一个连续的冻结期包含在一年内，且综合青藏高原独特的地质气候特征，将当年的7月1日至次年的6月30日定义为一个完整的冻融循环年[30]。定义方法参考Wang等(2019)，土壤冻融过程关键参量有冻结开始时间、融化开始时间和冻结持续时间[4]。在一个冻融循环，从7月1日开始，为了避免随机变化和小扰动，定义连续三天土壤温度小于0℃，三天的第一天与7月1日之间的天数，称为冻结开始时间；在冻结开始时间之后，连续三天土壤温度大于0℃，定义三天的第一天与7月1日之间的天数，称为融化开始时间；冻结开始时间与融化开始时间之间的天数为冻结持续时间[31]。
3. 研究结果
3.1. 验证CLM5.0土壤温度数据适用性
表2为不同土壤深度CLM5.0土壤温度资料与观测数据对比，计算得到的相关系数均在0.9以上，部分地区计算得到的均方根误差偏大，平均偏差正常。表3为CLM5.0土壤温度资料与站点5 cm土壤温度观测数据对比，计算的相关系数、均方根误差和平均偏差，相关系数大多在0.8以上，也有部分地区计算得到的值偏大，在亚东地区计算得到的均方根误差和平均偏差合理，表现最好。总体来看，模拟CLM5.0土壤温度数据适用性较好。
Table 2. Comparison between CLM5.0 soil temperature data and observational data
表2. CLM5.0土壤温度资料与观测数据对比
| 站点名 |
时间 |
|
5 CM |
10 CM |
20 CM |
40 CM |
80 CM |
| CST03 |
2014/10/22~2015/6/30 |
R |
0.95** |
0.94** |
0.90** |
0.79** |
0.58** |
| RMSE/℃ |
2.05 |
2.39 |
2.97 |
4.12 |
5.28 |
| Bias/℃ |
0.39 |
0.58 |
0.20 |
0.23 |
−0.34 |
| CST04 |
2014/10/22~2015/6/30 |
R |
0.96** |
0.95** |
0.92** |
0.91** |
0.82** |
| RMSE/℃ |
1.84 |
2.23 |
2.70 |
2.88 |
3.94 |
| Bias/℃ |
−0.23 |
−0.61 |
−0.42 |
−0.39 |
−1.08 |
| NST01 |
2014/10/22~2015/6/30 |
R |
0.96** |
0.95** |
0.93** |
0.84** |
0.59** |
| RMSE/℃ |
1.91 |
2.19 |
2.83 |
3.69 |
5.78 |
| Bias/℃ |
−0.01 |
−0.68 |
−1.51 |
−0.82 |
−2.48 |
| NST22 |
2014/10/22~2015/6/30 |
R |
0.93** |
0.90** |
0.93** |
0.82** |
0.63** |
| RMSE/℃ |
2.43 |
2.97 |
2.44 |
3.84 |
5.02 |
| Bias/℃ |
0.29 |
0.13 |
0.14 |
−0.10 |
−0.36 |
| NST21 |
2014/10/22~2015/6/30 |
R |
0.94** |
0.93** |
0.89** |
0.80** |
0.63** |
| RMSE/℃ |
2.22 |
2.47 |
3.04 |
3.91 |
5.05 |
| Bias/℃ |
−0.37 |
0.04 |
0.11 |
−0.23 |
−0.57 |
| NST25 |
2014/10/22~2015/6/30 |
R |
0.95** |
0.94** |
0.89** |
0.79** |
0.60** |
| RMSE/℃ |
2.14 |
2.51 |
3.24 |
4.21 |
5.35 |
| Bias/℃ |
−0.32 |
−0.45 |
−0.53 |
−0.62 |
−1.03 |
注：**表示计算的相关系数通过99%的信度检验。
Table 3. Comparison between CLM5.0 soil temperature data and 5 cm soil temperature observational data at stations
表3. CLM5.0土壤温度资料与站点5 cm土壤温度观测数据对比
|
称多 |
曲麻莱 |
雁石坪 |
玛沁 |
天峻 |
那曲 |
纳木错 |
亚东 |
| R |
0.92** |
0.94** |
0.94** |
0.93** |
0.95** |
0.93** |
0.92** |
0.94** |
| RMSE/℃ |
2.70 |
3.88 |
2.89 |
3.86 |
2.65 |
4.32 |
4.28 |
2.37 |
| Bias/℃ |
−0.63 |
−2.67 |
−0.09 |
−0.67 |
−0.24 |
−3.28 |
2.95 |
−0.38 |
注：**表示计算的相关系数通过99%的信度检验。
3.2. 青藏高原土壤冻融循环过程关键参量的空间分布特征
根据以上定义，得到图2江河源区1979~2018年土壤冻融循环冻结开始时间(a)、融化开始时间(b)和冻结持续时间(c)对应的空间分布图。通过分析该区域冻融情况及其空间分布规律，图2(a)中江河源区冻结开始时间，大致为100~140天，集中在10月~11月份。图2(b)中融化开始时间为230~250天，大部分在3~4月。图2(c)中冻结持续时间大范围为100~150天，青藏高原部分边界地区冻结持续时间可能偏大或偏小。总体来说，江河源区具有西北向东南，冻结开始时间增长，冻结延迟，融化开始时间减短，融化提前，冻结持续时间减少的总体趋势。在江河源区北部，有部分地区与周围有显著区别，是因为与周围地区相比，该地区海拔高度更高，土壤温度更低，所以冻结开始时间更长，融化开始时间更短。
Figure 2. Spatial spatial distribution of soil freeze start time (a), Thaw START Time (b), and freeze duration (c) in the source region of the Yangtze River (The red rectangle indicates the range of the source region)
图2. 江河源区土壤冻结开始时间(a)、融化开始时间(b)和冻结持续时间(c)的空间分布(红色矩形表示江河源区的范围)
图3为1979~2018年江河源区土壤冻融循环关键参数冻结开始时间、融化开始时间、冻结持续时间的变化趋势。图3(a)和图3(d)中，冻结开始时间十年变化率大于0，延迟平均值为2.1天/十年，年变化趋势线呈上升趋势。图3(b)和图3(e)中，融化开始时间十年变化率小于0，延迟平均值为3.4天/十年，年变化趋势线呈下降趋势。图3(c)和图3(f)中，冻结持续时间十年变化率小于0，延迟平均值为2.5天/十年，年变化趋势线呈下降趋势。
总体来说，除极少数地区外，大多数网格的模拟数据得到的冻结开始时间显示统计上的延迟趋势，融化开始时间和冻结持续时间显示统计上的提前趋势。
3.3. 青藏高原土壤冻融循环关键参量与气温和海拔高度的相关性分析
图4和表4进一步探究对江河源区土壤冻结开始时间、融化开始时间和冻结持续时间与年平均气温和海拔进行相关分析。在图4中，土壤冻融循环关键参数与气温存在显著相关性，且相关性通过置信度99%的显著性检验。图4(a)年平均气温与冻结开始时间呈现明显的正相关关系，当年平均气温上升1℃，
Figure 3. Spatial distribution of the trends in soil freeze start time (a), thaw start time (b), and freeze duration (c), and annual variation of freeze start time (d), thaw start time (e), and freeze duration (f) in the source region of the Yangtze River (The black dashed lines indicate the linear trend lines)
图3. 江河源区土壤冻结开始时间(a)、融化开始时间(b)和冻结持续时间(c)的变化趋势空间分布和冻结开始时间(d)、融化开始时间(e)和冻结持续时间(f)的年变化(黑色虚线表示线性趋势线)
冻结开始时间延迟3.8天；图4(c)和图4(e)年平均气温与融化开始时间和冻结持续时间呈现明显的负相关关系，当年平均气温上升1℃，融化开始时间和冻结持续时间提前6天和9.3天。表4为土壤冻融过程关键参量与气温、海拔的相关系数，冻结开始时间与年平均气温相关系数为0.34，融化开始时间与年平均气温相关系数为0.96，冻结持续时间与年平均气温相关系数为0.83。
除年平均气温外，土壤冻融循环三个关键参数与海拔也进行了相关分析。冻结开始时间、融化开始时间和冻结持续时间与海拔相关性均通过置信度99%的显著性检验，表4中相关系数也较高。在江河源区，海拔升高，图4(b)中冻结开始时间提前，图4(d)中融化开始时间延迟，图4(f)中冻结持续时间推迟。但随着海拔升高，其余气象要素和非气象要素也可能会发生变化，要考虑是否有其他要素的变化对结果造成了影响。
图5为江河源区年平均气温与土壤冻融循环关键参量的相关系数的空间分布图。图5(a)江河源区范围内，年平均气温与冻结开始时间大范围呈现正相关，在江河源区西北部相关性性最明显。图5(b)年平均气温与融化开始时间在江河源区相关性系数偏小，相关性不明显。图5(c)冻结持续时间与年平均气温大范围呈现负相关，江河源区西北部冻结持续时间与年平均气温的负相关系数较大，相关性明显，东部相关性也较为明显，中部一般。冻结持续时间的变化是由冻结开始时间和融化开始时间变化共同作用的结果，冻结持续时间得到的相关系数图与冻结开始时间有一定类似[4]。冻结开始时间、融化开始时间、冻结持续时间与年平均气温相关系数在江河源区中部都不是很明显。
Figure 4. Correlation analysis between soil freeze start time (a), thaw start time (c), and freeze duration (e) with annual mean temperature, and correlation analysis between soil freeze start time (b), thaw start time (d), and freeze duration (f) with elevation in the source region of the Yangtze River
图4. 江河源区土壤冻结开始时间(a)、融化开始时间(c)和冻结持续时间(e)与年平均气温的相关分析以及江河源区土壤冻结开始时间(b)、融化开始时间(d)和冻结持续时间(f)与海拔的相关分析
Table 4. Correlation coefficients between key parameters of soil freeze-thaw process and temperature, elevation
表4. 土壤冻融过程关键参量与气温、海拔的相关系数
|
冻结开始时间 |
融化开始时间 |
冻结持续时间 |
| 年平均气温 |
0.34** |
−0.96** |
−0.83** |
| 海拔 |
−0.07** |
0.63** |
0.50** |
注：**表示计算的相关系数通过99%的信度检验。
3.4. 青藏高原土壤冻结持续时间变化趋势与海拔高度和纬度的关系
图6为近地表冻结持续时间变化趋势与海拔高度和纬度的关系。图6(a)中，江河源区冻结持续时间随海拔高度十年变化率，先增大后缩小，变化率在2500~3000 m达到一个峰值，每十年缩短6.73天，在3000米~4500米变化率较为稳定，缩短率在5.88天/十年附近，然后下降。根据图6(a)，认为近地表冻结持续时间变化趋势与海拔高度相关性一般。
Figure 5. Spatial distribution of correlation coefficients between annual mean temperature and soil freeze start time (a), thaw start time (b), and freeze duration (c) in the source region of the Yangtze River
图5. 江河源区年平均气温与土壤冻结开始时间(a)、融化开始时间(b)和冻结持续时间(c)的相关系数空间分布
图6(b)中，研究范围为北纬30度~37度，冻结持续时间十年缩短率30˚到31˚增加到5.89天/十年，31˚到33˚减小，33˚到34˚又增加6.25天/十年然后变化趋向较为稳定。根据图6(b)，近地表冻结持续时间变化趋势与纬度相关性一般。
(a) (b)
Figure 6. Relationship between trends in near-surface soil freeze duration and elevation (a) and latitude (b) (The numbers above the bar chart represent the number of simulated grid points in each elevation range)
图6. 近地表土壤冻结持续时间的变化趋势与海拔(a)和纬度(b)的关系(柱形图上方的数字是每个海拔范围的模拟网格数)
4. 结论
本文利用CLM5.0模拟资料的土壤温度与观测数据对比，分析了冻结开始时间、融化开始时间和冻结持续时间的时空变化特征和变化趋势，并研究了土壤冻融循环关键参量与气温、海拔的相关性。
1) CLM5.0模拟资料的土壤温度在江河源区适用性良好。
2) 江河源区冻结开始时间，集中在10月~11月份，冻结开始时间增长。融化开始时间大部分在3~4月，融化开始时间减短，冻结持续时间减少，呈现西北向东南变化的趋势。近四十年，江河源区呈现冻结开始时间延迟、融化开始时间提前、冻结持续时间提前的变化趋势。
3) 土壤冻融循环关键参数与气温和海拔存在相关性。1979~2018年，气温与冻结开始时间呈正相关，与融化开始时间和冻结持续时间呈负相关。在江河源区，海拔升高，冻结开始时间提前，融化开始时间延迟，冻结持续时间推迟。
4) 近地表冻结持续时间变化趋势随海拔高度先增加后减小，随纬度先增加再减小再增加然后趋向稳定，其相关性还需要更深一步的研究。