青藏高原平均海拔在4 000 m以上， 是世界上海拔最高、 地形最为复杂的高原， 被称作地球“第三极”^［1］。由于青藏高原特殊的地理位置、 高大地形的动力作用和地面热源的强烈热力作用， 使其形成了独特的大气环流模式， 对我国、 甚至对全球气候都产生了深远的影响^［2］。青藏高原地区广泛发育着多年冻土， 面积约为1.05×10⁶ km^2［3］， 与北极地区的多年冻土相比， 青藏高原的多年冻土具有地温高、 活动层厚度大及热稳定性差等特点^［4］。作为全球气候变化的指示器与放大器， 青藏高原对气候变化十分敏感^［5］。在全球变暖背景下， 多年冻土活动层厚度及热状况的变化将在一定程度上对区域的水文循环产生影响， 且活动层厚度的增加与土壤温度的升高也将导致土壤蒸发量增大， 从而改变云量， 进一步影响整个区域的能量平衡^［6］。活动层与多年冻土的热状况控制着多年冻土区的生态过程^［7］， 随着气候变暖， 多年冻土退化将加速青藏高原土壤有机碳的分解， 并释放温室气体到大气， 增强气候变暖趋势， 形成对气候变化的正反馈效应^［8］， 活动层内部的土壤温湿度分布也将发生变化， 导致区域植被退化^［9］。多年冻土区工程建设中常见的路基融沉和冻胀问题也与活动层内部的热状况密切相关。因此， 研究青藏高原活动层土壤热状况的动态变化对深入了解高原土壤的热力作用及其对气候变化的响应、 生态环境的影响、 高原地区的工程建设与维护等方面都有重要意义^［10］。多年冻土区地-气系统之间的能量和水分交换通过位于多年冻土之上的活动层来实现^［11］。土壤的冻融深度是陆面与边界层大气水热交换的客观表现， 是下垫面水热综合作用的结果， 在某种程度上， 活动层厚度也是下垫面热状况的重要特征^［2，11］。土壤温度的变化是下垫面热量储存和释放的重要体现， 也是下垫面热状况变化的重要标志^［10］。有研究表明， 青藏高原下垫面热量的储放使东亚天气系统的槽脊配置发生了变化， 最终触发天气系统的变化^［4］。目前， 针对青藏高原土壤热状况的研究主要分为两类： （1）通过实测数据分析土壤热状况的变化特征； （2）利用模型对土壤温度进行模拟， 以分析或预测活动层热状况对气候变化的响应。多年冻土区气候条件恶劣， 交通不便， 野外调查十分有限。模型通过统计分析或数学方程来描述各种物理过程， 为研究活动层土壤热状况的变化提供了有效方法。

FEFLOW（Finite Element subsurface FLOW system）是由德国WASY（Wasserwirtschaftsplanung und Systemforschung）公司开发的一款地下水数值模拟软件， 灵活的有限单元剖分方法和强大的数值模拟算法使其能对流场、 溶质运移和热量传输等多种情景进行合理的模拟与预测^［12-13］。FEFLOW 7.1中加入了piFreeze模块， 本构关系描述了热容的温度依赖性， 包括相变潜热、 渗透率和流体黏度以及导热系数， 考虑了热量输送过程中温度变化导致的水分相变以及伴随的土壤水热性质的变化， 可用于模拟冻融过程中土壤温度的动态变化^［14］。基于改进的FEFLOW， Anbergen等^［15］模拟了冻结过程中土壤温度的变化过程， 其完美拟合了McKenzie等^［16］给出的Lunardini经典解析解。Magnin等^［14］模拟预测了西欧勃朗峰地区小冰期末期以来至21世纪末的多年冻土退化情况。Langford^［17］模拟了加拿大西北地区斯科蒂河流域的活动层及下覆多年冻土的温度变化情况， 模拟效果良好。下文所提到的FEFLOW模型均为加入piFreeze模块后的FEFLOW模型。

目前， 在对活动层和多年冻土的温度进行模拟时， 常将较易获得的气温作为模型的上边界条件。Luo等^［18］在不考虑地表特征的热力影响情况下， 分别采用气温、 陆面温度、 地表温度作为GIPL模型的驱动变量， 对黄河源区的土壤温度进行了模拟。结果表明， 地表温度驱动模拟效果最为理想， 更适合作为土壤温度模拟的驱动变量。在多年冻土区， 地表温度的变化受到海拔以及包括植被、 土壤质地、 太阳辐射等在内的表面特征因素的影响， 具有复杂性和不确定性^［19］。BP（back-propagation）神经网络作为信息处理的一种新方法， 具有自适应性、 学习能力和极强的非线性处理能力^［20］。目前， 已经有诸多学者将BP神经网络模型用于青藏高原冻土区的相关研究， 包括冻土层上水水位^［21］、 冻土路基变形^［22］、 土壤冻结温度与未冻水含量^［23］、 径流的变化预测^［24］。本研究将利用BP神经网络模型对研究区的地表温度进行模拟预测， 为FEFLOW模拟和预测活动层不同深度土壤温度做好准备。

有研究表明， 过去几十年中青藏高原多年冻土活动层厚度不断增加， 土壤温度也逐渐升高^［5，25］。1981 - 2010年， 青藏高原活动层厚度变化量为-1.54 ~ 2.24 m， 平均以1.29 cm·a^-1的速率逐渐增厚^［25］。多年冻土退化将对青藏高原地区水热循环乃至整个生态系统产生影响， 模拟和预测青藏高原多年冻土活动层的热状况对气候变化的响应是非常必要的。为此， 本研究在已有野外观测数据的基础上， 利用BP神经网络模型和FEFLOW模型对青藏高原风火山研究区活动层各深度土壤温度进行模拟重现， 并检验两模型的有效性； 模拟预测不同升温情景下活动层各深度土壤温度的动态变化过程， 分析气候变化对不同深度土壤冻融特征及活动层厚度的影响， 这将为多年冻土区气候变化、 地-气系统之间的物质和能量交换、 生态学和工程建设的研究提供有用的信息。

如图1所示， 研究区位于青藏高原风火山左冒西孔曲小流域内（93°03′ ~ 92°50′E、 34°40′ ~ 34°48′N）， 左冒西孔曲流域为长江源区北麓河的重要支流， 该区域属于典型的连续多年冻土。在本研究中， 活动层温度监测点位于左冒孔河左岸的典型高寒草甸区域， 该区域多年冻土发育良好。研究区海拔为4 680 ~ 5 360 m， 属高原大陆性干旱气候， 年平均气温为-5.2 ℃， 年平均降水为290.9 mm， 但年平均蒸发量高达1 316.9 mm， 年平均相对湿度为57%， 多年冻土厚度为50 ~ 120 m， 活动层厚度变化范围为0.8 ~ 2.5 m^［21］。该流域内无冰川与常年积雪分布， 河流以大气降水和积雪融水以及地下水补给为主。植被以高寒草甸为主， 河谷和山坡底部有高寒沼泽分布。高寒草甸植被群落以高山嵩草（Kobresia pygmaea）、 矮嵩草（K. humilis）和线叶嵩草（K. capillifolia）等为主^［21］。

采用冻土工程国家重点实验室研发的热敏电阻传感器以及美国Fluke公司生产的多用表（Fluke180系列）， 计算地表及地表以下5 cm、 10 cm、 15 cm、 20 cm、 40 cm、 80 cm、 160 cm、 320 cm处的土壤温度。传感器的测量范围为-40 ~ 50 ℃， 测量精度为0.05 ℃。2006年1月1日至2008年12月31日， 每3 h记录一次活动层各深度的温度， 将日平均观测值作为一天的温度数据进行计算， 其中320 cm处的土壤温度自2006年9月2日开始观测并记录。在试验场地附近建立了一个微型气象站（美国ONSET公司HOBO）， 用以记录降水、 气温、 风速和风向以及总辐射数据。各层土壤的容重、 粒径组成见表1。

BP神经网络是一种按误差逆向传播算法训练的多层前馈网络， 由输入层、 一个或多个隐藏层、 输出层及各层中相互连接的节点组成（图2）。在训练学习的过程中， BP神经网络通过反向传播误差信号来不断调整网络的权值和阈值， 使网络实际输出与期望输出的均方根误差最小， 是目前研究最多、 应用最为广泛的人工神经网络模型之一^［20］。在本研究中， 基于BP神经网络建立地表温度预测模型的过程主要包括以下步骤： 样本构建、 网络训练、 验证、 测试。

神经网络建模的首要之处在于输入变量的选取， 它提供了被模拟系统的基本信息^{［21，24］}。考虑到多年冻土区数据匮乏， 采用统计学方法确定神经网络模型中对于地表温度预测具有显著影响的输入变量的滞后数^［21］， 以缩短试错过程， 提高模拟精度。本文中采用偏自相关函数（partial autocorrelation function， PACF）确定对于地表温度预测具有显著影响的输入参数的滞后天数， 地表温度时间序列的偏自相关系数见图3。由图可知， 当滞后天数为7天时， PACF仍显示出显著的相关性， 此后控制在置信区间内。通过对比发现， 滞后7天输入的模拟精度与滞后3天输入的模拟精度相比， 并无显著提高， 最终选取当天及向前三天的数据作为模型的输入， 分别代表模拟地表温度当天及前一天、 前两天与前三天的输入参数， 用t、 t-1、 t-2、 t-3表示。

图4展示了气温与地表温度的变化过程， 两变量的变化趋势几乎完全一致， 地表温度随气温的升高和降低而升高和降低， 相关系数高达0.94（P<0.01）。因此， 可以将气温作为影响地表温度的输入变量。图5展示了总辐射与地表温度的变化过程， 两变量的变化趋势基本一致， 相关系数为0.61（P<0.01）。因此， 将总辐射也作为影响地表温度的输入变量。分别用A_t 、 A_t-1、 A_t-2、 A_t-3与R_t 、 R_t-1、 R_t-2、 R_t-3表示输入的当天、 前一天、 前两天与前三天的气温与总辐射。理论上， 神经网络模型输入参数的种类越多， 越有利于提高模拟的精度^{［21，24］}。考虑到多年冻土区数据获取较为困难， 且有时增加输入参数反而会降低模拟的精度^［21］， 本研究中仅选取气温与总辐射作为模型的输入参数。

在对BP神经网络模型进行训练和验证前， 首先要确定模型的拓扑结构。本研究中隐藏层的数目为1层， BP神经网络的拓扑结构为8-1-1。参考下式^［22］， 确定隐藏层的节点数。

式中： L为隐藏层节点数； m为输出层节点数； n为输入层节点数； a为0 ~ 10之间的常数。

输入层中的8个节点分别为A_t 、 A_t-1、 A_t-2、 A_t-3、 R_t 、 R_t-1、 R_t-2、 R_t-3， 输出层的1个节点即表示所模拟的地表温度ST_t。利用Matlab R2017a中的newff函数创建BP神经网络， 将训练样本带入网络中进行训练， 建立地表温度预报模型。模型的训练期为2006年1月1日至2007年12月31日， 验证期为2008年1月1日至2008年12月31日。表2给出了训练期气温、 总辐射及地表温度数据的统计参数， 其中气温和地表温度数据呈较低程度负偏态分布， 总辐射数据呈较低程度正偏态分布， 训练期此三组数据离散程度较小， 代表性好。

FEFLOW模型中热量传输计算基于能量守恒方程^［26］。

式中： ε为有效孔隙度； ρ_lc_l、 ρ_sc_s分别为液相和固相物质的容积热容量； q_i 为水流渗流速度； x_i 、 x_j 为空间变量； λ_i，j 为导热系数； T为温度； Q_T为热量的源或汇项。FEFLOW中只考虑一种液相和固相物质， 即液态水与土壤颗粒。本研究忽略了水流对温度场的影响， 仅考虑垂向上的热传导， 因此式（2）左侧第二项为0。

FEFLOW 7.1中的piFreeze模块将冰部分加入固相物质， 修改了式（2）中固相物质的导热系数与体积热容的计算。此时， ε_a+ε_w+ε_i+ε_r=1， 其中ε_a、 ε_w、 ε_i、 ε_r分别表示空气、 水、 冰和土壤颗粒的体积分数。piFreeze建立了冰与液态水之间的联系， 即未冻水含量与总的含水率之间的关系， 也称作冻结函数^［14］。

式中： ρ_w、 ρ_i分别为液态水和冰的密度。冻结函数F是与温度有关的一个函数， 有

式中： ［T_l， T_h］=［T₀-∆T/2， T₀+∆T/2］， 其中T₀为冻结温度。FEFLOW-piFreeze中假设当温度低于冻结温度时， 液态水并不立刻冻结成冰， 而是在给定的温度范围［T_l， T_h］内逐渐冻结。当温度低于T_l仍有部分液态水存在， 称作残余含水率φ_res。将冻结函数考虑在内， piFreeze对式（2）进行了以下修改。

固相物质导热系数改为

固相物质容积热容量改为

式中： L_f为冰形成时的潜热。FEFLOW模型的输入参数主要包括土壤剖面的温度和含水率的初始分布、 地表温度、 降水、 各层土壤水热参数。

本研究模拟了2006年9月2日至2008年12月31日活动层各深度土壤温度的动态变化过程， 其中2006年9月2日至2007年12月31日为模型校准期， 2008年1月1日至2008年12月31日为模型验证期。以2008年为基准， 分别按照0.02 ℃·a^-1、 0.048 ℃·a^-1、 0.07 ℃·a^-1的气温升高速率， 对2009年1月1日至2058年12月31日的活动层温度变化进行了预测。将模拟域从地表至地表以下600 cm深度分成10层， 代表性深度分别为0 cm、 5 cm、 10 cm、 15 cm、 20 cm、 40 cm、 80 cm、 160 cm、 320 cm、 600 cm。土壤剖面的温度和含水率的初始分布采用观测孔的观测值， 600 cm处土壤温度设为定值-3.5 ℃， 部分参数采用模型提供的默认值， 具体参数设置见表3。

模型的模拟性能可以通过一些统计参数进行评估。本研究采用相关系数（R）、 均方根误差（RMSE）、 平均绝对误差（MAE）、 纳什效率系数（NSE）对BP神经网络地表温度预测模型和FEFLOW模型的模拟性能进行评估。各统计参数的计算公式为

式中： n为观测样本数； M_i 、 M¯分别为模拟值及其平均值； O_i 、 O¯分别为观测值及其平均值。当R=1， RMSE=0， MAE=0， NSE=1时， 认为模型的模拟性能达到最佳。

以气温和辐射作为输入参数， 利用BP神经网络模拟预测得到的地表温度与观测值的动态变化过程如图6所示。可知， 地表温度模拟值的动态变化趋势、 峰谷值都与观测值基本一致， 很好地反映了地表温度一年内的动态变化过程。对模型验证期地表温度预测模型的模拟性能进行了统计分析， 其中模拟值与预测值的MAE、 RMSE分别为1.78 ℃、 1.41 ℃， 相关系数R高达0.98（P<0.01）， 模型的纳什效率系数NSE高达0.96。纳什效率系数常用来评估模型的模拟性能， 一般认为当0.6≤NSE<0.8时， 模拟结果是可以接受的； 当NSE≥0.8时， 可认为模型的模拟能力极佳^［21］。因此认为， 本文中由气温与总辐射作为驱动变量建立的地表温度预测模型可以很好地对地表温度进行模拟预测。由于多年冻土区地表温度监测点较少， 通过所建立的地表温度预测模型， 利用气象资料驱动， 可以很好地对本研究区过去与未来的地表温度进行模拟预测， 这将为下一步研究活动层不同深度土壤温度的预测和模拟奠定良好的基础。

在基于BP神经网络模型成功模拟预测地表温度的前提下， 利用改进的FEFLOW模型来模拟预测研究区活动层不同深度土壤温度的动态变化过程。图7是不同深度土壤温度模拟值与观测值年内的动态变化过程。从图中可以看出， 不同深度土壤温度的模拟值与观测值动态变化趋势一致， 二者吻合得很好， 模拟值能很好地反映活动层土壤剖面温度的动态变化过程。表4是不同深度土壤温度模拟的效率参数。各深度模拟值与观测值之间的相关系数R均高于0.94（P<0.01）， NSE值均高于0.93。模拟结果说明， FEFLOW模型对于模拟不同深度土壤温度的动态变化模拟精度较高， 效果良好。0 ~ 20 cm各层土壤温度模拟效果差别不大， RMSE值在1.73 ~ 1.76 ℃， MAE值在1.40 ~ 1.42 ℃； 自20 cm以下， 各层土壤RMSE值均小于1.24 ℃， MAE值均小于0.93 ℃， 即随着深度的增加， RMSE与MAE基本呈减小趋势， 说明随着土壤深度的增加模型的模拟精度逐渐提高。周剑等^［7］根据各深度土壤含水率的年内变化特征， 将风火山研究区土壤冻融过程划分为冻结发生期（9月下旬至11月中旬）、 完全冻结期（11月下旬至次年3月下旬）、 融化发生期（4月上旬至6月中旬）、 完全融化期（6月下旬至9月中旬）。本研究参照上述结论， 计算求得四个阶段各深度土壤温度平均RMSE值分别为1.26 ℃、 1.48 ℃、 0.94 ℃、 1.11 ℃， 可以看出模型在冻结发生期与完全冻结期的模拟精度均低于融化发生期与完全融化期， 这可能与模型计算的冻结期土壤水热参数存在误差有关， 从而导致了模型的不确定性。然而从整个模型的效率参数来看， FEFLOW模型能用于模拟多年冻土区土壤温度的动态变化过程， 且精度较高， 模拟效果好。

图8中（a）、 （b）、 （c）分别表示在0.02 ℃·a^-1、 0.048 ℃·a^-1、 0.07 ℃·a^-1升温速率下， 2008 - 2058年每隔10年的土壤剖面温度分布。从图中趋势线（黑色实线）可以看出， 活动层厚度（定义为0 ℃等温线所对应的最大深度， 图中为红色虚线对应的最深处）随着气温的升高逐渐增厚， 且升温速率越大， 活动层增厚速度越快。2058年时， 在0.02 ℃·a^-1、 0.048 ℃·a^-1、 0.07 ℃·a^-1升温速率下， 活动层厚度分别增加19.4、 51.8、 64.7 cm， 年平均增速分别为0.39 cm·a^-1、 1.04 cm·a^-1、 1.29 cm·a^-1。由图8可知， 随着气温的升高， 活动层各深度土壤开始融化的时间不断提前（逐渐融化阶段的0 ℃等温线不断前移）， 而开始冻结的时间不断向后推迟（逐渐冻结阶段的0 ℃等温线不断后移）， 说明土壤处于完全冻结状态的历时在不断缩减。气温升高速率越大， 开始融化时间越提前， 开始冻结时间也越向后推迟， 土壤完全冻结时间也越短。为了探究不同深度土壤对气温升高的响应， 本文选取10 cm、 80 cm、 160 cm三个深度， 分别代表浅、 中、 深层土壤， 对三个深度土壤温度与冻融特征变化进行了分析。

图9中（a）、 （b）、 （c）分别表示三种升温情景下， 预测得到的浅、 中、 深层土壤年平均温度。由图可知， 随着气温的升高， 浅、 中、 深层土壤温度也逐渐升高， 且升温幅度随土壤深度的增加而逐渐减小。以0.02 ℃·a^-1升温情景为例， 50年后浅、 中、 深层土壤温度分别升高了0.82 ℃、 0.47 ℃和0.35 ℃， 浅层土壤年平均升温速率（0.016 ℃·a^-1）与气温升高速率（0.02 ℃·a^-1）十分接近， 而深层土壤升温速率（0.007 ℃·a^-1）约为气温升高速率的1/3， 说明浅层土壤温度对气温变化的响应最为显著。在升温条件下， 不同深度土壤冻融过程的变化特征也印证了上述结论。由表5可知， 同一升温情景下， 随着土壤深度的增加， 冻结延迟时间与融化提前时间的变化幅度逐渐减小。以0.02 ℃·a^-1升温情景为例， 融化提前时间浅层（11 d）是深层（5 d）的2倍， 冻结延迟时间浅层（6 d）是深层（2 d）的3倍， 浅层土壤冻融过程对气温变化的响应最为强烈。对比各升温情景下融化提前时间与冻结延迟时间的变化可知， 气温升高对土壤融化过程的影响大于对冻结过程的影响（开始融化提前时间约为开始冻结延迟时间的2倍）。气温升高幅度越大， 土壤冻融过程变化的程度越剧烈， 但是随着气温升高速率的增大， 浅、 中、 深层土壤冻融过程变化的差异逐渐减小。

多位研究者采用不同的模型模拟了青藏高原多年冻土区的土壤热状况^{［4，27-28］}。在SHAW模型中， 辐射、 降水、 地表反射率、 植被生长信息等输入变量的偏差都可能会对土壤热通量的计算产生影响， 进而对土壤温度的模拟产生影响。此外， 该模型需要输入模拟开始和结束时的温湿度分布， 无法对气候变化情景下土壤温湿度的变化进行预测。Gao等^［27］指出， 复杂的地形、 植被覆盖、 土壤水分分布造成的辐射和地表反射率估计的不确定性导致GBEHM计算地表能量平衡时的不确定性， 进而使土壤温度的计算产生偏差。多位学者利用GIPL2模型成功模拟了多年冻土区土壤热状况的变化^［4，18］， 但GIPL2模型要求分别输入各层土壤冻结和融化期间的导热系数、 容积热容量等土壤热参数， 而FEFLOW可以根据冻融过程中的含冰率和未冻水含量， 自动计算土壤导热系数与容积热容量。表6统计整理并将本文模拟结果与其他模型模拟结果进行了比较（由于模拟深度不同， 取平均RMSE）。可知， FEFLOW模型平均RMSE小于另外三种（SHAW、 GBEHM、 GIPL2）模型， 模拟效果最好。综合来看， 以BP神经网络模型输出为输入的FEFLOW模型输入参数少、 模拟精度高， 更适用于多年冻土区土壤热状况的模拟。

在本研究中， FEFLOW模型对浅层土壤温度的模拟精度要低于深层土壤。原因可能来自多个方面： （1）研究点位于高寒草甸区， 植被根系深度约30 cm， 植被根部有机质的富集以及地表植被枯落物改变了浅层土壤的热性质， 而FEFLOW无植被输入选项， 势必会影响浅层土壤温度的模拟精度。（2）由表1可知各深度土壤粒径组成不同， 其热力性质必存在差异。本研究土壤热参数参考了文献［29］， 将各深度土壤作均质处理， 因此可能导致模拟出现偏差。由图8可知， FEFLOW模型在冻结阶段的模拟精度要低于融化阶段。Qin等^［30］在模拟青藏高原活动层厚度与多年冻土热状况时， 对土壤含水率、 土壤热参数（如融化和冻结阶段的导热系数、 热容）以及驱动数据进行了不确定性分析， 指出活动层厚度对冻结阶段的土壤导热系数最为敏感。在本研究中， 冻结和融化阶段温度模拟效果的差异可能与FEFLOW在计算冻结与融化阶段的水热参数存在误差， 造成模型的不确定性有关。土壤结构参数、 水热动力学参数的准确描述是多年冻土区活动层土壤冻融过程中水热传输模拟的关键， 未来可考虑适用于青藏高原多年冻土区的土壤水热性质参数化方案， 以提高模型的模拟精度。

在全球气候变暖大背景下， 青藏高原多年冻土大面积退化， 土壤温度不断升高， 多年冻土上限不断加深， 导致活动层厚度逐渐增大^［6，27］。本研究使用FEFLOW模型对三种不同升温情景下风火山区域活动层热状况的动态变化进行了单点尺度上的数值模拟预测。由图9可知， 随着气温的升高， 活动层厚度逐渐增厚， 0.02 ℃·a^-1、 0.048 ℃·a^-1、 0.07 ℃·a^-1三种升温情景下， 50年后（2058年）活动层厚度分别增加了19.4、 51.8、 64.7 cm， 年平均增速分别为0.39 cm·a^-1、 1.04 cm·a^-1、 1.29 cm·a^-1。张中琼等^［31］模拟预测了不同升温情景下青藏高原活动层厚度变化情况， 0.023 ℃·a^-1升温情景下活动层增厚速率为0.50 ~ 1.25 cm·a^-1； 0.052 ℃·a^-1升温情景下为0.50 ~ 1.37 cm·a^-1。Pang等^［32］模拟结果表明， 在0.018 ℃·a^-1升温情景下， 青藏高原活动层增厚速率为0.27 ~ 0.81 cm·a^-1； 0.04 ℃·a^-1升温情景下为0.81 ~ 1.89 cm·a^-1。非连续型多年冻土的平均温度偏高， 对于气温变化更为敏感^［32］， 本次模拟结果比上述结果稍微偏小， 可能与研究区位于青藏高原腹地， 属于连续型多年冻土有关。Rasmussen等^［33］模拟了升温情景下格陵兰岛东北部肯伯格山谷的活动层厚度变化， 在0.03 ℃·a^-1升温情景下， 活动层增厚速率为0.67 cm·a^-1， 比青藏高原活动层增厚速率小， 青藏高原地区多年冻土对气候变化更为敏感。

张明礼等^［34］等模拟了气温升高情景下青藏高原北麓河地区土壤温度的变化过程， 结果表明气温升高对土壤温度的影响随深度的增加逐渐减弱。刘光生^［28］指出土壤温度与气温的相关性随深度的增加而逐渐降低， 土壤温度对气温的响应程度随深度的增加而逐渐降低。在本研究中， 通过对比0.02 ℃·a^-1、 0.048 ℃·a^-1、 0.07 ℃·a^-1三种升温情景下不同深度土壤温度的变化， 发现土壤温度升高的幅度均随着土壤深度的增加而逐渐减小， 不同升温情景下温度的差异随深度增加而减小， 与上述研究结论相符。

气温升高将对土壤的冻融过程产生显著影响。Wang等^［35］利用开顶式气室（OTCs）升高近地表层气温， 评估了气温升高对土壤热状况的影响。结果表明， 近地表气温的升高延长了土壤的解冻时间， 延缓了秋季冻结过程的发生， 缩短了连续冻结的持续时间， 且这种影响随着土壤深度的增加而逐渐减小。刘光生^［28］指出， 随着气温的升高， 完全冻结期间地温变化幅度趋缓， 持续冻结时间逐渐减小； 未冻结期间地温变化速率逐渐增大， 持续解冻时间逐渐增大。气温升高导致的土壤融化时间提前、 冻结时间推迟， 增加了植被生长的有效积温并延长了植被生长期， 有利于植被在冬季的生长和生存^［35］。本研究中， 在气温升高情景下， 各深度开始融化时间提前， 开始冻结时间延迟， 且土壤的开始融化时间的提前天数、 开始冻结时间的延迟天数以及土壤处于完全冻结状态缩短天数均随着土壤深度的增加而逐渐减小， 与上述研究结论相符。

利用青藏高原风火山地区2006 - 2008年气象资料建立了BP神经网络地表温度预测模型， 对地表温度进行了模拟重现， 并对0.02 ℃·a^-1、 0.048 ℃·a^-1、 0.07 ℃·a^-1三种升温情景下的地表温度进行了模拟预测。根据风火山地区的气象与地温、 水分等观测资料， 利用FEFLOW模型在单点尺度上对活动层内各深度土壤的温度进行了模拟。结果表明， 两个模型的模拟精度高， 模拟效果好， 可用于青藏高原多年冻土活动层土壤温度动态变化过程的预测与模拟， 并对模拟结果进行了分析讨论。得到如下结论：

（1） 利用较易获得的气象资料所建立的BP神经网络地表温度预测模型与FEFLOW模型能够很好地模拟多年冻土区土壤热状况的动态变化过程。虽然FEFLOW模型对浅层土壤温度的模拟精度低于深层土壤、 冻结过程的模拟精度低于融化过程， 但模拟的不同深度土壤温度变化的趋势、 峰谷值都与观测值基本一致， 较好地体现了土壤温度的季节变化。通过对比FEFLOW模型与其他模型（SHAW、 GIPL2、 GBEHM）在青藏高原多年冻土区模拟的土壤温度结果， 发现以BP神经网络模型输出为输入的FEFLOW模型输入参数少、 模拟的各深度土壤温度更接近实际观测值。因此， 对青藏高原土壤热状况的模拟更具有优势， 能够成为模拟多年冻土区土壤热状况的理想模型。

（2） 在0.02 ℃·a^-1、 0.048 ℃·a^-1、 0.07 ℃·a^-1三种升温情景下， 50年后（2058年）风火山地区活动层厚度分别增加了19.4、 51.8、 64.7 cm； 各深度土壤温度均有不同程度升高， 开始融化时间提前， 开始冻结时间向后延迟， 不同深度土壤冻融过程对气温升高的响应差异随着升温幅度的增大而减小。其中， 浅层土壤对气温变化的响应最为显著， 升温幅度最大， 冻融起始时间变化幅度也最大， 且这种变化随着土壤深度的增加而逐渐减小。同时， 气温升高对土壤融化过程的影响大于冻结过程， 各深度土壤融化提前时间约为冻结延迟时间的2倍。

近几十年来， 青藏高原暖湿化趋势显著^［36］。气温的升高伴随着降水的增加， 降水可降低表层土壤温度， 减小土壤热传导通量， 有利于冻土的热稳定性^［36］， 也有研究认为降水的增加会使得活动层厚度增大^［37］， 有关降水对活动层热状况的影响还存在争议。本研究中， 模拟未来气候变化对活动层热状况（包括地表温度）的影响时， 仅考虑了气温的变化， 忽略了降水的影响， 这势必会对模拟结果产生影响。未来预测气候变化对活动层热状况的影响时， 应同时考虑降水和气温两种因素， 才能更好地理解、 评估和预测多年冻土区活动层热状况对气候变化的响应， 以及可能对地-气系统之间物质和能量交换过程产生的影响。