三江源地区季节冻土时空格局及影响因子

颜玉倩, 祁栋林,沈晓燕, 曹晓云,乔 斌,何彩青,赵梦凡,*

1 青海省气象科学研究所,西宁 810001 2 青海省防灾减灾重点实验室,西宁 810001 3 青海省气象台,西宁 810001

三江源地区位于青海省南部、青藏高原中部,作为青藏高原的腹地和主体,是长江、黄河和澜沧江三大河流的源头,面积达3.3×105km2[1]。其地理位置特殊、自然资源丰富、生态功能重要,是天气系统的上游区、气候变化的敏感区及生态环境的脆弱区,是我国生态环境安全和区域可持续发展的重要生态屏障[2—3]。三江源地区孕育着非常丰富的冻土资源,按地理位置划分属于青藏高原冻土区,该冻土区是世界中、低纬度地带海拔最高面积最大的冻土区,其范围北起昆仑山,南至喜马拉雅山,西抵国界,东缘横断山脉西部、巴颜喀拉山和阿尼马卿山东南部[4]。近年来,在全球气候变暖的大背景下,冻土的变化显得越来越重要,冻土的研究也越来越受到广大学者的关注。

季节冻土被定义为：冬天冻结而夏天融化的岩土层,它包括多年冻土区的活动层和非多年冻土区的土壤季节冻结层[5]。学者们利用遥感数据和数值模拟的方法估计了青藏高原季节冻土面积约为1.45×106km2,占高原面积的 56%—57.5%,其主要分布在34°N以南地区[6—7]。季节冻土最大冻土深度是一个既响应土壤又响应大气的独特指标[8],其冻融过程也影响着地表能量和水分交换[9]。大量研究表明[10—13]：近年来青藏高原气温升高、降水增多,向暖湿化发展,尤其温度的变化较大,其增温水平是全球平均水平的2倍(冬季高达0.3℃/10a—0.5℃/10a)。在青藏高原暖湿化的大背景下,季节冻土的变化主要呈现出最大冻土深度变浅、冻结日数缩短、冻结期缩短、融化期延长和范围萎缩等变化特征[14—17],这些变化无疑会对地下水循环、生态系统、岩土工程、基础建设乃至区域的可持续发展产生重要影响[18]。冻土变化造成的这些影响结果之间又相互作用,改变了区域的水文地质和水文条件,使植被逆向演替植物群落组成发生相当大的变化,而这种改变经常会导致草地生态系统结构和功能的改变,从而影响整个区域的生态过程[19]。此外,冻土退化使活动层和新形成的融区中的有机物质含量增加,使原来冻结在多年冻土中的碳暴露在地-气间碳循环过程,经微生物降解而释放温室气体到大气,从而使大气中的温室气体增加,进而使气候进一步变暖,形成这样一种正反馈的机制[20—21]。

以往的研究中,学者们多聚焦于多年冻土变化及分布,对于长时间序列、区域性的季节冻土的变化特征关注较少,并且在探讨研究季节冻土和气候因子时,只考虑了温度和降水或者单一因子的影响,对其他气候因子有所忽视[22—24],也未曾从大尺度天气背景场出发,探讨大气-冻土之间的相互关系。季节冻土的变化是一个复杂的过程,受局地因子的影响较大[25],本文从三江源地区1981—2020年40年季节冻土的时空分布特征出发,在温度、降水、湿润指数、≤0℃负积温、≤0℃负温日数、地表感热通量、地表潜热通量7个影响土壤、大气热力状况因子中寻找最能影响三江源季节冻土最大深度变化的气候因子,揭示土壤-大气互相影响的机制,为合理配置三江源地区资源进行农业生产、牧业发展和经济建设提供依据。

1 研究数据与研究方法

1.1 研究数据

三江源区的海拔高度和站点分布如图1所示,三江源地区位于青海省南部,青藏高原腹地,研究区内行政区域包括玉树、果洛、海南、黄南4个藏族自治州的16个县和格尔木市的唐古拉乡[3,26]。其中有气象站点23个,但由于位于三江源西部的沱沱河站(站号：56004)和五道梁站(站号：52908)无冻土数据,故采用21站的数据进行分析。由于气象观测中无潜热通量和感热通量的实际观测,尽管近年来随着野外站点的建立和观测设施的布设,已有一些实测数据可用,但仍然存在年限短、站点极少的问题,用再分析资料代替是一个比较折中的办法,已有大量研究指出[27—29],欧洲中心的再分析数据与其他数据相比,青藏高原地区的年际、年代际变化更接近台站资料,因此本研究采用欧洲中心ERA-Interim再分析数据。本研究所用的数据有以下几类：

图1 三江源地区海拔高度及站点分布示意图

(1)1981—2020年三江源区21个气象站季节冻土最大冻土深度(Maxmum Frozon Soil Depth,MFSD)(月值,单位：cm)；

(2)1981—2020年三江源区21个气象站温度(日值,单位：℃)；

(3)1981—2020年三江源区21个气象站20—20时降水(月值,单位：mm)；

(4)1981—2019年ERA-Interim地表潜热通量、地表感热通量再分析数据(月值,空间分辨率0.25°×0.25°,单位：W/m2)；

(5)1981—2019年ERA-Interim500hPa、100hPa位势高度(月值,空间分辨率1°×1°,单位：gadpm)300hPa温度(月值,1°×1°,单位：K)。

由于季节冻土最大冻土深度的变化是土壤和大气相互作用的体现,在选择气候因子时尽可能的考虑了表征热力状况的因子,主要有：温度、降水、地表感热通量、地表潜热通量、湿润指数、≤0℃负积温、≤0℃负温日数7个,其中前4个为单因子,后3个为复合因子。

(6)湿润指数

湿润指数是包含了蒸发(潜在蒸发)和降水在内的一个综合指数,潜在蒸发由 Thornthwaite方法求得[31]。湿润指数公式如下[31]：

(1)

式中:H为地表湿润指数;P为月降水总量；Pe为最大潜在蒸发量。月最大潜在蒸发量蒸发计算公式:

(2)

式中:Pe为可能蒸散量(mm/月)；Ti为月平均气温(℃)；I为年加热指数。

各月加热指数Ii由(3)式计算：

(3)

年加热指数I由公式(4)计算：

(4)

变量A由(5)式计算：

A=0.49+1.792×10-2H-7.71×10-5I2+6.75×10-7I3

(5)

负积温,即年内温度≤0℃温度的累加：

T=∑TkTk≤0℃

(6)

负温日数,即年内温度≤0℃日数的累加：

D=∑kTk≤0℃

(7)

1.2 研究方法

本研究采用的主要方法有以下几种：

(1)旋转经验正交函数分解(Rotated Empirical Orthogonal Function,REOF)

旋转经验正交分解是在经验正交分解(Empirical Orthogonal Function,EOF)的基础上展开的,EOF又称为主分量分析,是统计学上两种不同的提法,主要作用是将原变量的场进行时空分离,即分成空间的荷载向量(也称为特征向量,Load Vector,LV)和时间系数(也称为主分量,Principal Component,PC),构成为数不多的几个模态,每个典型模态都尽量多的反映原始场的信息,从而达到简单表达整个区域的变率结构,但EOF也有其局限：一是不能清晰的表示不同地理区域的特征,二是存在一定的取样误差,但REOF就可以克服以上EOF所带来的问题,旋转后的典型空间分布结构清晰,较好的反映不同地区的变化[32]。本研究采用REOF方法来揭示季节冻土最大冻土深度的时空分布特征。

(2)多元线性回归

多元线性回归指有2个或2个以上的自变量的线性回归方程,其原理是最小二乘法。本研究中因子选取的方法即是多元回归标准化系数的贡献率法,其因子的贡献率用R表示[33]：

(8)

其中ai为多元回归的标准化系数。

另外,本文采用相关分析法、合成分析法等常统计规方法探讨三江源地区季节冻土与气候因子之间的相互关系。

2 结果与分析

2.1 三江源地区季节冻土的时空分布特征

为分析季节冻土的空间分布特征,图2给出了1981年—2020年4年平均季节冻土MFSD的空间分布图,由图可以看出：季节冻土MFSD在200cm以上的站点有2个,占总站数的9.5%,MFSD在100cm—200cm的站点有13个,占总站数的61.9%,MFSD在100cm以下的站点仅有6个,占总站数的28.6%,即三江源区的大部分站点的最大冻土深度都在100cm以上,21站平均MFSD为136.66cm。黄河源区MFSD随海拔高度的减小递减,表现出较好的海拔高度特征,长江源和澜沧江源区更多的表现出纬度特征,即随纬度的较小,MFSD随之减小。

图2 1981—2020年三江源地区 季节冻土最大冻土深度(单位：cm) 季节冻土最大冻土深度变率示意图

从三江源地区季节冻土MFSD变率空间分布图中可以看出,三江源地区季节冻土MFSD近40年整体不同程度的减小趋势,即冻土厚度减薄,其递减率最高的站为曲麻莱站(站号：56021站),其递减率高达1.67cm/a。仅3个站出现最大冻土深度增加的情况,这3个站分别是：贵南站(站号：52955),玉树(站号：56029),河南(站号：56065),其中河南站有过迁站历史,对季节冻土深度的研究造成一定的影响,其结果需要进一步讨论。但这3个正变率值都较小。21站平均MFSD递减率为0.51cm/a,可以与之比较的是在已有的研究中,黄河源区从1961—2014年,季节冻土最大冻土深度的递减率为0.31cm/a,略微低于三江源地区最大冻土深度的递减率[17]。

对青藏高原站点季节冻土(标准化)进行旋转经验正交函数分解,由于前两个载荷向量(Rotated Load Vector,RLV1和RLV2)所占的方差贡献较大(解释方差分别为33.89%,10.02%),图3给出了REOF前两个载荷向量的空间分布特征,第一载荷向量在三江源地区的特征值均为负,表现出明显的全区一致型,代表整个三江源地区季节冻土变化的一致性的特点,这也是三江源季节冻土的主要分布型,由于已对原始数据做了标准化处理,其值大小反映了空间上季节冻土最大冻土深度的相对大小,重点反映出中部和北部的变化特征：即出现两个大值中心,分别是以三江源中部长江源中段和三江源东北部黄河源后段。结合第一载荷向量时间系数(Rotated Principal Component,RPC1),RPC1表现出明显的下降趋势,季节冻土深度逐年减小,这与上述图2所得出的结果也非常吻合：对三江源地区21站做1981—2020年平均,其与RPC1的相关系数高达0.89。并且时间系数的变化图中(图4)中可以看到,以2003年为界,2003年前(包括2003)RPC1为正值,2003年前(包括2003)RPC1为负值,可以解释为季节冻土最大冻土深度在2003前减薄,但最大冻土深度在其40年均值以上,2003年后继续减薄,在40年均值以下。减小的趋势在近4年表现的尤为明显。

图3 1981—2020年三江源地区季节冻土最大冻土深度REOF 第一模态、第二模态空间场

图4 1981—2020年三江源地区季节冻土最大冻土深度旋转正交分解第一模态 第二模态时间系数

第二载荷向量空间分布特征依旧为全区一致性,全区基本保持负位相,但在东部表现出小范围的正位相特征,整体呈现从南向北依次递减的特点,反映了南部的变化特征。第二特征向量时间系数RPC2虽然也表现出下降趋势,但是这种下降趋势并没有第一特征向量时间系数明显,其趋势系数只有-0.002,所以RPC2更多反映的是季节冻土在不同时间段内发生的变化,从RPC2的SG五点平滑函数来看RPC2经历了正-负-正-负-正交替变化,结合第二模态空间分布,在三江源地区季节冻土最大冻土深度表现出负-正-负-正-负的波动变化的特征。

2.2 三江源地区季节冻土气候因子分析

上述季节冻土的时空分布中可以看出三江源地区季节冻土层减薄,大量的研究表明冻土退化是由于下垫面的热力状况变化所引起的,这也是最主要的影响因素,但季节冻土的变化受到包括气候、植被、水文、人类活动等许多因素的影响[6,34],就单单气候因素而言就有：气温、降水、积雪、蒸发等等,这些因素间还有互相的影响及反馈,造成冻土变化研究的复杂性,本研究单从气候角度,寻找了7个表征热力状况的气候因子进行验证,这7个因子中有4个单因子：温度,降水,地表感热通量,地表潜热通量；3个复合因子：湿润指数,负积温,负温日数。湿润指数是一个既包含了温度又包含了降水,也考虑了蒸发又能反应土壤湿润程度的一个综合性指标,对研究冻土的变化有非常好的指示作用；考虑选择负积温和负温日数则是为了验证“冰冻三尺,非一日之寒”,是否季节冻土最大深度变化与之确实相关?

图5给出了各因子近40年空间分布,温度的空间分布特征与最大冻土深度的分布非常类似,温度的最低值点与最大冻土深度的极大值点位置相同,其空间分布也表现出从中部向四周扩散的逐渐升高的特点,三江源地区近40年平均气温1.83℃,温度最高站点尖扎站(站号：52963)与温度最低站清水河(站号：56034)平均温度相差12.63℃。表现出三江源地区气温分布的差异性,结合图1的高程分布图,可以明显的看出气温的差异主要受地形影响。三江源地区降水分布呈现出从西北向东南依次递增的特点,且三江源地区的降水量在整个青藏高原明显偏少[35]。最大降水量出现在三江源地区东南角久治站(站号：56067)。作为表征下垫面热力状况的重要因素[36],三江源地区地表感热通量和地表潜热通量的空间分布二者均表现出从东向西依次递减的空间分布型,不同的是,地表感热通量的最大值在三江源地区东北部,这与气温的的最高值位置相同,而地表潜热通量的最大值位于三江源地区的东南部,这与降水的最大值的位置是相同的。这是因为感热通量本身就是指由于温度变化而引起大气与下垫面发生的湍流形式的热交换,地表或大气在加热过程中,其相态没有发生变化,故温度和地表感热通量的关系密切,而地表潜热通量主要由水的相变产生,也被定义为大气和下垫面水分的热交换,因此降水和地表潜热通量的关系更为密切,这一点许多研究都有印证[3,27,29,37—39]。

图5 年平均温度、降水、地表感热通量、地表潜热通量、湿润指数、负积温、负温日数空间分布

三江源地区湿润指数的空间分布其最大值中心分别是清水河站和久治站,这里比较有意思的是,这两个站一个是温度的最低值中心,一个是降水的最大值中心。并且在黄河源的后段是湿润指数的低值区,表明该地气候干燥。负积温和负温日数者在空间分布及极大值中心的位置都非常类似,极大值中心位于中部的清水河站和南部的杂多(站号：56018)和囊谦(站号：56126)站,二者的相关系数高达0.98,即说明负温日数的增加就能造成负积温的增大。三江源区负积温和负温日数具有巨大的差异性,其最高值达-2209℃、218d,最低值低至-346℃、91d,造成巨大差异,三江源地区平均水平为-1076℃、156d,平均负温日数占全年的42.7%。

上述的7个与大气和土壤相关的热力因子中哪些是最能影响季节冻土最大冻土深度的因子呢？将MFSD作为因变量,其他7个因子作为自变量,通过比较回归方程的因子的标准化系数的贡献率及偏相关系数,确定影响最大冻土深的高影响因子。最终结果如表1所示：其中湿润指数的标准化系数贡献率最大,由大到小依次为：湿润指数>温度>降水>负积温>地表潜热通量>地表感热通量。排前三位的因子系数的贡献值比较大,三者累积可达到86.23%,且偏相关系数绝对值的大小也可进一步证明,故认为湿润指数、温度、降水是影响季节冻土最大冻土深度的主要因素。负温日数和负积温之所以剔除,可以解释为：最大冻结深度取决于进入土壤的年度能量收支,所以可能导致单独负积温关系不好。因而很明显看出,温度才是影响季节冻土的最主要的因素。将地表感热通量和地表潜热通量剔除是因为二者不MFSD的直接影响因素,在上文中也提到,地表感热通量是先与温度相互影响、地表潜热通量是先与降水相互影响,因而热通量的变化对MFSD具有一定的滞后性,因而被剔除。

表1 多元回归标准化系数及贡献率及偏相关系数

图6给出了最大冻土深度第一模态时间系数与三个所选择出来的最能影响三江源地区季节冻土最大冻土深度的因子的相关系数分布场,由图可得：RPC1与温度的相关系数分布呈现出东部和西部较大,中部地区较小的分布特征,但三江源地区均通过了α=0.1的显著性检验,表明整个三江源区MFSD与温度的相关关系都非常好。RPC1与降水的相关系分布中基本是负相关,东北部为通过α=0.1的显著性检验的区域。RPC1与湿润指数的相关系数场的分布于降水的类似,但在值上略有差异,通过α=0.1的显著性检验的区域比较少,只有三个站点包含其中。将三个因子相关系数分布通过显著性检验的区域叠加,得到高影响因子的的关键区,但在叠加的过程中,湿润指数通过显著性检验的站点较少,为后续寻找典型的高低值年份带来不利的影响,故将其范围调整至大于-0.15的区域,最终得到如图关键区,关键区中包含7个站点。

图6 RPC1与 温度、降水、湿润指数相关系数空间分布气候因子影响冻土的关键区

对关键区7个站点MFSD做平均处理后得到1981—2020年的时间序列,并对其做标准化处理(图7)。以±1个标准差(σ)为依据,挑选出关键区季节冻土的典型高值年和典型低值年,其结果如表2所示。典型高值年份有：1983、1984、1986、1993；典型低值年份有：1988、1994、1999、2010、2013、2017、2019。并且从图中可以看到,1986年到1988年最大冻土深度下降剧烈(两年间下降了近30cm),从前文的分析中可知,气温和降水是影响季节冻土最大冻土深度最主要的因素,且与最大冻土深度呈负相关,因而在1986—1988 年,温度、降水、湿润度三者均呈现上升趋势,在温湿协同作用下这两年的最大冻土深度呈现出明显的下降趋势。

图7 关键区7站季节冻土最大冻土深度、温度、降水和湿润指数时间变化特征

表2 季节冻土最大冻土深度典型高值年、低值年

2.3 大尺度天气背景对季节冻土的影响机制

从以上分析中可以看出,温度和降水是影响季节冻土的最直接因素,那么大尺度天气背景又是如何影响到温度和降水的分布的呢？本研究将从极涡和南亚高压两方面出发进行讨论。

极涡是影响我国乃至全球天气气候的最主要环流实体之一,也是冷空气活动的最主要标志,且该系统在500hPa等压面上最为强盛[40]。因此图8给出了典型高值年和典型低值年北半球500hPa位势高度场平均场和距平场的合成图。由图可以看出在典型高(低)值年,其距平场表现为负(正)距平,对应季节冻土更厚(薄)。关于极涡对我国温度和降水影响的相关研究很多,张恒德等[41]研究指出,极涡面积的大小与我国温度呈显著负相关,且与春夏秋冬四季相比较,全年的平均温度与极涡的相关性最好,即当极涡面积较大(小)时,温度较低(高),与我们之前的结论相对应。也有研究指出[42],北半球极涡面积指数与青海降水具有较好的负相关关系,即当极涡面积指数较大(小)时,青海省降水偏少(多),对应季节冻土较厚(薄)。

图8 典型高值年和典型低值年500hPa位势高度场平均场、距平场(单位：gadpm)合成

南亚高压是长期活动在青藏高原上空最稳定的最强大的高压系统,并与青藏高原发生强烈的陆气相互作用[43]。其活动对北半球大气环流具有重要作用。大量研究表明南亚高压的强度和位置对我国的旱涝分布有着显著影响。图9给出了季节冻土最大冻土深度典型高值年和典型低值年100hPa高度场和300hPa温度场的平均场和距平场合成。从图中可以看出,二者的南亚高压的中心强度均为1660dagpm,但在典型高值年,其范围极小,南亚高压的主体位势高度以1655dagpm为主；而典型低值年其范围较大,几乎包含了南亚的大部分地区,包括青藏高原在内。除此之外从温度场的合成中可以发现：在高层300hPa,南亚高压中心与温度的大值中心相配合,在典型高(低)值年,温度场的大值中心温度更低(高),分别是-32℃和-30℃,这与之前的关键区与温度的相关性分析的研究也比较吻合,即温度越低(高),季节冻土就越厚(薄)。其距平场也显示出,在最大冻土深度典型高值年,南亚高压表现为负异常,同时对应300hPa的温度更低,青藏高原地区温度负异常值在-0.5℃至-0.6℃；在最大冻土深度典型低值年,南亚高压表现为明显的正异常,并且与300hPa温度场的正异常中心相吻合。李跃清[44]指出,南亚高压的负(正)异常对应青藏高原东侧的干旱(洪涝)年,马振锋[45]等与之有较为类似的结论。本研究也充分说明了该结论：在典型高(低)值年,南亚高压偏弱(强),高原东侧三江源地区偏于干旱,降水偏少,而在典型区域的相关分析中我们知道关键区最大冻土深度与降水呈现反相关的关系,即降水偏少,季节冻土较厚。

图9 典型高值年和典型低值年100hPa位势高度场、300hPa温度场平均场及距平场合成

3 结论

本文利用三江源地区21个气象台站的最大冻土深度、气温、降水资料,欧洲中心再分析资料对三江源地区季节冻土MFSD时空分布特征进行了分析,并探讨了季节冻土与气候因子之间的相关关系,得到以下结论：

(1)1981—2020年三江源地区21个气象站40年平均MFSD空间分布呈现出以玛多站和清水河为最大值中心,向四周递减的分布特征,中心最大值为218.85cm,21站平均MFSD深度为136.66cm。40年平均变率空间分布和旋转正交经验分解第一模态均表现出近40年三江源地区季节冻土MFSD呈现明显的下降趋势,冻土层明显减薄,21站平均最大冻土深度递减率为0.51cm/a。

(2)将最大冻土深度作为因变量,温度,降水,地表感热通量,地表潜热通量,湿润指数,负积温,负温日数7个表征热力状况的气候因子作为自变量,建立多元线性回归模型,通过多元回归方程的标准化系数的贡献率,选择出最能影响季节冻土的气候因子,分别是：气温、降水和湿润指数。

(3)通过分析高相关因子与最大冻土深度的相关系数分布场,找到三江源区气候因子影响季节冻土的关键区：三江源东北部地区,包括7个气象站点。并找出三江源地区季节冻土最大冻土深度的典型高值年分别是1983年、1984年、1986年、1993年和典型低值年分别是1988年、1994年、1999年、2010年、2013年、2017年、2019年。

(4)通过对500hPa位势高度场典型高值年典型低值年合成分析,季节冻土典型高(低)值年,500hPa位势高度场负(正)异常,对应三江源地区应季节冻土更厚(薄),南亚高压负(正)异常,南亚高压范围偏小(大),强度偏弱(强),温度场的大值中心温度更低(高),对应三江源地区应季节冻土更厚(薄)。