高原寒区地表土壤冻融作用的空间参数化表征研究

高会然， 许 冲， 张万昌， 易亚宁， 肖子亢

（1. 应急管理部 国家自然灾害防治研究院，北京 100085； 2. 中国科学院 空天信息创新研究院，北京 100094；3. 复合链生自然灾害动力学应急管理部重点实验室，北京 100085）

0 引言

在全球气候持续变暖的背景下，冰冻圈领域已成为研究热点［1-4］。当前，冰冻圈研究不仅关注其在气候系统中的作用，还更加关注对人类社会可持续发展的现实影响和潜在威胁［5-7］。作为冰冻圈的重要组成部分，冻土的变化或退化对生态环境系统、人类生活和生产安全等方面的直接或间接风险，这些研究方向也越来越受到学者的关注［8-9］。

当前，冻土在全球范围内正呈区域性退化趋势，主要表现为多年冻土面积缩小、活动层加厚，季节冻土层变薄、冻结时长缩短等［10-13］。与多年冻土变化相比，季节冻土冻融循环过程周期短、范围广，对陆地表面与大气之间的物质与能量交换、陆地表面景观格局的影响剧烈而深远［5-6，14-16］。现有研究表明，季节冻土对地表景观的影响机制主要体现在冻土冻融的风化作用、冻胀和融沉作用以及冻结滞水作用［17］。一方面，季节冻土是深度不断变化的隔水层，并通过影响区域土壤水特性而形成特殊的水文过程和生态过程［18-19］。另一方面，在地表季节冻土冻融过程中，山体边坡土壤结构、坡体的稳定性、空隙水压力、土体摩擦力等物理特性发生不同程度的变化，为浅层地质灾害的发育提供致灾条件和驱动力［20-23］。因此，摸清地表季节冻土的分布格局及其时空变化对开展寒区自然科学研究、保障或维护人类的生产活动和生存环境安全均具有重要的研究意义。

传统的冻土研究多采用地面观测的方法，但是传统方法在数据实时获取、信息快速提取以及特征参量空间离散化的研究需求下面临重大挑战［24-26］。20世纪80年代以来，国内外学者利用遥感技术进行了大量的区域或全球尺度的冻土时空动态遥感监测研究，并取得一系列数据或技术成果［27-32］。但是，遥感冻土监测受限于当前遥感技术水平，低空间分辨率的遥感冻土数据产品在中小区域尺度上的实际应用仍十分有限。随着冻土物理学的发展以及对冻土系统中的水热传输特性和机制的深入研究，国内外学者发展出了多种基于物理机理的冻土水热传输过程数值模型。例如，Flerchinger 等［32］构建的一维水热耦合过程模型SHAW（Simultaneous Heat and Water）模型，Wood等［33］提出的适用于大区域尺度水文过程模拟的VIC（Variable Infiltration Capacity）模型，陈仁升等［34］建立的嵌入分布式水文过程模型的高寒山区水热耦合模型（DWHC），Gao等［31］基于水热耦合原理研发了一套空间全分布式冻土水热过程数值模型（Fully Distributed Frozen Soil Hydrothermal Processes Integrated Modeling System，FFIMS 模型），以及近十年来提出或改进的多种考虑冻融循环的陆面过程模型，如GEOTOP 模型、CLM 模型、UW-VIC 模型等［35-38］，为冻土系统动态过程研究提供了有效的方法和手段。

当前，基于物理机理的冻土水热传输过程数值模型研究取得了较大的进展，部分陆面过程模型通过模型集成的方式，具备了精细刻画地表土壤冻融循环多参量和多过程的能力［39-40］。但是，目前冻土变化监测与模拟研究成果在陆面环境要素响应和冻土灾害防治等领域的应用能力仍处于较低的水平。一方面，高时空分辨率冻土特征参量数据获取难度相对较大，限制了大范围长时段冻土冻融循环及其冻融作用的定量分析；另一方面，地表土壤冻融作用作为冻土水热过程与陆面景观格局演变相互联系的重要纽带，其表征指标识别与空间参数化的定量研究不足。区域尺度冻土冻融作用的参数化表征及其进一步在灾害风险评估与防控领域中的应用研究仍是空白。

因此，为了精细刻画全球气候变化下的季节冻土水热过程，系统性表达冻土冻融循环的时空变化特征以及冻融作用特征，本研究以地处青藏高原东南缘的高山峡谷区及其周边地区为研究区，综合考虑气象、植被、积雪与土壤等多个陆面过程以及各个过程中的水热动态平衡，在基于冻土水热耦合系统性理论的空间全分布式的冻土水热耦合过程数值模型FFIMS 模型框架下，进行了适用于青藏高原高海拔冻土区的冻土水热过程数值建模。在研究区开展长时段（2010年8月1日至2020年7月31日，共10 个水文年）逐日冻土水热过程模拟，获取研究时段内高精度的冻土系统空间分布式水热过程参量，并着重分析研究区冻土环境特征及其时空变化特征。在此基础上，结合冻土边坡稳定性理论，创新性提出冻土冻融作用的空间参数化方法，定量研究地表土壤冻融循环对土体结构稳定性的影响，为区域尺度冻土相关研究和冻融灾害风险评估与灾害防控等研究提供思路借鉴和数据及技术支撑。

1 理论与方法

1.1 总体方法框架

空间分布式的季节冻土特征参数是开展地表土壤冻融作用空间参数化的数据基础，本研究首先利用现有的基于物理机制的分布式冻土水热过程数值模型，获取研究区长时间序列的冻土特征参数数据集，包括植被冠层、积雪、土壤剖面等冻土系统多要素、多过程空间水热参量。然后，利用GIS空间分析、数据挖掘等数据处理技术，分析近十年青藏高原东南缘季节冻土冻融循环及其变化特征，揭示气候变化下的复杂地理环境中的冻土时空演化规律。最后，根据土壤剖面水热过程的多种冻土特征参量，结合土体抗剪强度和冻土边坡稳定性相关理论，提出可以表征冻土冻融作用的空间参数化方案。主要研究框架如图1所示。

图1 本研究的主要方法框架示意图Fig. 1 Overall framework of the integrated approach adopted in this study

1.2 分布式冻土水热过程数值建模

作为冻土过程数值模型基础理论的水热耦合原理考虑了气象、植被冠层、积雪与土壤等多个过程以及各个过程中的水热平衡［31，34，41-42］，各个过程中的水热平衡方程见式（1）～（5）。

（1） 植被冠层

植被及其冠层设计为单层单节点结构，同时考虑了植被生长和季节变化特征，冠层能量与水平衡方程如式（1）和式（2）所示。

式中：ρa、ca和T分别为大气密度（kg·m-3）、大气比热容（J·kg-1·℃-1）和冠层温度（℃）；ke为冠层热量传导系数（m2·s-1）；Hl为冠层向大气传输的热量（W·m-3）；ρv为蒸气密度（kg·m-3）；El为冠层蒸散发量（mm）。

（2） 积雪

将积雪设计为单层双节点的结构，考虑积雪的累积和升华以及水热在雪层中的传导过程，积雪能量平衡方程可用式（3）表示。

式中：ρsp、wsp和ksp分别为雪密度（kg·m-3）、体积含水量（m3·m-3）和热传导系数（W·m-1·℃-1）；ci和ρl分别为冰的比热容（J·kg-1·℃-1）和水密度（kg·m-3）；Rn为积雪下的净辐射通量（W·m-2）；Lf和Ls分别为融化潜热与升华潜热（J·kg-1）；qv为水汽通量（kg·s-1·m-2）。

（3） 土壤

设计为不同网格不同深度的多层土壤结构，考虑了冻融过程对水热通量的影响，土壤能量和水平衡方程可分别表示为式（4）和式（5）。

式中：Cs为土壤比热容（J·kg-1·℃-1）；ks为土壤热传导系数（W·m-1·℃-1）；K为土壤导水率（m·s-1）；ψ为土壤水势（m）；U为土壤水通量的源/汇项（m3·m-3·s-1）。

在FFIMS 模型的水热平衡方程的求解过程中，将植被冠层、积雪和土壤剖面分为有限层数，每层用一个节点表示，每个节点中的各水热分量的存储量以各节点所在层的土壤厚度为准。利用隐式有限差分方程组求解每个节点以及节点间的能量和水的平衡方程［43］。模型中各个子过程的理论和方法涉及大量参数和变量，包括一些可调参数、经验常数和物理常数等。其中，可调参数主要是指受研究区域环境条件影响的参数，需要根据研究区地形、气象、生态等条件调整。本研究建立的模型中主要的参数取值及其说明见表1。

表1 FFIMS模型可调参数及其在本研究中的取值Table 1 Variable parameters of FFIMS model and their values in this study

FFIMS 模型的驱动数据、输入和输出均为空间分布式的数据，本研究中所有的数据均统一采样为1 km 空间分辨率，模型模拟的时间步长为24 h（逐日模拟）。其中，模型驱动数据和输入数据主要包括气象驱动数据（气温、降水、风速、相对湿度、大气压强、日照时长等要素）、土壤属性数据（土壤类型、土壤深度、土壤机械组成、有机质含量等属性）、地表覆被属性数据（叶面积指数、叶面特征、植被高度、干生物量等）、土壤温湿度初始状态数据。模型输出数据主要是可以表征冻土系统特征或状态的变量，包括冠层水热参量（冠层辐射、冠层截留等）、积雪水热参量（雪面辐射、积雪厚度、雪水当量等）、土壤剖面水热参量（包括土壤剖面温度、土壤剖面含水/冰量、地表热通量等）。

1.3 冻融作用参数化表征方法

地表土壤的冻融作用直接影响或引起土体分裂以及土壤颗粒的重新组合，从而影响土体的抗剪强度，导致土体从稳定状态转向相对不稳定状态，威胁地区的边坡稳定、工程安全等［44-45］。但是，冻融作用对土抗剪强度以及相关指标的影响程度目前仍无定论。根据本研究获取的精细冻土特征参量，结合现有的理论研究基础，可以实现定量化表征冻土冻融作用对土体抗剪强度指标的相对损伤程度。因此，本研究以土壤剖面含水量和土壤剖面含冰量两个关键冻土特征参量为主要参数，构建土体抗剪强度的损伤系数。

土壤内摩擦角和黏聚力是表征土壤抗剪强度的两个关键指标。已有研究表明，土体抗剪强度随着土壤含水率的增加而降低，其中土壤内摩擦角随土壤含水率增加而减小，黏聚力随之多呈现先增大后减小的趋势［46-47］。土壤在冻结状态下，土体中既有固体冰也有液态水，此时的土体抗剪强度主要由土体固有特性和固体冰以及两者间的相互作用决定的［48-50］。一般情况下，土壤含冰量越高，土体抗剪强度越大。因此，本研究以统计时段内的土壤剖面含冰量的变化量为依据，将土体抗剪强度损伤系数表示为e的指数函数，见式（6）。

式中：f为土体抗剪强度损伤系数（0≤f＜1），f值越小，表示冻融作用对土体抗剪强度的损伤越大；Δθice为某时段内土壤含冰量的变化量（m3·m-3），计算方法见式（7）；k为关于土壤含水量的系数，计算方法见式（8）。

式中：D为统计时间段（d）；θice，d为d时刻的土壤剖面含冰量（m3·m-3）；θw，d为d时刻的土壤剖面含水量（m3·m-3）；θ'D为D时段土壤剖面含水量的平均值（m3·m-3）。

式（6）～（8）表示，随着冻土融化，土壤含冰量降低，土壤含水量上升，土壤含冰量的变化量增加，f值随之下降，意味着冻融作用对土体稳定性的损伤加大。土壤含水量通过影响土壤黏聚力从而对土体抗剪强度产生影响，关于土壤含水量的调节系数k决定了f值变化的程度。当土壤开始融化至融化初期，在土壤含水量低于θ'D时，土壤含水量缓慢增加，k值不断增大，f值下降缓慢甚至有所增加，此时土体抗剪强度由土壤水和固体冰以及两者间的相互作用等条件决定；当土壤含水量高于θ'D时，k值开始减小，f值降低，土壤冻融作用对土体抗剪强度的损伤增大。另外，由于土壤抗剪强度随冻融循环次数的增加呈先减小后趋于稳定的趋势，因此，多个水文年或多个时段计算的f值不能直接累加计算。

2 研究区与数据

2.1 研究区概况

以青藏高原东南部的高山峡谷地区为研究区（90.67°～104.2° E、28.6°～31.46° N），总面积为4.12×105km2（图2）。研究区地势起伏剧烈，自东向西分别为四川盆地西部、青藏高原东南缘高山峡谷地带以及青藏高原腹地南部，海拔范围为302～6 436 m（平均约4 000 m）。研究区中部大面积的高山峡谷区历史构造运动强烈，地质构造发育，且岩土体破碎严重，是各类地质灾害的高易发和高风险区，对当地生态环境安全和人类生产建设安全具有潜在威胁。

图2 研究区位置和基本地理环境Fig. 2 Location and basic geographical environment of the study area： location （a）， terrain and distribution of observation stations （b）， soil types （c）， and land use/cover types （d）

21 世纪初，研究区多年冻土和季节冻土覆盖全域，其中金沙江以西主要为大片岛状多年冻土，以东为季节冻土和短时冻土。受气候变暖的影响，青藏高原总体变暖的趋势十分明显［51-52］。近40 年来，青藏高原升温的速率比全球同期升温速率高2 倍左右［2］，青藏高原东南缘是青藏高原地区升温最显著的区域之一［53］。因而，该地区多年冻土和季节冻土都发生了剧烈的变化和退化，对高原地表景观形态和地表自然过程产生了深远的影响。在此背景下，探索和查明高原寒区季节冻土的分布格局及其时空变化，是为寒区基础自然科学研究和保障人类生产活动和生存环境安全提供科学依据的重要手段。

2.2 基础数据

本研究收集了DEM 高程数据、土壤质地、土地利用类型等基础资料，其中，DEM 来源于USDS EROS 数据中心（https：//srtm. csi. cgiar. org/srtmdata/）发布的90 m 空间分辨率的SRTM DEM。土壤类型空间分布数据来源于全国土壤普查办公室编制并出版的《1∶100 万中华人民共和国土壤图》，基本单元为土壤亚类（研究区共有79 类，土壤深度约0.4～1.5 m），包含16 种土壤属性数据项，如土壤深度、土壤机械组成、土壤有机质含量等，土壤属性数据库可在线访问（http：//vdb3. soil. csdb. cn）。土地利用类型数据来源于中国多时相土地利用现状数据库。土壤类型和土地利用类型数据均下载自中国科学院资源环境科学数据中心（https：//www.resdc.cn）。FFIMS模型通过数据查找表的方式实现模型所需的土壤和植被相关参数的空间参数化。

2.3 气象驱动数据

（1）气温和降水

空间分布式的气温和降水数据是FFIMS 模型的主要驱动参数，本研究收集了研究时段内的1 km分辨率气温和降水同化数据产品（a high-resolution and long-term gridded dataset for temperature and precipitation across China， HRLT）［54］，数据来源于地球和科学环境科学数据库（https：//doi. pangaea.de）。HRLT 数据产品使用综合统计分析方法，结合高程及其地形参数、地形湿度指数等协变量，对中国国家气象数据共享网发布的0.5°网格逐日气象数据参数进行空间插值，得到全国长时段高分辨率的最高气温、最低气温和降水数据产品。最高温度的平均绝对误差（MAE）、确定系数（R²）和Nash建模效率（NSE）分别为1.07 °C、0.98 和0.98；最低温度的MAE、R²和NSE 分别为1.08 °C、0.99 和0.99；降水数据的MAE、R²和NSE 分别为1.30 mm、0.71和0.70。

（2）其他气象参数

其他气象参数包括相对湿度、潜在蒸散发、大气压强、风速、日照时长、0 cm 地表温度等，数据主要来源于中国国家气象数据共享网提供的研究区范围内的25 个气象站（图2）的逐日观测资料（http：//data.cma.cn）。采用自然邻域空间插值方法（Natural Neighbor， NN）［55］，将气象站点数据转换为1 km 分辨率的空间格网数据。由于上述气象参数具有空间大尺度的特点，利用有限的站点数据插值的方法是可行的。

2.4 站点观测数据

对于模型精度验证，本研究获取了位于青藏高原中东部那曲地区附近的土壤温/湿度监测数据集（SMTMN），该监测网络利用ECH2O EC-TM 电磁传感器获得了研究区56个站点（位于本研究区的站点为16个，图2）上的地表土壤湿度和温度的长时间逐日地面观测数据集，数据采集时段为2010 年8 月1 日至2016 年12 月31 日。该数据集通过利用土壤质地和土壤有机碳含量对土壤含水量和土壤温度观测记录进行了校准和质量控制处理，具有较高的精度［56］。

2.5 遥感数据产品

本研究收集了两种地表土壤水分遥感数据产品，与FFIMS 模拟的土壤含水量进行对比分析。一个是风云三号B星携带的微波湿度计获取的逐日地表土壤湿度数据集（以下简称“FY-3B”），下载自风云卫星数据中心（http：//satellite. nsmc. org. cn），该数据集的空间分辨率为0.25°，时间范围为2011年1月1 日至2018 年12 月31 日。另一个是融合了AMSR-E（NASA National Snow and Ice Data Center Distributed Active Archive Center）、AMSR2（Global Change Observation Mission）和SMOS（Soil Moisture and Ocean Salinity）土壤湿度数据产品的逐月地表土壤湿度数据集（以下简称“AAS”），下载自国家青藏高原科学数据中心（https：//data. tpdc. ac. cn），该数据集的空间分辨率为0.05°，时间范围为2011年1月至2018年12月。

2.6 土壤温湿度初始条件

由于模型的初始土壤水热条件是多层的结构，而土壤温湿度数据通常为地表单层土壤的状态。因此，在模型输入数据预处理模块，预设了两种土壤剖面一维温度场和湿度场的计算方法。通过地表单层土壤的温湿度数据，即可近似估算土壤垂直剖面上的温度和水分分布。土壤剖面一维温度场和湿度场的计算方法分别见式（9）和式（10）［5，57］。

式中：Tz为土壤深度为z时的土壤温度（℃）；z0为地表土壤层的厚度（m），z0

3 结果与讨论

3.1 冻土水热过程数值模拟结果与验证

FFIMS 模型模拟的部分结果展示见图3，空间分布图包括土壤剖面平均温度、土壤剖面平均含水/冰量以及地表热通量，曲线图为逐日的区域平均土壤温度、含水/冰量。

图3 FFIMS模型部分模拟结果的时空分布Fig. 3 Spatiotemporal distribution of partial simulation results of the FFIMS model

陆面过程模型模拟的参量较多，而用于模型精度验证的实测数据类型相对较少。在模型模拟过程中，各参数之间相互联系又相互影响。一般而言，如果主要模拟参数的验证指标良好，则通常认为模型精度是可靠的。土壤剖面温度和含水量是冻土的关键特征变量，本研究主要利用青藏高原那曲地区16个野外站点的土壤温湿度观测数据，采用纳什系数（NSE）、均方根误差（RMSE）和决定系数（R2）指标，验证FFIMS 模型的精度，总体精度验证结果见图4，逐年精度验证结果见表2。FFIMS 模型模拟的土壤剖面温度的纳什系数和均方根误差分别为0.89（0.87～0.91）、2.38（2.33～2.74），土壤剖面含水量的纳什系数和均方根误差分别为0.59（0.44～0.74）、0.05（0.04～0.06），该结果验证了FFIMS 模型在刻画冻土水热过程中具有良好的精度。从R2指标的验证结果看，模型模拟的地表土壤温度和土壤含水量的时间变化趋势在验证时段内是十分一致的，证明了FFIMS 模型对冻土水热过程模拟的稳定性。

表2 FFIMS模型逐年精度验证结果Table 2 Annual precision verification results of the FFIMS model

图4 FFIMS模型在16个冻土观测站的精度验证结果：土壤温度模拟结果与实测值时间序列对比（a），土壤温度模拟精度（b），土壤温度模拟值与实测值的箱线图（c），土壤含水量模拟结果与实测值时间序列对比（d），土壤含水量模拟精度（e），土壤含水量模拟值与实测值的箱线图（f）Fig. 4 Overall accuracy verification results of the FFIMS model： comparison of time series between simulated soil temperature results and measured values （a）， simulation accuracy of soil temperature （b）， boxplot of simulated and observed soil temperatures （c）， comparison of time series between simulated soil water content results and measured values （d），simulation accuracy of soil water content （e）， and boxplot of simulated and observed soil water contents （f）

图5展示了FFIMS模型模拟的土壤表层水含量与现有数据产品在逐日和逐月尺度上的对比验证结果。可以发现，模型模拟土壤水含量的变化趋势与现有的土壤湿度遥感数据产品具有较好的一致性。但是，与逐日的遥感数据产品的剧烈变化幅度相比，模型模拟的地表土壤水含量表现出了更稳定的变化趋势。

图5 FFIMS模型模拟的土壤表层水含量与遥感数据产品对比验证结果Fig. 5 Comparison of surface soil water content simulated by the FFIMS model and remote sensing data products：daily FY-3B soil moisture data （a）， and monthly AAS soil moisture data （b）

3.2 冻土变化特征分析

（1）时间变化特征

本研究主要通过土壤剖面平均温度、土壤剖面平均含水/冰量和地表热通量等参数描述冻土系统水热过程，图6 展示了以上多种参量的逐日和逐年变化曲线，图中灰色区域的上下边界分别表示该参量在研究区范围内的最大和最小值，红色折线表示年均值。其中，地表热通量指地表土壤单位面积的热量交换（W·m-2），具有方向性，正值表示大气向土壤传递热量。

图6 冻土水热参量逐日和逐年时间变化趋势Fig. 6 Trend of daily and annual changes in the water and heat parameters of frozen soil： profile temperature （a），liquid water content （b）， ice content （c）， and surface heat flux （d）

在本研究时段内，研究区气温增幅明显，年平均气温自4.0 ℃上升到5.9 ℃，平均增幅为0.2 ℃·a-1。总的来说，从逐日的冻融过程看，地表土壤冻融循环中的各水热参量呈周期性的规律变化，逐年变化幅度不明显。从研究区平均的角度看，近十年，青藏高原东南缘冻土系统呈现明显的消退趋势，主要表现为土壤剖面温度上升、土壤剖面含水量增加等。其中，土壤剖面的年均温度从3.5 ℃上升到6.2 ℃，平均温度上升率达1.9 ℃·（10a）-1；土壤剖面含水量从0.2 m3·m-3上升到0.26 m3·m-3，平均上升率为0.07 m3·m-3·（10a）-1；土壤剖面含冰量呈现波动降低趋势，但变化幅度较低，总体趋势平稳；地表热通量则呈显著上升趋势，从2.9 W·m-2上升到12.0 W·m-2，平均上升率为6.8 W·m-2·（10a）-1。地表土壤热通量上升，说明大气向陆地输送的热量增加，且热通量数值为正，意味着土壤温度总体上升的趋势仍将持续。

（2）趋势率空间分布特征

空间倾向率方法是在普通线性倾向率方法的基础上，扩展至空间网格尺度上的方法，即在每一个空间栅格上建立相关变量的时间序列，对时间序列逐一进行线性拟合，获取线性拟合的斜率（趋势率）与显著性检验值，从而形成空间分布的趋势特征。本研究利用空间倾向率方法分析研究区冻土相关特征参数的时空演化规律，重点分析各水热参量变化的空间分布特征。根据逐年土壤剖面平均温度、含水/冰量、地表热通量空间分布数据进行空间倾向率计算，获得其时间变化趋势的空间分布，结果如图7所示。

从冻土水热参量变化趋势率的空间分布上看，各种参量的空间分布具有显著的异质性，空间分布特征受地形地貌、土壤质地等地表景观因素以及气候要素的影响。由于土壤温度是其他参量变化的主导因素，土壤温度与土壤含水/冰量等要素的趋势率在空间上的分布特征具有一定程度的相似性。其中，土壤温度变化趋势的空间分区特征最为明显，研究区西部地处青藏高原腹地的区域，除南部低洼地区，土壤温度总体呈显著上升的趋势；中部川藏高山峡谷地区地形多变，土壤温度变化的空间特征复杂；东部四川盆地的土壤温度呈不显著下降趋势或基本无变化，但与四川盆地交界的青藏高原东南边缘地带，是研究区内土壤温度上升幅度最大的区域。土壤含水量总体上呈现上升趋势，研究区大部分区域为显著上升趋势，局部为不显著上升，研究区南部洼地小面积的区域存在下降趋势。土壤含冰量总体呈现不显著下降的变化趋势，在空间分布上与土壤含水量上升趋势的分布特征较一致，其中呈显著下降趋势的地区多位于研究区东部的地势较低的沟谷地带。随着年均气温的升高，研究区地表土壤热通量呈逐年增加的趋势，且大部分地区较为显著，局部地区呈现下降的趋势。由于地表热通量随时间变化较剧烈，即使在同一天也可能出现较大的波动，因此其变化趋势的空间分布无明显规律性特征。

3.3 土壤冻融作用表征指数

根据FFIMS 模型模拟的精细冻土水热过程参量，按照本研究提出的冻融作用参数化表征方法，计算了研究区土体抗剪强度损伤系数。图8展示了不同冻融循环次数（1 个水文年为1 次冻融循环，图中以n表示）下的土体抗剪强度损伤系数空间分布。

图8 不同冻融循环次数下的土体抗剪强度损伤系数空间分布Fig. 8 Spatial distribution of the damage coefficient of soil shear strength under different freeze-thaw cycles

基于冻土冻融循环水热参量计算的土体抗剪强度的损伤系数，主要反映了区域地表土壤冻融作用对土体抗剪强度的相对损伤程度。损伤系数越低，表示损伤程度越大，边坡失稳发生冻融灾害的风险或概率（相对）加大。从时间变化上看，随着冻融循环次数的增加，冻融作用的损伤程度不断加剧，但其空间分布特征总体较稳定。从空间分布上看，土体抗剪强度损伤系数具有较强的地带性分布规律。高值区域集中分布在研究区东部和南部低洼峡谷地区，这些地区冻结天数较少甚至未发生冻结，土壤冻融作用相对较弱。无论在研究区西部的青藏高原腹地、中部的川藏高山峡谷地区，还是研究区东部低洼地区，均存在损伤程度较高（低值）的区域。其中，青藏高原东部边缘与四川盆地接壤地带有成片低值区域分布。受青藏高原东部边缘多年冻土退化程度较大的影响，该区域地表土壤冻融作用整体较强。

本研究以冻融循环5 次的计算结果为例，进一步结合地形、植被类型进行了统计分析（图9）。统计结果表明，地表土壤冻融作用对土体的损伤程度与高程有较显著的相关性，高程低于1 000 m 的区域损伤程度基本较小，而青藏高原腹地和高山峡谷区高程在5 000 m 左右的区域，土体强度受冻融作用影响最大。由于研究区内的人类活动多处于地势低洼的地区，因此土地利用类型为农田和建设用地的区域，地表土壤冻融作用的损伤程度较小。反之，高海拔山区的草地和林地，尤其是地表裸露的区域（盐碱地、戈壁、裸土/石等），土体抗剪强度在冻融循环的作用下损伤程度较大。另外，通过对比不同植被覆盖率下的损伤系数发现，植被覆盖度越高，损伤系数越高，土壤冻融作用的损伤程度越低。

图10 展示了研究区平均土体抗剪强度和土壤冰融化量随冻融循环次数增加的变化。其中，年均损伤系数和累积损伤系数的区别在于，计算累积损伤系数的基准年份是研究时段的起始年份。通过定量统计可以发现，累积损伤系数随冻融循环次数的增加而降低，但随冻融循环次数的增加，累积损伤系数降低的程度逐渐减小，直至平稳的趋势。这表明冻融循环导致土体损伤失稳的作用具有一定的限度，不会成为引起冻融灾害发生的主导因素。

图10 土体抗剪强度损伤系数和土壤冰融量随时间的变化Fig. 10 Variations of the damage coefficient of soil shear strength and melted soil ice with time

另外，需要特别注意的是，本研究计算得到的损伤系数并不能反映土体抗剪强度本身大小或边坡固有稳定性特点。例如，损伤系数较低的区域，表示冻融作用对土体的损伤程度大，但并非意味着边坡稳定性差或冻融灾害发生概率高，仅表示该区域地表岩土体受到的冻土冻融作用更强烈。除了冻融作用外，土体本身的抗剪强度对边坡稳定性起到决定性作用，冻融循环起到了弱化其相对强度的作用。

4 结论

针对当前分布式冻土过程模拟研究的不足，弥补区域尺度冻土研究在自然灾害风险评估与防控领域中应用研究的短板，本研究以青藏高原东南缘的高山峡谷区及其周边地区为研究区，重点关注季节冻土冻融循环及其水热传输过程的系统性表达，建立了适用于青藏高原高海拔冻土区的空间全分布式的冻土水热耦合过程数值模型（FFIMS 模型）。在获取高精度的冻土系统分布式水热过程参量的基础上，提出冻土冻融作用的空间参数化方法，得到以下主要结论：

（1）基于物理机制的FFIMS 模型对地表冻融循环中的冻土水热过程的模拟效果较好，冻土系统主要特征参量（土壤温度和土壤含水量）经过验证精度良好，可以为精细刻画空间分布式冻土过程提供有效的模型方法。

（2）随着气温升高，青藏高原东南缘季节冻土变化剧烈，空间异质性较强，季节冻土除了周期性冻融循环外，总体呈退化的趋势，表现为土壤温度显著上升、土壤含水量显著上升、土壤含冰量不显著下降等，为冻融作用下的岩土体结构的抗剪强度变化增加了更多不确定性因素。

（3）另外，还进行了地表土壤冻融作用空间参数化表征方法研究，创新性地提出了土体抗剪强度损伤系数，通过统计分析表明该指标在数值表达和空间分布等方面是合理可行的。利用高精度的冻土过程特征参量，可以有效进行冻土冻融作用的参数化表达。

本研究成果可以为冻土相关研究提供新的思路，同时，为寒区区域尺度下的相关基础研究和灾害风险评估与灾害防治等领域的研究提供冻土系统动态演化的数据和技术支撑。

致谢：本文使用的土壤温湿度观测数据和遥感反演的土壤水分数据下载自国家青藏高原科学数据中心（http://data.tpdc.ac.cn），土壤类型和土地利用类型数据下载自中国科学院资源环境科学数据中心（https://www.resdc.cn），气象驱动数据下载自地球和环境科学数据库（https://doi.pangaea.de）以及国家气象数据共享网（http://data.cma.cn），在此一并表示感谢。