Stefan方程在土壤冻融过程模拟中的应用

刘文惠， 谢昌卫， 刘海瑞， 庞强强， 王 武，刘广岳， 杨雨昆， 王 铭， 张 琪

（1.青海大学地质工程系，青海西宁 810016； 2.中国科学院西北生态环境资源研究院冰冻圈科学国家重点实验室藏北高原冰冻圈特殊环境与灾害国家野外科学观测研究站，甘肃兰州 730000； 3.青海大学生态环境工程学院，青海西宁 810016）

0 引言

多年冻土是指温度在0 ℃或低于0 ℃至少连续存在两年的岩土层[1］。多年冻土是冰冻圈的重要组成部分，其影响能量交换、水文过程和生态环境，进而影响全球气候系统。活动层是指覆盖于多年冻土之上的夏季融化冬季冻结的土层[1］。活动层的水热动态变化过程影响着冻土区水文和生态系统的生物、物理及地球化学过程。多年冻土与大气圈之间的相互作用主要是通过活动层中的水、热动态变化过程而实现的。多年冻土和活动层研究很早就受到了国内外学者的广泛关注，并加大了这方面的观测和研究力度。多年冻土与活动层被世界气候研 究 计 划（World Climate Research Programme，WCRP）列入“气候与冰冻圈计划”（Climate and Cryosphere，CliC）的主要观测研究内容之一。由国际冻土协会和加拿大地质调查局联合发起的全球多年冻土监测网络（Global Terrestrial Network for Permafrost，GTN-P）用于监测全球多年冻土地温。环北极活动层监测网络（Circumpolar Active Layer Monitoring Network，CALM）旨在对活动层厚度和热状况进行监测。活动层作为多年冻土与大气圈进行热、质交换的最主要场所和媒介。活动层厚度又是判断多年冻土退化最为直观的标志，其可以更好地反映多年冻土对气候变化的响应状况。随着全球变暖加剧，多年冻土退化严重[2-9］，尤其是活动层厚度增厚尤为明显[4-6］。因此，气候变化背景下的活动层冻融过程模拟、厚度制图和变化预测是研究冻土区生态环境、水文、工程以及碳循环的基础，也是目前冻土学领域的研究热点。活动层厚度可以通过钻探、坑探等方法直接测量，也可用测温法、地球物理勘探等间接方法来估计。但是，由于监测条件的限制和监测数据的有限，实地监测很难满足活动层厚度空间分布的模拟需要，大尺度活动层厚度空间分布以及长时间序列活动层对气候变化的响应还是得依靠模型来解决。

用于模拟活动层冻融过程和厚度的模型较多，主要分为两大类：经验半经验公式和以求解热传导方程为基础的数学物理方法[10-11］。由于数学物理方法的初始边界条件在时空上变化较大，因此，常用于计算单点的冻结融化深度。对于空间上高度非均匀的冻土空间的变化，经验和半经验型公式是一个较为实际的选择，尤其是半经验半理论的计算方案，既充分考虑了冻土的物理特性和过程，又可以用大尺度参量获取冻土在水平和垂直空间的变化，更适宜于冻土的实际研究[11］。目前，Stefan 方程是国内外用于计算多年冻土冻结和融化深度最常用的经验公式。它充分考虑了气候条件、土壤热属性和水分条件，形式简单，驱动参数少，模拟效果较好，既可以用于模拟单点的冻结融化深度，也可以较方便地模拟大尺度的活动层厚度空间分布。本文就不同修正形式的Stefan方程在国内外多年冻土活动层冻融过程模拟中的研究进展和应用中存在的问题进行了探讨，旨在为今后的研究提供参考。

1 Stefan方程介绍

斯蒂芬方程是奥地利科学家Josef Stefan 在研究北极地区湖冰形成过程时提出的[12］。假设冰体内热量传导非常迅速，并且冰体内的温度变化是线性的。当冰的表面温度低于相变温度时，冰体下部与湖水接触面处的温度则等于相变温度。在一个给定的时间内，冰体从下面湖水中得到的热量和其从表面排出的热量是相等的，这一关系可以表达为：

由此得出冰体厚度随时间的公式：

事实上，式（1）仅可以描述冰体内简化的热传导过程，对于常规的热传导过程，通常用一维傅力叶热传导方程描述。在给定的边界条件和假设下，Nemman 和Stefan 先后给出了这一非齐次二阶微分方程的解。当假设冻结体内温度梯度为线性的条件下，Stefan 解仍然可以简化成式（2）的形式。因此，一般认为Stefan 方程是一维热传导方程Stefan解简化后的公式表达。

Stefan 方程提出后的数十年里面，一般用于计算湖冰、海冰等冰体的厚度。1943年，Berggren[13］对式（2）中纯冰体的热容用土壤中冰体的热容代替，将Stefan方程应用到浅层土层冻结和融化过程的计算中。在理论上，冰-水相变伴随的潜热释放或者吸收要远大于干土本身热容的变化，因此，当土层内含水量较大时，将Stefan 方程应用到土壤中与应用到冰体内造成的差别并不大。由此，式（1）转化为式（3），即为广泛应用于冻土学中的Stefan方程的通用形式[14］。

式中：Z 为冻结/融化深度（m）；K 为导热系数（W·m-1·℃-1）；QL为土壤水分相变引起的潜热变化，QL=Lρ(ω −ωu)，L 为冰的融化潜热（3.3×105J·kg-1）；ρ为土壤的干容重（kg·m-3）；ω 为总的含水量（%）；ωu为未冻水含量（%）；DDF 和DDT分别为冻结指数和融化指数（℃·d），指地表温度日均值连续低于0 ℃的温度累计之和地表温度日均值连续高于0 ℃的温度累计之和；Tt为日平均地表温度（℃）；T0= 0。

Stefan 方程首次将地表（或者大气）温度的变化与冰层（或者土层）的冻结融化过程以简单公式的形式联系起来，极大地简化了土壤冻结融化深度的分析计算过程。因此，Stefan 方程在推出后被广泛运用到寒冷地区的工程建设中，并逐渐被应用到冰川消融、冻土变化研究中[15-18］。目前，在多年冻土活动层冻融过程的研究中，Stefan 方程仍然是应用最广泛的模拟计算方法。

2 Stefan方程在多年冻土中的应用

2.1 在均质土壤中的应用

2.1.1 引入N因子

Stefan 方程最初被运用到土层冻融过程时，为了避免因未考虑外部热量交换因素引起的误差，建议采用地表以下5~10 cm 的日平均温度作为温度驱动[19］。由于日尺度该深度的日平均温度获取困难，一般用日平均地表温度计算冻结-融化指数。相比地表温度或者一定深度内土层温度，气温更容易观测，实际应用中则大多采用气温来计算计算冻结-融化指数。这在一定程度上造成了额外的计算误差。针对这种状况，Carlson[20］最早提出N 因子，对Stefan方程做了进一步的修正：

式中：Nt为融化N 因子，是地面融化指数与气温融化指数之比；Nf为冻结N 因子，是地面冻结指数与气温冻结指数之比。从概念上讲，N 因子用一个简单的数值代表了特定下垫面地表能量通量的总和。一般情况下，融化N 因子大于冻结N 因子。不同下垫面类型的N因子不相同，但是，同一区域相同地表条件的N 因子相差不大，基本可以取固定值。从表1的统计结果可以发现，青藏高原的融化N 因子和冻结N因子普遍大于达拉斯加和加拿大等高纬度地区的融化N 因子和冻结N 因子；土壤颗粒越粗，融化N因子越大。裸地融化N 因子最大（＞1.28），其次是灌丛和苔原，河道和沼泽的融化N 因子最小（＜0.7）[31］。北极地区的N 因子值比较稳定，其年变化幅度小于10%[32］。山地多年冻土区的N因子值年际变化比较大，尤其是冬季随着积雪的年变化而变化较大[33］。

表1 融化N因子（Nt）和冻结N因子（Nf）值统计结果Table 1 The thawing N-factors（Nt）and freezing N-factors（Nf）

引入N 因子后的Stefan 方程被广泛用于更复杂地形和下垫面的大尺度活动层厚度空间分布及变化预测模拟。庞强强等[34］模拟了青藏高原多年冻土活动层厚度的空间分布状况。徐晓明等[35］得到青藏高原多年冻土活动层厚度平均为2.39 m，活动层厚度在羌塘盆地最小，在多年冻土区边缘、祁连山、西昆仑山、念青唐古拉山较大。气候变化条件下，活动层厚度呈整体增大趋势，1981—2010 年，活动层厚度的变化量为-1.54～2.24 m，变化率为-5.90～10.13 cm·a-1，平均每年变化1.29 cm。张中琼等[36］模拟了青藏高原活动层厚度空间分布状况并预测了未来不同气候情景下活动层厚度的变化情况。到2050年，A1B 情景下活动层厚度最大达10.20 m，增加约0.30~0.80 m；B1 情景下增加0.20~0.50 m；A2 情景下增加0.20~0.55 m。Klene等[37］完成了库帕河流域活动层厚度的空间分布制图，Shur 等[32］完成了俄罗斯多年冻土分布制图，张廷军等[38］模拟得到俄罗斯北极地区鄂毕河、叶尼塞河和勒拿河流域平均活动层厚度分别为1.80 m、1.70 m和1.66 m。

引入N 因子后的Stefan 方程的模拟精度也得到大大提高。Klene 等[37］进行库帕河流域活动层厚度空间分布制图时用地表温度作为温度驱动得到的模拟误差为14.5%，用3 年的平均N 因子和气温得到的误差为17.6%，用气温得到的误差为29.2%。同时，根据不同温度驱动估算马衔山2010年的活动层厚度，结果显示用地表温度作为温度驱动得到的活动层厚度为1.13 m，相对误差最小（4%），最接近实测值；用气温得到的活动层厚度为1.01 m，相对误差最大（13.6%），明显偏小；而用气温和N 因子得到的活动层厚度为1.16 m，偏大，但是误差要远远小于只用气温驱动得到的结果。地面温度是土壤通过地面与大气间热交换特征的综合指标，比气温更能体现冻土的热状况。而且，地表温度对外界条件反应更为迅速，能较为准确地反映冻土的上边界热状况，是众多冻土模型最佳的选择。因此，在完成大尺度活动层厚度空间分布模拟时，当地表温度缺失的情况下，为了尽可能减小模拟误差，可以考虑用N因子和气温作为温度驱动。

2.1.2 引入E值

Harlan 和Nixon 认为Stefan 方程可以表示为土壤特征与融化指数之间的线性函数[39］：

式中：E 为与土壤热参数、含水量和地表覆盖类型有关的比例参数，表示活动层的融化速率。已知任意点的实测活动层厚度和融化指数，可以得到E 值。已知E 值和其他土壤参数的情况下，便可得到其中的任一土壤热参数。Hinkel等[40］得到阿拉斯加北部森林多年冻土区1992 年和1993 年的“E”值分别为1.22和1.21。假定土壤孔隙度分别为0.5和0.6时，得到导热系数分别为0.126和0.151 W·m-1·℃-1。

已知E值和融化指数就可以得到活动层厚度的空间分布状况。Peng 等[41］得到黑河流域的E值范围为0.028~0.053，其中砂砾石E值最大，裸岩次之，荒漠的最小。并结合气温得到黑河流域2000—2008 的年平均冻结深度的空间分布范围为1.0~3.5 m。Shiklomanov 等[42-43］成了阿拉斯加库帕河1987—1999 年的年平均活动层厚度高精度制图（50 m 分辨率），并发现在年变化尺度上，活动层厚度与融化指数具有高度一致性。1989 年和1998 年的融化指数最大，对应最大的活动层厚度，1991 年和1996 年的融化指数最小，对应最小的活动层厚度。Brown等[44］和Hinkel等[45］基于极地活动层监测数据和气温融化指数做拟合分析得到，活动层厚度与融化指数具有很好的对应关系，阿拉斯加、加拿大北部、北欧地区和俄罗斯地区的一些监测点均显示活动层厚度在1998年达到最大值（对应最大融化指数），在2000年达到最小值（对应最小融化指数）。

2.1.3 加入降雨、地形因子

山地多年冻土的形成和发展与局地微气候和微地形（坡度和坡向）有很大的密切关系。基于此，一些学者在上式中加入了降雨和地形因子。Hinkel等[40］将阿拉斯加北方森林多年冻土区的实测活动层厚度和融化指数做最小二乘回归拟合后得到Ste⁃fan方程的另一个形式：

认为α项表示为降雨对活动层厚度的影响。暖季一天内不同时刻的降雨对活动层会产生不同的影响效果。日出前后，气温很低，此时形成的降水雨水温度较低（可称之为冷降雨事件），而形成于午后气温相对较高时的降水温度较高，可称之为暖降雨事件。其中，暖降雨雨水温度高，降落地面后会带入大量感热进入深层土壤，从而加速了活动层融化速率，使得活动层厚度增厚。但是，这个研究结果只限于阿拉斯加森林环境下泥炭层较厚的多年冻土区，其他地区没有尝试过。Shiklomanov 等[42］完成库帕河流域活动层厚度高精度制图时用降水代替含水量（该区域蒸发很小，假定降水全部用于下渗）发现活动层融化速率与前一年融化期末的降水有很大的相关性，1996 年融化期末降水最少，对应1997年的活动层融化速率最大。

Nelson 等[46］在上式中加入了地形因素，完成了阿拉斯加库帕河流域更详尽的活动层厚度分布制图（1 km 分辨率），得到北坡（阴面）的活动层厚度在减小（负值），南坡（阳面）的活动层厚度在增大。

式中：r 为地形因子引起的辐射因素，是个无量纲参数；Rs为斜坡上的潜在太阳辐射，与纬度、坡度和坡向相关；Rh为水平面上的潜在太阳辐射。

2.2 分层土壤中的应用

土壤不同于冰体那样是由均质成分构成，在不同深度岩土成分、结构以及水分条件等通常有较大的差异。将本来用于计算均质冰体冻融深度的Ste⁃fan方程应用到非均质的岩土中，必然会引起一定的计算误差。因此，不同时期的不同学者提出了一些将Stefan方程用于计算分层堆积土壤冻融深度的算法，如被广泛应用的J-L算法[47-49］、N-M算法[50］、分层总和法[51］和X-G算法[52］。

2.2.1 J-L算法

J-L 算法（也被称为St Paul 方程）由Jumikis 和Lunardini自1950年代提出，在工程计算中得到了广泛的应用，并被应用到许多大型的数学模型中。原理是通过计算冻融到第n层所需要的冻融指数来得到最大冻融深度，其推导过程见参考文献[47-48］。Woo 等[53-54］利用J-L 算法模拟了加拿大西北部马更些河流域的活动层厚度并分析了有机质层的影响。当泥炭层厚度为0.2 m 时，A2 情景下，苔原活动层厚度由0.68 m 增大到2100 年的0.96 m；森林活动层由1.35 m 增大到2100 年的1.92 m。当泥炭层厚度增厚到1 m 时，A2 情境下苔原的活动层厚度由0.39 m 增大到2100 年的0.50 m；森林活动层由0.65 m 增大到2100 年的0.75 m。同时，分析了土壤质地、含水量和温度等对该区域不连续多年冻土活动层厚度的影响[55］。Woo等[56］于2004将J-L算法做了改进后用于模拟双向的冻结融化过程，得到北美自北向南多年冻土区到季节冻土区6种下垫面的冻结融化过程。对比单向的冻结融化过程，双向的冻结融化过程大大提高了模拟精度，误差范围为0.16~0.58，而单向模拟的结果误差高达0.53~1.34。

2.2.2 分层总和法

分层总和法是通过分别计算冻融时各层土所消耗的冻融指数得到各层的冻结深度，最后各层的冻融深度之和为最大冻融深度。其推导过程见参考文献[51］。分层总和法的研究较少，应用没有得到推广。

Stefan 方程通用形式[式（3）］可以改写成如下形式：

从式（8）中可以看出，冻融指数可以看成是冻融深度的幂函数，即冻融指数与冻融深度的平方成正比，对于厚度为z = z1+ z2的均质土壤，如下关系式成立：

从式（9）可以看出，即使对均质土壤，要冻结/融化厚度为z 的土壤，其所需要消耗的冻结/融化指数要大于分别冻结/融化厚度为z1和z2的土壤所需要的冻融指数之和，因此，简单地利用分层总和法估算非均质土壤冻结/融化深度是错误的。

2.2.3 X-G算法

2013 年Xie 等[52］利用类比递推的方法，提出了基于Stefan方程计算分层土壤冻融深度的简易算法（X-G 算法）。X-G 算法能够计算由任意多层不同厚度的土层组成的土壤的冻融过程，是目前唯一能将斯蒂芬方程运用到非均质土壤冻融过程的算法。它基于不同土层的物理参数将斯蒂芬方程应用到估算任意多层、任意厚度的土壤冻结/融化深度，而不是将不同层土壤参数进行平均处理，这样避免了不必要的计算误差。X-G 算法原理简单，既可模拟单向冻结过程也可模拟双向冻结过程，目前在模拟国内多年冻土冻融过程方面取得了较好的进展。Xie 等[52］模拟得到马衔山、两道河、昆龙山垭口和西大滩的活动层厚度分别为1.12 m、1.23 m、1.72 m和1.86 m，与实测值的相对误差范围为4.2%~17%。刘文惠等[57］利用该算法模拟了马衔山2010—2013 年的活动层厚度，模拟值分别为1.06 m、99 m、1.16 m 和1.04 m，与实测值很接近，相对误差范围为1%~9%。

X-G 算法发表以后，在相关领域受到了广泛的关注，Kurylyk[58］认为X-G 算法在数学上是合理的，但不能很好解释土壤冻融过程中温度变化的过程，并认为N-M 算法能更好地将Stefan 方程应用到分层堆积的土壤。从非线性傅里叶热传导方程的斯蒂芬解入手，对Stefan 方程基本原理、X-G 算法、J-L算法和N-M 算法推导过程等进行了如下详细分析[59］：

分层堆积的土壤中冻结融化深度作为时间的分段函数，其函数形式为：

系数m 是土壤物理性质参数的函数。对上式求导数，得到：

可见其导数形式也是不连续的。分层堆积土壤中的冻融过程是时间的分段函数，并不是连续函数，从而最终证明基于连续函数推导得到，即建立在连续微分形式上的J-L 算法和N-M 算法是错误的，X-G 算法是目前唯一能将斯蒂芬方程运用到分层堆积土壤的算法。

3 与陆面和水文模型的耦合

目前的陆面过程模型和水文模型的研究对均匀、覆盖稠密的下垫面条件研究比较成熟，且较多，但忽略了冻土和积雪等复杂因素或仅做了十分简单的概化。考虑各圈层相互作用，发展多圈层综合陆面过程面模型和水文模型成为一种必要。多年冻土对地气之间能量交换、水循环有很重要的影响，活动层冻融过程的准确模拟可以很好地研究多年冻土区地-气-能-水的交换过程。因此，在陆面过程面模型和水文模型中耦合冻土模型，反映气候变化背景下土壤冻融过程中水热过程迁移等对陆地-大气水热交换过程和寒区环境水文过程的影响是目前大气和水文学者的关注重点。许多研究者也在陆面过程模型框架下发展了冻土参数化方案，发现在陆面水文过程模型中考虑冻土作用能显著地增强模型模拟能力[60-63］。不少研究人员尝试在分布式水文模型中添加冻土过程，以适应寒区环境水文过程模拟。研究表明，包含冻土冻融模块的模型能够成功模拟由冻土融化引起的春季洪水过程，而运用不具备冻土模拟功能的分布水文模型在冻土明显的寒区流域进行模拟时，不能准确捕捉融雪径流过程，由于缺失冻土模块而严重低估融雪期间的径流洪峰，但暖期径流量又明显偏高[64-67］。

Stefan 方程由于其原理简单、计算简便、所需要参数易于获得，在一些陆面过程和水文模型中得到应用。多年冻土与气候模式的耦合有两种方法，一种是将冻结融化过程融合到气候模式中；一种是基于气候模式的温度数据作为冻结融化过程的温度驱动，这种方法目前应用较多，比较成熟。Stendel等[68］基 于OAGCMs（ocean-atmosphere general cir⁃culation models）和Stefan 方程得到当前气候条件下（1961—1990 年）北半球多年冻土活动层厚度范围为0~1.2 m；IPCCA2情景模式下2071—2100年活动层厚度将增加1.1~1.6 倍，其中俄罗斯北极地区活动层厚度将增加30%~40%，增加幅度最大的是西伯利亚东北部和加拿大西部，其次是中国和蒙古北部。Li等[69］将Stefan 方程的常用表达式耦合到陆面模式SiB2 中发现，加入冻土参数化后的模拟精度大大提高，模拟结果误差不到9%。其中，模拟得到的表层和底部土壤湿度的绝对误差分别为0.020和0.013，远远小于没有加入冻土参数化的模拟结果。Yi等[70］考虑泥炭层、未冻水和冠层储热后将双向的J-L算法耦合到陆面模式CLM3中，对比没有加入J-L算法的模拟结果，改进后的CLM3 消除了冻结融化锋面在3、4月份的较大波动以及在1月和3月出现稽延期后较大的跳跃下降现象；并且在模拟土壤温度、土壤含水量和雪深方面的精度大大提高。Fox[71］将J-L算法融合到水文模型中去，基于土壤水热相互作用，分析了冻结融化过程对土壤水量平衡要素的影响，并通过敏感性实验分析得出冻融过程对土壤径流变化和预测极端径流事件具有重要的潜在影响。与最新版本的大尺度寒区水文模型（cold region hydro⁃logical model，CRHM）冻土模块相比，尽管耦合Ste⁃fan 方程后的水文模型在描述冻土活动层冻融过程对水文过程的影响方面有较大优势，但不能准确描述冻土融水过程、水分迁移入渗过程、冻土坡面汇流以及产汇流过程等，更不能定量化地表径流和壤中流的量及冻土融水对径流的贡献值。

4 存在的问题

4.1 模拟误差分析及可能的原因

土壤中的冻融过程是非常复杂的，许多因素起着关键的影响作用。Stefan 方程以一维热传导方程为理论基础，假设地面吸收的热量全部用于多年冻土的融化且温度变化是线性得出的，这显然与实际冻融过程不一致。因此，利用Stefan 方程模拟得到的冻结融化深度与实际冻融深度之间通常有或大或小的误差。一般认为，Stefan 方程忽略了感热变化，没有考虑外部热量交换和冻结岩层与下覆融土层的热量交换，最终导致计算结果偏大。但是，在实际的应用过程中却出现了模拟值偏小的情况，如蒙古北部地区和青藏高原的个别点（表2）。针对这种情况，可能的解释是同一土壤剖面内不同深度的土壤含水量、导热率和容重等参数存在较大的差异，将这些参数取平均值或者依照一层同性选取参数有可能拉低了整个参数取值，造成计算结果有一些偏小。

表2 基于Stefan方程的活动层厚度估算值小于实测值的结果统计Table 2 The relative errors between calculated active layer thickness and measured active layer thickness

冻融指数和导热系数作为Stefan 方程两个重要输入参数，其计算的准确性和获取方法差异影响模拟结果，进而引起一定的模拟误差。冻融指数在一定程度上可指示冻融作用的深度、强度及持续时间，其变化深刻影响冻融作用下形成的冰缘地貌和寒区地质环境[72］。冻融指数计算方法包括经典日计算方法、月平均方法和年振幅方法[23］。冻融指数计算时针对温度数据缺测时间长短采用不同的插补方式，缺测1天，选择其前后各一天数值取平均值插补；缺测2天，缺测第一天选择该日前两天的数值取平均，缺测第二天选择该日后两天的数值取平均；缺测超过2天但不超过一个月，选择前后各一年该月数据进行插值补充[73］。此外，不同学者针对不同研究区域的冻结期和融化期开始结束时间的界定方法不同。冻融指数不同计算方法和温度数据缺测不同插补方式得到的冻融指数有一定的差异，进而影响冻融指数的计算准确性。冻土导热系数反映了地层冻结过程中包括土中相变潜热影响的综合导热能力，直接影响冻结冷量在地层中的传递过程。多年冻土导热系数基于野外测试方法和计算模型得到。现存的适用于冻土区的导热系数计算模型多以一种或几种土壤条件为前提，或者多考虑局地因素影响。同时，多年冻土区土壤受冻融循环影响较大，多年冻土内部水热传输过程复杂，模型没有考虑未冻水含量、土骨架组成及冻土结构等对冻土导热系数影响，使得模型得到的导热系数精度不高、适用性具有局限性[74-75］。

4.2 没有考虑降水的影响

一般来讲，降雨增加，活动层含水量增加，土壤冻结需要释放更大的潜热，使得土壤冻结过程中的温度降低大大减小；同时，降雨减少，秋季冻结的冰含量减少，来年用于融化冰消耗的热量减少，更多的热量用于加热活动层，使得活动层融化速率增大。然而，Subin 等[76］却认为降雨会增大地表感热传递和土体融化潜热量，使得活动层和多年冻土温度升高。

降雨对活动层水热的影响比较复杂。不同区域的影响不同；同一区域不同降雨强度和降雨时长、一年中暖季降雨和冷季降雨以及一天中不同时刻的暖降雨和冷降雨均对活动层水热有影响且影响机理不同。国外有关这方面的研究仅在阿拉斯加北方森林有过。青藏高原降水较多，主要集中在5—9 月，东部降水多，西部降水少[77］。降水的这种时空分布不均匀性对活动层水热影响差异较大。目前，除了北麓河地区外，整个青藏高原上有关这方面的研究至今是空白。张明礼等[78-79］在北麓河的研究认为夏季短期、高频次降雨有减少地表净辐射、增加地表蒸发潜热、降低土壤表层温度的作用。李德生等[80］在北麓河地区发现暖季的高频率、小雨量降雨对活动层具有显著的冷却效果，且凌晨2 点左右的降雨产生能量变量很小，而14点左右的降雨产生的能量变量最大。Wen 等[81］通过冻土监测发现，北麓河夏季降雨入渗对流作用降低了地表温度梯度、减小土壤热通量，降雨增加减缓冻土的退化。因此，可以尝试在Stefan 方程中加入降雨后定量分析降雨尤其是极端降雨事件对活动层融化速率的影响。

4.3 含水量的敏感性

土壤含水量和未冻水含量作为Stefan 方程的重要输入项，考虑了冰水二相转换释放的潜热对融化深度的影响。当其他条件恒定，土壤含水量发生变化时，季节冻结和季节融化深度相应变化。Woo等[56］也发现相比干容重、有机质和矿物质含量，J-L算法对地表温度和土壤含水量更敏感，得到当地表温度升高1 ℃时，最大冻结深度减小了0.22 m，当土壤重量含水量分别增加50%和150%时最大冻结深度减小了0.06 m和0.24 m。但是，在少数研究中往往没有考虑土壤含水量和未冻水含量的影响。Ro⁃manovsky 等[18］忽略含水量后模拟阿拉斯加北极海岸平原自北向南的West Dock、Franklin 和Franklin Bluffs 三个点1987—1992 年的活动层厚度，得到的最大误差分别高达26%、3%和71%，存在明显高估的情况。Stefan 方程对含水量具有高度依赖性。在水分较充足的区域模拟得到的误差较小，在相对干旱的区域得到的误差较大。正如图1 所示，模拟值与实测值的相对误差随含水量的增加而降低。这符合Stefan 方程理论上的定义，Stefan 方程只适用于计算湖水（冰）的冻结（融化）厚度，如果用以计算土层的冻结或融化深度，则其计算结果的误差随土层含水（冰）量的减少而增大。

图1 青藏高原活动层厚度模拟值与实测值的相对误差与土壤含水量的关系Fig.1 Relationship between relative error and soil water content of permafrost on the Qinghai-Tibet Plateau

在多年冻土区活动层底部及上限附近最冷时期存在着较高的未冻水含量，这部分未冻水并没有参与实际的冻结融化过程。在计算中应该除掉这部分未冻水含量。Wilhelm 等[82］估算了南极半岛西部阿姆斯勒岛活动层厚度，由于该区域含水量比较少（重量含水量0~8.4%），忽略了含水量后得到的活动层厚度明显偏大。Uxa[83］对Wilhelm 等的模拟结果做了纠正，发现该区域多年冻土中还存在2%的未冻水，假定这部分未冻水参与冻融过程，得到的活动层厚度更接近实测值。基于X-G 算法不考虑未冻水得到马衔山2010—2013 四年的活动层厚度均小于实测厚度，但实际上马衔山多年冻土区活动层下部最冷时期始终保持0.1 m3·m-3的体积未冻水。将这部分体积未冻水转换成重量未冻水并假定其不参与融化过程，得到的活动层厚度为118 cm；而当以总的含水量为输入项开展模拟时，由于高估了冻融过程中冰水二相转化消耗的潜热，导致模拟深度减小，模拟的活动层厚度为104 cm，明显小于未冻水不参与的情况。不同土壤质地的未冻水含量是不相同的，同一土壤类型的未冻水含量随温度而变化。一般情况下，由于缺乏土壤水分的观测资料，尤其是未冻水含量的监测更少，很难满足冻土空间分布上的确定，在计算过程中未冻水含量常取固定值，必定使得模拟结果存在较大的误差和不准确性。

5 讨论与结论

针对Stefan 方程进行的多种改进措施虽然增强了方程对影响冻融过程的因素的包容性，但改进后的方程往往限制更多，不同地区应用时要引用更多的经验性因子。事实上，由于影响冻融过程因素很多，复杂的方程并不一定会取得更好的模拟效果。如俄罗斯冻土学家Kudryavtsev 于1974 年提出的Kudryavtsev 方程[84］是在冻土学中与Stefan 方程并列的计算活动层冻融深度的重要方法。该方程既可以模拟活动层厚度，也可以模拟多年冻土顶板温度。计算中不仅考虑了潜热，而且考虑了积雪、土壤的热传导和热容量效应，输入参数较多，更接近于真实情况。过去20 多年来，Romanovsky 方程被广泛用于北极圈活动层厚度和多年冻土顶板温度的模拟、不同气候情景模式下活动层厚度的预测[85-86］、活动层加厚引发的潜在危险评估[87-88］，以及全球变暖背景下俄罗斯北极湿地温室气体的排放评估[89］。在青藏高原地区也有多次应用，如Pang等[90］模拟了青藏高原活动层厚度的空间分布；Luo等[91］得到黄河源区多年冻土顶板温度和活动层厚度的空间分布。从现有的文献来看，一般认为Ro⁃manovsky 方程的模拟精度要好于Stefan 方程。如Romanovsky 等[18］、Shiklomanov 等[92］在模拟阿拉斯加北坡多年冻土活动层厚度时得到Kudryavtsev 模型模拟精度远远大于Stefan方程。但也有的研究认为Kudryavtsev 方程所需要参数太多从而导致了模拟误差。如Yin 等[93］发现Steafn 方程在模拟五道梁多年冻土活动层厚度时的误差（＜7%）要小于Kudry⁃avtsev 模型（2.8%~27.4%）。从理论上讲，Kudry⁃avtsev 方程可以看作是Stefan 方程的一种全面改进，Romanovsky 等[18］对Kudryavtsev 方程与Stefan方程之间的转换方法进行了较为详细的介绍，本文不再赘述。

实际上，多年冻土的冻融结过程是土壤水分相变的过程，当水分由冰相转化为水相时，相对应的土壤热容量变小，热传导加快，从而使得冻结锋面的位置发生变化。在较小的日时间尺度上，含水量的影响也许可以忽略，但在时间较长的年尺度上具有至关重要的作用。因此，在一些含水量较高的多年冻土区，不仅要考虑总含水量，还要考虑到未冻水含量的年变化。未冻水不仅直接影响Stefan方程的模拟结果，而且通过影响导热系数进而再次影响模拟结果。因此，Stefan 方程在冻土模拟中未来考虑和改进的关键是未冻水，如何通过大尺度参量来确定土壤未冻水含量的变化是一个基本而又很重要的问题，这也是目前Stefan 方程最急需考虑和解决的重中之重。同时，多年冻土作为冰冻圈很重要的主体之一，其与其他圈层的相互作用越来越受到重视。一些学者试图在寒区陆面模式和水文模型中考虑加入多年冻土的影响，但是仅仅将多年冻土的影响作为其中一个单独的子模块。如何将多年冻土的诸多参数真正融入陆面模型和水文模型中是目前亟待解决的问题。

Stefan 方程参数少、形式简单、模拟效果可靠，是活动层厚度计算中运用最广泛、使用最方便的公式之一。但在国内的应用相对较少，现有的研究大多只是将此公式简单应用于模拟多年冻土活动层厚度的空间分布状况。随着青藏高原地区多年冻土变化研究的深入，Stefan 方程的应用必将日趋广泛。本文简要介绍了Stefan方程的推导背景和在冻土研究中的一些改进，希望起到抛砖引玉的作用，未来相关学者可在此基础上，对这一历史悠久的冻融过程模拟计算方法开展更深入地研究，使其在青藏高原多年冻土研究中得到更广泛地应用。