青藏高原热喀斯特湖演化及其对多年冻土的热影响模型计算研究

尹国安， 牛富俊， 林战举， 罗 京， 刘明浩

（中国科学院西北生态环境资源研究院冻土工程国家重点实验室，甘肃兰州730000）

0 引言

多年冻土作为冰冻圈要素之一，因其对气候变化具有高度敏感性以及重要的反馈作用而倍受关注。青藏高原被称为“世界第三极”，是山地多年冻土分布最广泛的区域［1-3］。近年来，随着全球气候变暖，青藏高原多年冻土温度升高，地下冰融化，活动层厚度增加，对高原生态环境、地貌类型、气候以及人类工程活动等产生了深刻的影响［2，4-5］。多年冻土区地下冰融化，引起地表土层沉陷或坍塌，从而形成多年冻土区特有的热喀斯特地貌［6］。其中，由于热融沉陷并产生积水形成的湖塘，被称为热融湖塘或热喀斯特湖（thermokarst lake）。热喀斯特湖广泛分布于北极地区［7］及青藏高原多年冻土区［8-9］。由于水体的热作用，热喀斯特湖的形成改变了地表能量平衡过程，加快了湖底深部以及湖岸周围多年冻土的融化速度［10-11］，尤其是当湖水深度超过了冬季最大冻结深度，湖底的未冻结水将不断对其下部多年冻土产生热影响，最终形成湖底融区（talik）。湖底融区将为地表水与地下水交换提供通道，向大气释放大量储存在地下冰中的二氧化碳和甲烷气体［10］。同时，热喀斯特湖在向周围扩张的过程中，不断侵蚀湖岸，形成不同规模的坍塌，对附近重要工程以及生态环境带来潜在威胁。

在北极多年冻土区，研究人员通过遥感技术［12］、现场观测［13］以及数值模型计算［14-15］等手段，分析了热喀斯特湖的形成和演化过程，加深了对热喀斯特湖在冻土区环境变化和全球气候系统中作用的认识。关于青藏高原多年冻土区热喀斯特湖的形成和演化过程的研究起步较晚，由于气候不断变暖及人类活动频率增加，近年来青藏高原热喀斯特湖不管是从数量还是面积都在逐年增加［16］。青藏高原大多数热喀斯特湖呈椭圆形，深度从几十厘米到几米不等［8］。湖温及湖冰观测数据表明，当湖水深度超过0.7 m 时，在冬季将不能回冻至湖底部［17］。热喀斯特湖的形态与夏季风向与风速，以及太阳辐射角度有很大关系［8］。青藏高原热喀斯特湖竖向发展及横向扩张受到多种因素的影响，变化复杂，由于系统观测资料的缺乏，其动态演化过程的机理目前仍然不清楚。Lin 等［17］通过10 年的地温观测，提出热喀斯特湖分为形成、动态发展、稳定和最终干涸等四个演化阶段，但这一过程需要上百年甚至千年［18］。罗京等［19］、You等［20］利用地质雷达技术，分别探讨了热喀斯特湖对周边多年冻土的热影响。Ling等［18］通过建立数值计算模型，分析了青藏高原热喀斯特湖形成以后对湖底及周围多年冻土的热影响，指出湖底开放融区形成需要约800 年，随后将处于稳定状态。令锋等［21］建立二维扩张模型，研究了热喀斯特湖横向扩张速率对湖下融区发展的影响。Li等［22］通过耦合水热过程模型，研究了热喀斯特湖在气候变暖背景下湖底融区的变化过程。这些研究从现场观测和模型模拟角度促进了对青藏高原热喀斯特湖演化过程的理解［23］。热喀斯特湖的形成和演化与地下冰状况密不可分［10，17］，因此地下冰含量及其热状况也是研究热喀斯特湖变化过程中普遍关心的重要问题，然而目前少有模型研究考虑了地下冰特别是过剩冰融沉在热喀斯特湖竖向发展过程中的重要作用。因此，本文试图在已有研究基础上，建立一个耦合气候—湖塘—冻土融沉相互作用过程的一维模型，探讨气候变化背景下青藏高原热喀斯特湖的形成及演化过程，加深热喀斯特湖对地气相互作用系统影响的理解。

本文的主要研究目标如下：①建立一维热喀斯特湖竖向变化过程数学模型，包含气候、湖水、湖冰，以及多年冻土融沉等关键水热相互作用过程；②以青藏高原北麓河盆地典型热喀斯特湖持续监测资料为模型输入信息，分析不同地下冰状况在气候变化背景下对热喀斯特湖形成和演化过程的影响。

1 研究区概况

本研究选取的热喀斯特湖位于北麓河盆地。该盆地位于青藏高原腹地（34.8°N、92.0°E），平均海拔约4 600 m［图1（a）～（b）］。地表植被类型主要为沼泽草甸、高寒草甸、高寒草原和荒漠化草原。地表土层主要为第四纪冲洪积粉土和细砂，地下深层为泥岩和砂岩［24］。根据北麓河气象站观测，近20年来，该地区年平均气温为-4.5 ℃［图1（c）］，降水量为357 mm。北麓河位于连续多年冻土区，地下冰含量丰富，最大体积含冰量达到了70%，平均过剩冰含量达到了19%［25］。多年冻土年平均温度介于-1.5～0 ℃之 间，活. 层 厚 度 介 于1.4～3.4 m 之间［24］。该地区热喀斯特湖发育，且近年来数量和面积持续增加［16］。北麓河地区拥有丰富的观测场地和大量关于气候变化、多年冻土地温、热喀斯特湖的连续观测数据，能够对本研究提供数据基础。再分析数据［26］及第五次耦合模式比较计划（CMIP5）模型气候模式HadGEM2-ES［27］在不同气候情景下（RCP2.6、4.5、8.5，详见2.2.2 节）预测表明，北麓河地区从2000 年开始至2100 年将会有较明显的增温过程，到21 世纪末，增温幅度最高将会达到6.9 ℃［图1（c）］。

2 数据与方法

2.1 现场观测数据

本研究的观测气象数据主要来自北麓河气象站［图1（b）］，该气象站记录了近20年来北麓河地区的主要气象要素变化，包括2.0 m 高度处的气温变化（芬兰Vaisala 公司HMP45C_L11）、太阳辐射（日本EKO 公司MS-102型辐射表）、降雨等。气象站自2001年12月开始观测，每30分钟记录一次数据。

图1 研究区基本情况Fig. 1 Basic information of the study area：permafrost distribution on the Qinghai-Tibet Plateau（QTP）［1］and location of Beiluhe（a），picture of Beiluhe field observations（b）and annual mean air temperature in history period［26］and projections in the 21th century［27］of Beiluhe area（c）

关于热喀斯特湖水热状态数据来源于北麓河A湖塘（BLH-A）的长期监测数据［28］。BLH-A 湖中央位置布设有温度传感器（热敏电阻电缆，精度为±0.05 ℃）以及湖冰厚度记录传感器（武汉大学WUUL-I 超声测距器，精度为1‰）［28］，自2010 年开始对湖底温度及湖冰状况进行监测，每3小时记录一次数据。同时，在BLH-A周边高寒草甸场地，一个深度为15 m的钻孔对多年冻土地温进行观测，0～10 m深度每隔0.5 m 布设一个热敏电阻（冻土工程国家重点实验室制作并校准，精度为±0.05 ℃），10～15 m深度每隔1 m 布设一个热敏电阻。数据记录设备为Campbell公司生产的CR3000数据采集系统，地温数据采集从2012年开始，每4小时采集一次。

本文最终计算了以上记录数据的日平均值，用于模型的检验和验证。

2.2 数值计算模型

图2为描述热喀斯特湖形成过程的示意图。根据已有研究［10，17］，热喀斯特湖的形成和演化主要经历三个阶段，即初始状态为富冰多年冻土［图2（a）］，在经历气候变暖过程后，富冰多年冻土温度升高，地下冰特别是过剩冰融化［图2（b）］并向地面排水，地表沉陷形成凹地，形成初始热喀斯特湖，随着地下冰不断融化，融水不断向上排出，湖塘深度加深［图2（c）］。本研究对该过程进行了数学建模。

图2 热喀斯特湖形成过程模型Fig. 2 A model depicting the formation process of thermokarst lake：initial state of ice-rich permafrost with excess ice（a），ice-rich permafrost is warming and the water is draining（b），subsided cells are added to lake depth（c）and schematic illustration of an ice-rich soil layer of thickness（Δd）which is composed of excess ice，pore ice and soil matrix（d）

2.2.1 控制方程

图2（a）中，多年冻土中的土壤温度（T）随时间（t）的变化由伴有相变的热传导方程描述［29］。

式中：C（z，T）为土壤容积热容量（J·m-3·℃-1），根据式（2）计算；k（z，T）为土壤导热系数（W·m-1·℃-1），根据式（3）计算；z为土层深度（m）；ρw为水的密度（kg·m-3）；Lw为水的相变潜热（3.34×105J·kg-1）［30］；θu为土壤未冻水含量，根据式（4）计算［31］。

式中：θj为各组分的体积含量；cj为各组分的比热容；下标j表示土壤中各组分，即土壤矿物质、水分和孔隙冰；a、b为与土壤性质有关的经验常数；θw为土壤融化时的体积含水量。根据已有研究［18，21，30］，模型中用到的各物理参数如表1所示。

表1 与土壤、水和冰有关的物理参数Table 1 Physical parameters related to soil，water and ice

图2（b）～（c）中，由于气温变暖，多年冻土厚层地下冰融化并假设全部排水至地表，同时，地表沉陷并形成湖塘，这个过程中，地表下沉总量可以用式（5）表示［32］。

式中：Δz为地表沉陷量（m）；Δd为富冰冻土层厚度（m）；θi为总含冰量（过剩冰与孔隙冰体积含量之和，即θi=θex+θp）；ønat为土壤孔隙度，饱和状态下ønat=θp。图2（d）为富冰多年冻土的组成示意图，式（5）的详细推导过程可参考文献［32］。

当湖塘形成以后，应用Liston 等［33］建立的一维水—冰—积雪模型模拟湖冰及湖水中的传热过程。在夏季，湖水温度主要受水的热传导及太阳辐射的影响，其控制方程如式（6）所示。

式中：Cw为水的容积热容量（J·m-3·℃-1）；Tw为湖水的温度（℃）；h为辐射在湖水中所传输的深度（m）；q为传入水中的辐射总量（W·m-2），与地表总辐射Q（W·m-2）有关。

式中：α为水面反照率（0.06）；η为消光系数（0.6 m-1）［33］。地表总太阳辐射量可以根据现场气象站获得（详见

2.2.2 节）。在冬季，湖塘水体表面形成冰层后，湖水中只考虑热传导。由于青藏高原热喀斯特湖较浅（＜3 m）［8］，同时已有研究表明，这些湖塘具有较好的混合度［11］，因此湖冰的形成与消融控制方程如式（8）所示。

式中：ρi为冰的密度（kg·m-3）；A为对流传导系数（0.56 W·℃·m-2）［33］；Tf为水的冻结温度（℃）；Tw为水温（℃），由式（6）得到；Ts为表面温度（℃）；hi为冰的厚度（m）。

2.2.2 模型求解与边界条件

模型耦合了冻土—湖水—湖冰的水热过程，热传导偏微分方程［式（1）和式（6）］采用有限差分方法求解，时间步长为1 d，求解时间为1 000 a（即365 000 d），空间步长设置为垂直剖面10 m 以上取0. 05 m，10 m 以下取1.0 m，垂直剖面总长度为1 000 m。求解过程中，每一个时间步长模型将会检测求解域的材料属性（水体、土壤、冰层），当检测到冰层融化，土壤层沉陷，模型自动更新求解域材料并匹配相应的求解方程。模型的所有求解过程采用Python 语言编程完成。

模型的上边界条件为北麓河地区的地表日平均温度与太阳辐射日均值。由于北麓河气象站观测时间较短，因此本文应用Yang等［26］研发的再分析资料重建北麓河地区长时间序列气温度数据（1990—2000 年），该数据与气象站资料具有较好的吻合度，并符合正弦函数变化规律（图3）。最终，模型上边界条件为符合相应正弦函数的气温与太阳辐射值。

图3 模型上边界历史气候数据的重建Fig. 3 Reconstruction of historical climate data for upper boundary of the model：daily air temperature（a）and solar radiation（b）（A sinusoidal fit to the data was used）

由北麓河气象站［图1（c）］观测数据表明，2000年以前年平均气温处于较稳定状态，年平均气温为-4.5 ℃，2000—2100年，在不同气候情景下，研究区气温呈现不同升温模式。因此，本研究设置上边界条件如下：在1 000 a 里，前500 a 气温处于较为稳定状态（图3）；在500～600 a，气温升高，本研究选取极端最高升温情景（HadGEM2-ES在RCP 8.5情景）下热喀斯特湖的演化过程，即100 a 间升温6.9 ℃［图1（c）］；600～1 000 a，气温再次处于平稳状态。下边界条件为1 000 m 处的地热梯度，设置为0.1 W·m-2［34-35］。

初始时假设热喀斯特湖未形成，模型的初始条件通过设置天然钻孔场地记录的不同深度土壤年平均值为各深度处的初始值，再通过上边界条件（0～500 a 气温）反复迭代计算直到温度场达到稳定状态，稳定状态下的地温场即为模型的初始条件，然后加入不同深度的热喀斯特湖，利用边界条件求解相关控制方程。

2.2.3 土壤物理参数与实验设计

根据北麓河热喀斯特湖周边的地质资料［24］，计算垂直区域可以分为3个子区域（表2），土壤类型及物理性质基于钻孔岩芯资料获取，即0～2 m为砂土，土壤体积含水量为0.30；2～10 m 为粉质黏土层，体积含水量为0.35，在2～10 m区间，含有富冰冻土层，最大体积含冰量达到了0.65，即过剩冰体积含量达到了0.3；10 m以下为砂岩层，体积含水量低于0.1。根据已有模拟研究［36］，与土壤未冻水含量有关的经验参数a、b取值如表2所示。

表2 土壤地层情况及相关模型参数Table 2 Soil texture and the related model parameters

为了模拟热喀斯特湖对近代气候变化的响应过程，本文设计了4 种模拟方案。考虑4 种不同深度的热喀斯特湖对气候变暖的响应，即湖塘深度分别为0.5 m、1.0 m、1.5 m 和2.0 m。这4 种深度的热喀斯特湖下伏多年冻土过剩冰含量相同，但分别位于不同深度（图4）：3.5～6.0 m、4.0~6.0 m、4.5～6.0 m和5.0～6.0 m。

图4 不同深度湖塘地下冰的初始状态Fig. 4 Initial states of underground ice for the lake：0.5 m deep（a），1.0 m deep（b），1.5 m deep（c）and 2.0 m deep（d）

2.3 模型验证

本文主要关注不同深度热喀斯特湖对气候变暖的响应过程，因此模型的验证主要从两个方面考虑：第一，验证模型对湖水及湖冰的模拟能力，主要是通过对比模型对湖塘底部温度的模拟以及湖冰的生长与消融过程；第二，验证模型对多年冻土温度的模拟能力。以上验证的数据来源于对北麓河典型热喀斯特湖BLH-A 的关键参数及其周边多年冻土温度的长期监测。验证的主要指标有确定系数（R2）、绝对误差（ME）以及均方根误差（RMSE）。

3 结果与分析

3.1 模型验证结果

图5展示了本文所发展的模型模拟结果与实测值的对比。在多年冻土耦合湖模型中，湖底部温度是一个非常重要的变量，它反映了水体与下伏沉积物之间的传热过程。本研究中，模型计算的湖底温度与测量值具有较好的吻合度［图5（a）］。在夏季，由于强风扰动作用，测量值要稍微高于模拟值，而在冬季，模型并未考虑冰面积雪对太阳辐射的遮挡作用，因此模拟值高于实测值。同时，这也使得模型对冰厚度的模拟值要高于实测值［图5（b）］。图5（c）～（d）为模型对多年冻土地温的模拟结果，模拟地温曲线与实测地温曲线基本重合，确定系数超过了0.9，绝对误差为0.2 ℃。通过以上对比可见，本文发展的模型对热喀斯特湖水热状态以及多年冻土地温变化具有较好的模拟能力。

图5 模型验证Fig. 5 Model validation：comparison between modeled and measured temperatures at bottom of the BLH-A Lake（a），comparison between modeled and measured lake ice thickness（b），comparison between modeled and measured ground temperatures at different depths on August 15，2015（c）and comparison between modeled and measured annual mean ground temperatures at different depths（d）

3.2 热喀斯特湖变化及其对多年冻土的影响

图6为模型计算的热喀斯特湖深度以及湖冰厚度随时间的变化规律。图6（a）表明，在模型计算的前500 a 里，气候稳定，只有0.5 m 深的热喀斯特湖保持稳定，地下冰没有融化产生融陷，其他深度的热喀斯特湖在开始的50 a里，深度不断增加，直至其下伏多年冻土中的过剩冰不再融化为止，最终，初始深度分别为1.0 m、1.5 m、2.0 m 的热喀斯特湖深度分别增加至1.4 m、1.8 m、2.2 m。同时图6（b）表明，初始深度为0.5 m、1.0 m 的热喀斯特湖在冬季湖水能够回冻至湖底，而初始深度为1.5 m、2.0 m的热喀斯特湖最大湖冰厚度为1.65 m，即冬季湖水不能回冻至湖底，而湖底年平均温度将超过0 ℃。因此，大于1.5 m 深度的热喀斯特湖在冬季仍然为热源，能够对下伏多年冻土产生热影响。实际上，当湖深超过0.5 m，湖底部就将会形成融区，如图6显示，在模型计算初期（100 a以内）融区厚度达到了20 m［图7（b）~（d）］，湖底年平均温度更是超过了4.0 ℃（图8）。在模型计算的500 a里，0.5 m热喀斯特湖下伏多年冻土各个深度处的年平均温度保持基本平稳状态，没有明显的增加［图7（a）、图8（a）］，而其他深度的热喀斯特湖底部多年冻土地温（小于25 m）在前100 a 里有一个缓慢的增温过程，同时在50 a 处，地温有一个跳跃（图8 虚线框所示），这是因为当年湖水未能回冻至湖底所致。

为了进一步查明在气候明显变暖情境下（21 世纪末期增温6.9 ℃）热喀斯特湖的变化过程，图7～图8绘出了在500～1 000 a所有深度湖塘及其下伏多年冻土所经历的增温过程。500～600 a 期间，0.5 m 湖塘下部过剩冰融化，湖塘急剧加深直至1.0 m，并在随后的时间里保持不变［图6（a）］。其他深度湖塘由于过剩冰在0～500 a 期间已经融化，因此，湖塘深度保持不变。最大湖冰厚度在500～600 a 期间变化明显。0.5 m 深度湖塘最大湖冰厚度增加至1.0 m，即冬季湖塘仍然能够回冻至湖底。其他深度湖塘最大湖冰厚度逐渐减小，直至1.1 m。如图7 所示，不同深度湖塘底部年平均地温在500～600 a 期间明显增加。0.5 m 深度湖塘底部年平均地温在随后的400 a 里逐渐增加，底部形成融区，但多年冻土并未完全融化。1.0 m 湖塘在700 a 时，底部多年冻土完全融化，形成贯穿融区，而其他深度湖塘底部多年冻土在600 a后将会完全融化消失。

图6 湖深与最大湖冰厚度的模拟结果Fig. 6 Simulation results of lake depth（a）and maximum lake ice thickness（b）

图7 不同深度湖塘在模拟期间的年平均地温剖面Fig. 7 Annual mean ground temperature profiles during simulation period for the lake：0.5 m deep（a），1.0 m deep（b），1.5 m deep（c）and 2.0 m deep（d）

图8 不同深度湖塘在模拟期间的年平均地温变化Fig. 8 Variation of annual mean ground temperatures during simulation period for the lake：0.5 m deep（a），1.0 m deep（b），1.5 m deep（c）and 2.0 m deep（d）

4 讨论

本文以青藏高原北麓河多年冻土区热喀斯特湖为背景，研究了4 种不同深度热喀斯特湖对气候变化的响应过程以及对下伏多年冻土的热影响。热喀斯特湖的形成及发展受到多重因素的影响且变化复杂［17］，以往的多数研究往往只关注了热喀斯特湖的横向扩张［15，22，37-38］，对热喀斯特湖的竖向发展研究较少，特别是热喀斯特湖变化与地下冰的关系仍然不清楚。本文基于一维热传导数值计算并耦合了地下过剩冰融沉过程，发展了热喀斯特湖变化模型，量化了湖深、湖冰深度与地下冰融化的关系。

模拟结果与北麓河盆地典型热喀斯特湖观测值具有较好的吻合度。该模型将有利于理解青藏高原富冰多年冻土区地貌演化过程及在气候变化背景下的变化规律。该模型并未考虑积雪对湖冰的影响，同时，该模型只考虑了一维热传导，忽略了侧向对流以及周边冻土对湖的热影响。因此，未来将耦合更多的物理过程，比如加入地表能量平衡过程以及侧向传热过程。

目前，关于青藏高原热喀斯特湖演化规律研究重点关注特定深度湖塘对多年冻土的水热影响［18，21-22］，简化了水体这个上边界条件（即简化为湖底温度），且并未考虑纯冰层的热状况及其融沉效应。本研究详细考虑了不同深度湖塘水体的热变化规律，同时也考虑了地下冰的融沉效应。本研究结果显示，在近代稳定气候背景下，2.0 m 深度湖塘底部的融区在500 a 时深度超过了40 m，600 a 后消失，这比同期其他模型［18］模拟值要高。模拟结果表明，当湖水深度超过1.0 m，湖底年平均温度要高于气温，其热源效应明显，因此，青藏高原的较深热喀斯特湖的热源效应要比以往预估更高，同时，伴随气候变化，热喀斯特湖深度不断加深，其热源效应也随之加强。

与极地地区相比［14，39］，青藏高原多年冻土具有较高的地温以及较薄的过剩冰层，因此热喀斯特湖的竖向发展受到相应的约束，但水体对下伏多年冻土的影响却更大。本文研究结果表明，0.5 m 热喀斯特湖在气候变暖过程中，下伏多年冻土将持续退化形成融区，而大于1.0 m 热喀斯特湖底部将会形成开放融区。因此，青藏高原地区热喀斯特湖底部多年冻土将会对气候变化极为敏感。

5 结论

（1）热喀斯特湖具有热源效应，但不同深度的热喀斯特湖对地下多年冻土热影响差异较大。0.5 m深热喀斯特湖将会处于稳定状态，而深度超过1.0 m的热喀斯特湖即使在目前稳定气候状态下（0～500 a期间），地下冰持续融化，100 a内底部将会形成超过20 m 深的融区。深度超过1.5 m，冬季湖冰将不能回冻至湖底。

（2）在气候变暖影响下（500～600 a 期间），热喀斯特湖热源效应加剧，即使0.5 m 深浅湖底部多年冻土发生快速退化，较深湖塘底部融区快速增加，在100 a 内（500～600 a 期间）多年冻土融化消失，形成开放融区。