气候温升背景下青藏工程走廊带多年冻土热融蚀敏感性分布预测研究

崔福庆， 刘志云， 张 伟， 王 伟， 陈建兵

（1.长安大学地质工程与测绘学院，西安 710054； 2.中国交通建设股份有限公司第一公路勘察设计研究院有限公司高寒高海拔地区道路工程安全与健康国家重点实验室，西安 710065）

随着青藏工程走廊带内一系列重大线性工程的建设落成，以及青藏高速列入《国家高速公路网规划》和《国家公路网规划（2013—2030）》，如何保障寒区重大工程安全性和稳定性、提高工程构筑物服役性能、降低冻融循环和冻土退化导致的工程地质灾害，是青藏高原多年冻土研究和工程建设中面临的巨大挑战. 在人为活动或工程影响下，将会使多年冻土由天然上限处向下加速融化. 地温条件、地表条件和土壤含冰状态不同将导致冻土对外界热量扰动的敏感性响应过程不同，从而使得季节融化深度产生差异，一般将其定义为冻土热融蚀敏感性，通常用季节融化深度与潜在季节融化深度的比值来表示冻土热融蚀敏感性［1］：

式中：STe为冻土热融蚀敏感性；X为季节融化深度（m）；Xp为潜在季节融化深度（m）.

多年冻土的空间分布特征是工程建筑物的安全稳定性重要影响因子［2-4］，其热融蚀敏感性可视为冻土环境对于外部传热过程的响应速率，是描述多年冻土特征的关键指标之一. 因此，开展青藏工程走廊带内多年冻土热融蚀敏感性的分布特征研究，既是冻土环境监测与保护的前提，也是指导寒区工程构筑物设计、建造及后期养护科学决策所必需的先决条件之一［5-8］.

多年冻土与外部环境之间有着复杂的非稳态能量交换过程，同时与大气环境、地表植被覆盖情况、海拔、太阳辐射等多个因素直接关联，因此青藏工程走廊带内多年冻土热融蚀敏感性的分布特征具有很强的地域性和特殊性［9-11］. 吴青柏等提出冻土热融蚀敏感性和冻土年平均地温具有很强的相关性，但同时也和冻土的季节融化深度（活动层厚度）相关，所以如何获取二者的分布是确定热融蚀敏感性分布的关键. 早期多年冻土年平均地温大区域分布特征研究主要采用野外布设传感器监测方法或现场人工钻探方法［12］，但青藏高原幅员广阔，加之高寒缺氧的恶劣气候条件，使得走廊带全线的地温监测或钻探代价高昂无法实现. 近年来，随着遥感技术的不断发展和研究手段的增加，以及相关数据获取途径的增多，基于数字高程的局部因素对冻土分布的定量化影响研究获得了快速发展［13-14］. 与此同时，基于GSM-R和GPRS通信系统的冻土地温自动检测系统也陆续投入使用，从而实现了冻土地温变化的实时监控与掌握［12，15-16］. 在多年冻土活动层研究方面，张中琼等［17］运用基于融化指数的Stefan 公式，计算和预测了不同温升速率下（0.04 ℃/a、0.02 ℃/a和0.017 ℃/a）青藏高原冻土活动层厚度的变化特征；Chaves等［18］对南极洲的Keller半岛的表层地温进行了监测，获取了气候温升条件下活动层厚度和冻土热流的年际变化特征，证明了极地冻土热收支与大气热收支的强烈正相关联系；Wu 等［19］对2002—2012年北麓河地区10个野外观测点的地温数据进行了分析，结果表明该地区冻土活动层厚度平均增加速度为4.26 cm/a，其主要诱因为夏季降雨量的增加；张明礼等［20］开发了基于水分和能量平衡的冻土水-汽-热耦合数值模型，分析了降雨增加和温升对冻土活动层厚度和热稳定性的影响；Du等［21］依据现场实测结果，提出地表浅层土介电常数可用于预测所有类型冻土的活动层厚度并给出了幂次拟合公式；李韧等［22］根据唐古拉观测场太阳辐射和冻土地温监测数据，分析了地表能量收支的季节变化规律以及其与冻土活动层厚度的非线性变化过程.

本文依据青藏工程走廊带冻土年平均地温和活动层厚度的野外监测数据，结合MODIS 和SRTM-DEM遥感数据，建立了多年冻土热融蚀敏感性预测模型，并采用开放系统地气耦合模型对未来气候温升模式下多年冻土年平均地温和活动层厚度变化进行数值研究，建立二者与气候温升的相关关系，从而对青藏工程走廊带未来多年冻土热融蚀敏感性分布特征进行预测.

1 热融蚀敏感性预测模型

1.1 冻土年平均地温拟合公式

冻土地温受地表植被覆盖情况、高程、地表温度、降雨和大气热收支等多个因素影响. 本文依据野外监测数据，以上述影响因素为变量，建立了青藏工程走廊带冻土年平均地温的多元回归拟合公式［23］. 在进行地温计算时，首先采用多年冻土辨识Logistic模型对走廊带多年冻土分布进行辨识，冻土分布的概率计算式如下：

式中：P为多年冻土存在的概率；w0为拟合公式常数项；wi为与冻土分布有关的环境、地形等因子的量化指标；αi为回归方程系数.

利用公式（2）就可以建立多年冻土存在与否与外部环境之间的联系. 公式（2）中各参数及其回归系数参见文献［23］. 表1是采用多年冻土Logistic辨识模型与实测结果对比情况，表中切割值为0.7，即概率P大于0.7判断为冻土，否则判断为融土. 由表1可以看到，此模型对于非多年冻土和多年冻土的预测准确率分别达到88.2%和93.3%，基本满足辨识精度要求.

表1 多年冻土Logistic辨识模型与实测结果对比Tab.1 Comparison between prediction of permafrost logistic identification model and field observation

而后，对走廊带内127个地温监测点的年平均地温观测数据和地温相关影响因子进行分析拟合，得到青藏工程走廊带冻土年平均地温多元线性回归模型拟合公式：

式中：φ为纬度；H为高程；φ′为等效纬度；N为植被指数大于0 的算术平均值. 上述对应位置数据可通过NASA发布的2000—2016年MODIS遥感数据和SRTM-DEM数字地形数据产品获得.

1.2 冻土活动层厚度拟合公式

大量研究表明［19-21］，多年冻土活动层厚度与冻土年平均地温密切相关，且受地理区域、地表性状、坡度等因素影响. 本研究基于青藏工程走廊带内多年冻土活动层厚度的371个实测数据及对应位置的地温和遥感数据，建立了青藏工程走廊带多年冻土活动层厚度多元线性回归模型：

式中：h为冻土活动层厚度；B为拟合公式常数项；xi为与冻土分布有关的环境、地形等因子的量化指标；βi为回归方程系数. 不同地温条件下的具体数值参见表2.

表2 活动层厚度拟合公式系数Tab.2 Coefficients of the active layer depth fitting formula

1.3 冻土热融蚀敏感性多元线性回归模型

热融蚀敏感性反映了冻土自身对于外部气候环境和工程构筑物水-热-力影响的响应速率. 已有研究表明，热融蚀敏感性与多年冻土年平均地温和季节融化深度等参数均有一定联系［1］. 表3为26组冻土年平均地温和活动层厚度实测数据与冻土热融蚀敏感性的相关性分析结果.

表3 冻土年平均地温、活动层厚度与热融蚀敏感性参数相关性Tab.3 Partial correlation analysis among the permafrost thermal thawing sensitivity，the mean annual ground temperature and active layer depth

由表3结果可知，冻土年平均地温和活动层厚度两参数与冻土热融蚀敏感性显著相关，对其进行多元线性回归，得到冻土热融蚀敏感性多元线性回归模型如下：

图1为26个实测点基于公式（5）所得的热融蚀敏感性值与基于Kudryavtsev方法计算的热融蚀敏感性数值的对比情况. 由图1可知，本研究所提出的预测模型计算结果与文献［1］理论公式方法计算值具有较好的一致性. 预测模型的R2=0.935，说明了本模型准确性与可靠性.

图1 冻土热融蚀敏感性预测模型计算值与文献［1］理论公式计算值比较Fig.1 Comparison between prediction values of permafrost thermal thawing sensibility and calculation results of empirical formula in Ref［1］

2 多年冻土年平均地温与活动层厚度数值预测模型

2.1 数值计算模型

通过建立统一的地气耦合数值模型，研究气候温升背景下多年冻土年平均地温和活动层厚度随时间的变化规律. 计算模型由空气环境和天然地层构成（图2）. 为使计算结果更具一般性，天然地层简化为单一土性，同时考虑到青藏工程走廊带表层土的分布特征，天然地层土性依据青藏公路沿线野外勘察结果取表层分布最为广泛的砾石土［24］，天然土层的厚度取20 m. 为减少模型上部壁面边界对于空气区湍流计算的影响，空气层的高度取20 m. 此外，为避免入口效应对计算结果的影响，计算模型长度取120 m.

气象观测资料表明，青藏高原5 m高度年平均风速约为5.2 m/s［25］. 因而对于冻土地层内外部的耦合传热过程采用二维非稳态、湍流模型进行数值求解. 在数值计算中，冻土地层上部空气环境视为自由流体，并将空气考虑为密度不变的不可压缩气体，采用标准κ-ε模型进行湍流计算［26］：

式中：κ为湍流脉动动能；ηt为湍流脉动所造成的动力黏度；σκ为脉动动能的Pr数；Gκ为湍流脉动动能产生项；ε为湍流耗散率；σε为湍流耗散的Pr数；c1和c2为经验常数.

冻土地层的主要传热方式为热传导，本文采用显热容法处理相变热对路基传热的影响，其控制方程为［27］：

式中：C*和λ*分别为各层固体材料的等效体积热容和等效导热系数.

图2 数值计算模型示意图Fig.2 Schematic of numerical model

2.2 边界条件与计算参数

为反映一年之中风向变化对于冻土地层耦合换热的影响，模型空气区域左右边界根据计算时刻的不同，分别取为速度入口边界条件和充分发展的自由出流边界条件，风速v的取值按照其距地面高度的不同，设为与地面高度呈指数变化的函数形式：

模型的上边界取温度边界类型，壁面温度Tenv设为时间周期性变化函数：

式中：T0为年平均气温；g(t)为年地温增幅；A为气温振幅；φ0为初始相位. 冻土地层左右两侧取绝热边界条件. 为模拟地热对于冻土层的热效应，冻土地层下表面设0.06 W/m2的热流边界条件［28］. 同时，为了体现出土层上表面对太阳辐射的吸收、地表与天空背景的长波辐射换热及地表水分蒸发导致的散热等过程对于冻土地层温度场的影响，故将地层表面设为厚度0.1 m的耦合内热源项的换热壁面. 通过对各类换热过程分析可知，地表耦合换热壁面源项值QS的计算式为：

式中：α为吸收系数；Qsun为投射到地表的总太阳辐射；Qr为地表对环境的长波辐射热损失；Qe为地表通过水分蒸发带走的热量. 上述各项以及风速和气温的取值或计算方法具体见文献［10］和［24］.

采用FLUENT6.3.26软件对气温增幅为0.022 ℃/a条件下青藏高原多年冻土年平均气温和活动层厚度变化情况进行了数值模拟. 数值计算基于二维非稳态湍流模型，并采用隐式算法进行数值求解. 空气的定性温度取0 ℃，冻土地层的物性取温度相关的阶梯函数. 为体现外部环境周期性变化对于定图地层耦合换热的影响，环境温度、风速、太阳辐射、蒸发量等均设为随时间年度周期性变化函数. 计算时间设为50 a，计算残差的收敛标准统一取10-5.

3 计算结果与分析

3.1 冻土年平均地温与活动层厚度变化规律

基于所建立的地气耦合数值模型，对气候温升条件典型类型冻土地层的年平均地温与活动层厚度进行了模拟. 图3 为年平均气温-0.5～-5.5 ℃、年度温增0.022 ℃/a 工况条件下，冻土年平均地温未来50 a 的变化情况. 由图3可以看到，开始阶段（第0年），年平均地温与年平均气温基本相当，以后随着气温的逐年增加，冻土地层的年平均地温也随之增加. 且高温冻土（T0＞-1 ℃）受气温年际增加的影响较小，低温冻土（T0≤-1 ℃）受气温年际增幅的影响相对较大. 年平均气温-5.5 ℃工况下，其年平均地温增加0.77 ℃.

图3 不同年平均气温条件下冻土地层年平均地温逐年变化情况Fig.3 Variations of annual average ground temperature of permafrost foundation under different annual average air temperature conditions

表4 为不同年平均气温条件下冻土活动层厚度在2036 年和2066 年的变化情况. 由表可知，随着气温年际上升（0.022 ℃/a），多年冻土活动层厚度也随之增加. 在年平均气温-1 ℃工况下，2036年的活动层厚度相对于2016 年将增加0.54 m，2066 年将增加1.35 m；而在-5.5 ℃工况下，2036 年和2066 年的冻土活动层厚度将仅分别增加0.09 m和0.43 m. 这说明随着冻土地温的增加，冻土地层对于外界环境热扰动的响应也在加强，外部换热对于冻土地层的影响深度随之增大. 此外，还可以看到高温冻土的活动层厚度在2036年或2066年的增加量将远大于低温冻土，表明低温冻土的热融蚀敏感性要小于高温冻土.

表4 不同年平均气温条件下冻土活动层厚度的变化Tab.4 Variation of active layer depth under different ground temperature 单位：m

3.2 青藏工程走廊带冻土年平均地温分布预测

根据冻土年平均地温的数值计算结果，可求得2036年和2066 年冻土年平均地温增幅与年平均气温的线性拟合公式：

式中：a和b分别为线性拟合公式系数，2036年预计为-0.133 9和-0.147 2，2066年预计为-0.023 7和-0.253 7.

将此拟合公式代入公式（5），并结合遥感MODIS 和SRTM-DEM 数据产品以及多年冻土Logistic 辨识模型［23］，进行栅格运算后，可得到未来20 a和50 a青藏工程走廊带内多年冻土年平均地温分布（图4）. 由图4可知，青藏工程走廊带内融区和高温冻土区主要存在于河流、谷地和盆地等区域，且随着气温的逐年增加，低温冻土（-7～-1.5 ℃）的区域逐渐变小，中高温冻土区域（-1.5～0 ℃）和融区增大：至2036年时，低温冻土、中高温冻土和融区的比例将分别为45.6%、30.4%和24%；至2066年时，低温冻土面积比例将由当前的58.9%降低为28.2%，中高温冻土比例将由当前的30.8%上升为47.9%，说明受全球气候变暖的影响，青藏工程走廊带内多年冻土将逐年退化.

图4 未来20年和50年后青藏工程走廊带内冻土年平均地温分布Fig.4 Annual average ground temperature distribution of Qinghai-Tibet engineering corridor in the next 20 and 50 years

3.3 多年冻土热融蚀敏感性分布特征

同理，基于冻土活动层厚度的数值计算结果，亦可求得2036年和2066年冻土活动层厚度变化量与年平均气温的多项式拟合公式. 在此基础上，综合MODIS和SRTM-DEM遥感数据产品、多年冻土Logistic辨识模型，利用公式（3）和（4）计算每个图像栅格未来20 a和50 a的年平均地温和冻土活动层厚度值，然后利用多年冻土热融蚀敏感性预测模型进行栅格运算，可求得青藏工程走廊带多年冻土热融蚀敏感性分布预测图.

依据文献对于热融蚀敏感性的划分标准［29］，分别将所得结果依次分为不敏感型、弱敏感型、敏感型和极敏感型4类. 对于敏感系数区间分别为≤0.54、0.54～0.66、0.66～0.8和＞0.8，结果如图5所示. 由图5可知，随着气候逐步变暖，昆仑山、可可西里、风火山和唐古拉山等地热融蚀敏感性为不敏感型冻土将大幅退化为弱敏感型冻土，而楚玛尔河平原、北麓河盆地、乌丽盆地、通天河、沱沱河等地的敏感型冻土也将急剧退化为大片极敏感型冻土. 各类型冻土在2036 年和2066 年的比例变化如表5 所示. 可以看到，随着气候不断变暖，青藏工程走廊带内冻土的热融蚀敏感型冻土比例将大幅增加. 至2036年时，不敏感型、弱敏感型、敏感型和极敏感型4类冻土比例将分别为14.8%、14.4%、25.2%和19.2%，而至2066年时，极敏感型冻土比例将增长1 倍以上，由当前的16.5%上升到36.7%，且敏感型和极敏感型冻土将占整个走廊带内多年冻土区的78%以上，而不敏感型冻土将下降为不足当前的一半. 还可以看到，走廊带内极敏感型冻土的增加比例将随时间而加速增长，相应的非敏感型冻土将加速蜕变为敏感型冻土. 可以认为，随着环境温度不断增加，人类工程活动的不断加剧，青藏工程走廊带内多年冻土将变得更为脆弱.

图5 青藏工程走廊带多年冻土热融蚀敏感性分布预测Fig.5 Distribution prediction of permafrost thermal thawing sensibility in QTEC

表5 青藏工程走廊带多年冻土热融蚀敏感性分布比例预测Tab.5 Prediction of the proportion distribution of permafrost thermal thawing sensibility in QTEC 单位：%

4 结论

青藏工程走廊带多年冻土热融蚀敏感性分布特征研究，可为走廊带内冻土环境监测与保护、寒区工程构筑物设计和建造提供科学依据. 本文通过对现有冻土年平均地温和活动层厚度野外监测数据的分析，结合MODIS和SRTM-DEM遥感数据，建立了多年冻土热融蚀敏感性预测模型，并采用开放系统地气耦合模型对气候温升模式下多年冻土年平均地温和活动层厚度变化进行数值研究，给出了青藏工程走廊带内多年冻土区热融蚀敏感性分布预测图. 研究结果表明：

1）走廊带内冻土年平均地温越低，其上升幅度随气候温升越大，而活动层厚度将随地温和气温的升高而增大，在年平均气温-1 ℃和-5.5 ℃工况下，年平均地温和活动层厚度增幅分别为0.023 ℃/1.35 m 和0.77 ℃/0.43 m.

2）融区和高温冻土区主要分布在青藏工程走廊带沿线的河流、谷地和盆地等区域，且随着气温的逐年增加，至2066年低温冻土区域比例将减少52.1%，高温冻土区域和融区面积比例总计增大74.7%.

3）随着气候温升，走廊带内多年冻土的热融蚀敏感性将大幅增加，且极敏感型冻土的增加比例将随时间而加速增长，非敏感型冻土将加速蜕变为敏感型冻土. 至2066年时，极敏感型冻土比例将由当前的16.5%上升到36.7%，且敏感型和极敏感型冻土将占整个走廊带内多年冻土区的78%以上.