青藏高原热喀斯特湖与多年冻土的相互作用

柯贤敏,欧安锋,王 玮,牛富俊,高泽永

(1. 长安大学水利与环境学院,陕西 西安 710054；2. 长安大学旱区地下水文与生态效应教育部重点实验室,陕西 西安 710054；3. 中国科学院西北生态环境资源研究院冻土工程国家重点实验室,甘肃 兰州 730000)

青藏高原是中国多年冻土的代表,影响着中国乃至全球的气候变化[1]。近年来,在人为活动和气候变暖的影响下,青藏高原多年冻土逐步退化[2- 3]。多年冻土退化最直接的体现是冻土活动层厚度增大和热喀斯特湖数量及面积的增大[4- 5]。热喀斯特湖作为重要的热源,对区域水文循环、碳循环以及生态环境[6- 7]具有重要影响。目前,国外开展的研究主要集中于热喀斯特湖对地貌、区域环境及碳循环的影响[8- 9]。青藏高原热喀斯特湖的研究主要集中于水热过程、分布特征、湖水水文地球化学特征及其对工程稳定性影响方面[10- 13],关于热喀斯特湖与多年冻土相互作用的研究相对较少。深入了解热喀斯特湖与多年冻土间的相互作用,有助于更好地理解多年冻土区的生态水文过程。

长期以来,学者们多通过遥感影像和现场监测手段分析热喀斯特湖的发展过程及其与冻土间的水热过程[14]。然而,热喀斯特湖与周围多年冻土的长期相互作用十分复杂,涉及热传导、水分迁移、冰水相变等过程,短期的观测数据无法揭示和预测热喀斯特湖演化过程及其与冻土间的相互作用[15]。因此,数值模拟技术成为研究冻土水热变化、热喀斯特湖发育的重要手段。国外针对多年冻土区的水热耦合模型主要有2个研究方向:一个是只考虑冻土内部的水热变化[16- 17],另一个是考虑在大气- 植被- 土壤整个系统中的水热变化,再加上降水和融雪等过程[18- 19]。目前为止,针对青藏高原热喀斯特湖开展的模拟研究相对较少。林战举等[20]通过构建热传导数值模型,预测了热喀斯特湖对多年冻土的影响,结果表明湖的温度对湖周围冻土影响显著；尹国安等[21]通过构建一维热传导模型,模拟了不同深度的湖在气候变暖背景下对多年冻土的影响；令锋等[22]利用热传导模型计算了热喀斯特湖不同扩张速率对多年冻土退化的影响,结果表明湖的扩张速率越大,对冻土退化影响越大；纪舒婷[23]利用数值模型模拟了冻土退化对湖泊储量的影响,指出冻土退化增加了地下水入湖量；杨振等[24]利用移动网格技术构建了热喀斯特湖演化的热传导数值模型,分析了热喀斯特湖深度和半径对冻土的影响,并得出恒定的热喀斯特湖边界会高估对湖底冻土的影响、低估对湖附近冻土的影响。

上述研究多是基于热传导理论开展的,忽略了热对流过程对冻土退化的贡献,这些方法只适用于地下水流速很小的环境。青藏高原多年冻土分布格局和动态变化与地下水密切相关[25],仅考虑热传导过程对冻土的影响是不全面的。为此,本文以北麓河盆地典型热喀斯特湖区域为研究对象,结合野外监测与数值模拟手段,对热喀斯特湖与多年冻土间的相互作用展开研究,分析热喀斯特湖的水均衡变化特征,探讨地质环境和气候变暖对热喀斯特湖和生态环境发展方向的影响。

1 研究区与研究方法

1.1 研究区概况

北麓河研究区位于青藏高原北麓河盆地(图1),属于青藏高原亚寒带半干旱气候。研究区年平均气温、降水量、蒸发量分别为-4 ℃、274.70 mm和1 421.40 mm。降水通常以雨、雪和冰雹的形式发生于6—8月。研究区内地层起伏不大,自上而下依次为第四系全新统冲洪积层、上第三系湖相沉积层。全区分布厚度约为30～100 m的多年冻土,活动层厚度一般为1.60～3.40 m。地下水的排泄方式主要有蒸发、以泉的形式溢流于地表和向热喀斯特湖及冻土层下水补给等。本文选取的热喀斯特湖(BLH- B)位于北麓河盆地内部,面积约为7 hm2,最大深度为2 m(图1),湖周围下垫面类型为高寒草地。

图1 BLH- B湖温度监测点布设Fig.1 Distribution of temperature monitoring sites in BLH- B

1.2 野外试验和数据

采用MRDT605型温度测量仪对湖中和湖附近区域温度进行监测。湖东北、西南和东南侧各设置了3个监测点(图1),监测深度为地面之下0.10 m、0.25 m、0.50 m、1.00 m,监测频率为4 h/次；在湖中布设了3个温度监测点(图1),用于监测湖底之下的地温。湖底温度和湖周围地温(地面之下0.10 m)分别用于模型湖底边界和地表边界的赋值。此外,湖周围地温监测数据也用于对计算模型进行验证。

1.3 水- 冰- 热耦合模型

1.3.1 主要控制方程

SUTRA能够解决饱和非饱和条件下的溶质运移和地下水流动等问题。然而,SUTRA原始程序并不能对多年冻土区地下水的冻结和融化进行刻画。为此,McKenzie等[26]对SUTRA程序进行了改进,得到能够模拟饱和条件下地下水冻结与融化过程的SUTRA- ICE程序,该程序同时考虑了热传导和热对流过程。SUTRA- ICE程序的地下水流动控制方程为

(1)

SUTRA- ICE程序中能量传递的实质是冰、孔隙介质和液体中的能量平衡,其方程可以表示为

(2)

1.3.2 模型结构和边界条件

根据模拟区实际地形坡度(约0.87°),利用水- 冰- 热耦合计算程序,采用柱坐标系建立二维剖面模型。模型横向由外流域边界到湖泊中心(X=0)的500 m径向距离确定,纵向由深度确定,模型径向截面如图2所示。模型在垂向上共划分了3大层,自上而下地层渗透率依次为5×10-13m2、1×10-13m2和1×10-14m2。模型中间层和底层共划分了49小层,单层厚度为10 m。为了更好地展示活动层的变化,将顶层进一步划分为100层,单层厚度为0.10 m。此外,为了分析热喀斯特湖周围冻土的变化情况,将湖及湖岸附近20 m区域进一步划分为35列,每列宽2 m；其他区域划分为43列,列宽为10 m。模型厚度为H=2πX。模型的顶部设定为空气热力学边界,其温度为大气温度,压强为0 Pa；模型底部设定为地热通量边界(0.06 W/m2)；两侧为零通量边界,不与外界发生水热交换。模型西北角设定1个宽50 m、深2 m的湖,利用温度和压强边界来刻画湖泊。

图2 模型结构、网格剖分及边界条件概化Fig.2 Figure showing the structure,grid,and boundary condition of model

2 结果及分析

2.1 模型验证

构建多年平均气温(-3 ℃)和地热通量共同作用下的稳定流模型,得到稳定的多年冻土结构。稳定流模型的温度和压强计算结果作为非稳定流模型的初始条件。基于稳定流模型,在模型西北角增加1个深2 m、宽50 m的湖泊,根据2018—2019年热喀斯特湖和地表温度监测数据分别对模型地表边界和湖底边界进行赋值,并利用温度监测数据对计算模型进行验证。实测温度与计算温度组成的坐标点均位于1∶1直线附近,相关系数(R2)为0.97～0.99,均方根误差(ERMS)为0.48～1.31 ℃(图3),计算结果与实测结果拟合较好。因此,该模型可用于模拟热喀斯特湖的排泄量及其与多年冻土间的相互作用。模型参数如表1所示。

图3 实测温度与计算温度的对比Fig.3 Comparison of the measured and calculated temperature

表1 模型最终参数表

2.2 多年冻土及活动层变化特征

为了揭示热喀斯特湖长期作用下多年冻土和活动层的变化特征,在稳定流模型的基础上构建了长达1 000 a的非稳定流数值模型。该模型设置了12 000个时间步长,每个时间步长代表1个月。地表和湖底边界温度根据多年平均温度进行循环(图4)。当地层冰饱和度小于0.05时(无颜色覆盖区域),认为地层未发生冻结；模型中的剩余液态水饱和度设定为0.05,即冻结的冻土中也有少量水分存在,故代表多年冻土存在的冰饱和度范围为0.05～0.95。冰饱和度的增加(由红色区域向蓝色区域变化)代表了地层冻结程度增强(图5)。

图4 热喀斯特湖和大气温度变化Fig.4 Temperature variation of thermokarst lake and air

湖水持续向多年冻土传递热量,引起湖周围冻土逐步退化,多年冻土面积逐渐减小(表2)。第200 a,湖底冻土已经开始退化但尚未融穿；而热喀斯特湖持续作用400 a,湖底形成了贯穿融区。贯穿融区形成前,地下冰限制了地下水的流动,地下水流速缓慢,热传导占主导作用，此时热量主要在垂向上进行传递,垂向上的冻土退化速度快于水平方向；贯穿融区形成后,热对流占主导作用,热量伴随着地下水的流动进行传递。在湖附近区域,多年冻土底界面因不断受到地下水流的热侵蚀而抬升，同时,地下水流速和湖的排泄量有所增大,而湖获得的补给量减小。地表的排泄量及其获得的补给量在贯穿融区形成前后基本稳定,在模拟期末出现增长，由于热喀斯特湖的影响范围有限,距离湖较远的区域多年冻土保持稳定；贯穿融区形成后,地下水流速增大,携带热量的地下水对多年冻土的热侵蚀作用增强,活动层厚度增大了0.23 m。整个模拟时段内,活动层地下水保持近水平流动。在多年冻土的影响下,冻土层下水、冻土层上水、冻土层间水以及湖水间的水力联系微弱,冻土层下水在静水压力作用下自东南向西北流动,流速自下而上逐渐增大。对比第200 a和第1 000 a计算结果可以得出,贯穿融区形成后,水流通道恢复,冻土层下水的流向发生了改变,各水体间的水力联系增强,形成了新且稳定的地下水循环系统。

图5 不同时期多年冻土、流速矢量及活动层分布计算结果Fig.5 Distribution of permafrost,active layer,and velocity vector in different periods

表2 不同时段模型计算结果统计表

图6截取了模型第400～420 a的计算结果,来反映湖边观测孔处(X=60 m)的活动层冻融循环过程。地表长期处于升温、降温循环过程中,故形成了厚度为3.35 m的活动层。11月至次年5月,地表温度低于0 ℃,活动层处于冻结状态；而6—10月,地表温度高于0 ℃,活动层处于融化状态。整个模拟时段内年际间地表温度未升高,故活动层厚度保持稳定。地表热量向深部传递的过程中逐渐被消耗,故活动层深部处于融化状态的时间比浅部短,而处于冻结状态的时间比浅部有所延长。此外,随着深度的增加,深部多年冻土对地表温度的响应也表现出一定的滞后性。

图6 湖边观测孔处活动层厚度计算结果Fig.6 Active layer thickness of lakeside observation well

2.3 热喀斯特湖水均衡计算结果分析

图7 热喀斯特湖水均衡计算结果Fig.7 Results of thermokarst lake water balance

多年冻土和热喀斯特湖的水热过程势必会引起地表边界和湖底边界的水均衡发生变化。湖底边界获得的补给量(进入湖底边界的水量)及其自身的排泄量(湖底边界流出的水量)决定了湖底边界的水均衡状况。如图7(a)所示,热喀斯特湖获得的补给量与其自身的排泄量同步变化。随着湖底温度的升高,流入湖底的冻土层上水水量增加,湖与冻土层下水间的水头差减小,湖的排泄量减小，为了维持湖水位保持不变,湖获得的补给量也同步减小；反之,湖水温度降低过程中,湖的排泄量及其获得的补给量均增加。在1个温度周期内,湖表现为负均衡,即获得的补给量小于自身的排泄量。湖向周围冻土传递热量需要一定的时间,从而导致湖的补给量和排泄量滞后于其温度的变化。热喀斯特湖在整个模拟时段内也表现为负均衡,即一直处于排水状态(图7(b))。模型第285～388 a计算结果表明:随着冻土含冰量的降低,地下水流动通道恢复,热喀斯特湖的排泄量逐渐增大,其获得的补给量减小；热喀斯特湖持续作用388 a时,湖底形成贯穿融区,热喀斯特湖的净排泄量达到最大；贯穿融区形成后,冻土层上水、冻土层下水以及湖水间的水力联系增强,形成新的地下水循环系统；随后,由于湖底边界的压强是恒定的,故热喀斯特湖的排泄量开始减小并逐渐恢复至初始状态。

3 讨 论

3.1 不同地质环境下的热喀斯特湖发展方向

地层渗透率影响着地下水的径流速度和湖的排泄量。图8(a)显示了不同地层渗透率(方案1—方案4的地层渗透率分别为1×10-12m2、1×10-13m2、1×10-14m2和1×10-40m2)地层环境下热喀斯特湖净排泄量的变化特征。随着地层渗透率的增大,地下水流动受到的阻滞作用减小,热喀斯特湖的排泄量逐渐增大。地层渗透率很小时(1×10-40m2),地下水近似处于静止状态,湖的排泄量接近于0 kg/s。湖底形成贯穿融区时,湖水将持续向深层地下水排泄,湖的排泄量有所增大,引起湖水进一步流失。当湖获得的补给量得不到保证时,热喀斯特湖将逐渐排干,最终形成新的冻土。北极地区地层岩性以砂卵石为主,地层渗透率大,热喀斯特湖底冻土融穿后,湖的排泄量增大,而其获得的补给量有限,故一些热喀斯特湖逐渐干涸消失[27]。本次模拟的青藏高原热喀斯特湖区域,地层岩性以砂岩和泥岩为主,地层渗透性弱,湖的排泄量相对较小。虽然湖底形成贯穿融区后,湖的排泄量有所增大,但增大程度有限。此外,气候变暖加速冻土退化,湖获得来自地下水的补给量增大,降水量也能得到保证,故热喀斯特湖逐渐扩张[28]。

图8 不同渗透率地层环境下湖和地表排泄量计算结果Fig.8 Results of the lake and ground- surface discharge in different permeability environments

3.2 不同地质环境下的生态环境效应

天然条件下,降水量、蒸发量、地表径流量与地下水径流量处于动态均衡状态。当降水入渗补给量大于维持地表水位所需补给量时,地下水完全饱和处于蓄满产流状态,形成地表径流；而当降水入渗补给量小于维持地表水位所需补给量时,潜水面将会下降。不同的地质环境下,冻土活动层厚度及地表的排泄量(地表边界流出的水量)也有所不同。随着地层渗透率的增大,地表的排泄量也逐步增大(图8(b)),从而需要获取更多的补给量才能维持地表水位；当地层渗透率较小时(1×10-14m2和1×10-40m2),地表的排泄量接近于0 kg/s，此时,很小的补给量便能维持地表水位。

地表排泄量会直接影响地表生态环境的发展。地质环境对生态环境的影响主要表现在以下2个方面:① 渗透率较大时,地表排泄量大,用于维持地表水位的补给量也越大。当地表获得的补给量小于该值时,地下水位将不断下降，植被生长与冻土水热过程的平衡受到干扰,本就十分脆弱的植被生态系统被破坏,不仅使地表侵蚀容忍量骤降,而且生态环境将逐渐向荒漠化的趋势发展。② 渗透率较小时,地表排泄量较小,用于维持地表水位的补给量也较小。当降水入渗补给量大于该值时,地下完全饱和,处于蓄满产流状态并形成地表径流,生态环境将向着沼泽湿地方向发展。对于青藏高原而言,地表朝着沼泽湿地的趋势发展,有利于发挥高原湿地的生态蓄水、水源补给和气候调节等重要生态功能,使生态系统良性循环。

3.3 气候变暖条件下热喀斯特湖与多年冻土相互作用

在气候变暖的影响下,多年冻土不断退化,对冻土水文过程及高原生态环境产生了重要影响[29]。为探究气候变暖背景下热喀斯特湖与多年冻土的相互作用,在前述模型基础上构建了湖和地表升温方案。升温方案中,湖和地表边界温度在多年平均温度(图4)的基础上,以0.03 ℃/a升温直到167 a后保持不变,温度整体升高5 ℃。由图9(a)和图9(b)可以看出,第100 a,模型范围内存在高含冰量冻土；而到第400 a,仅模型中部残留了小范围的低含冰量(冰饱和度小于0.06)冻土,表明地表和湖温度的升高加速了多年冻土的退化。对比第100 a和第400 a计算结果可知,冻土退化使水流通道得到恢复,地下水流速增大；同时,地下水流向也发生了变化。模拟中后期,深部冻土完全退化,地表仅在冬季发生冻结,此时形成了季节冻土(图9(c))。模型第50～70 a间,在地表升温的作用下,活动层厚度呈现逐渐增加的变化趋势,由3.81 m增加至4.50 m(图9(c))。由于受到湖的侧向热侵蚀,活动层之下的冻土逐步退化,融区上界面不断上移；最终,深部冻土全部退化,地表形成季节冻土。气候变暖加速冻土的退化,将加速释放储存在冻土中的碳,从而加剧全球气候变暖,进一步加速冻土的退化。此外,在气候变暖的背景下,冻土退化可能导致潜水面下降,青藏高原生态环境更容易向荒漠化方向发展。

图9 升温方案模拟结果Fig.9 Simulated results of the temperature rise scheme

4 结 论

本文以青藏高原北麓河盆地典型热喀斯特湖区域为研究对象,建立了考虑热传导和热对流过程的“水- 冰- 热”耦合数值模型,并对其进行了验证。模拟了热喀斯特湖与多年冻土的相互作用关系及热喀斯特湖的排泄量,并对不同地质和气候变暖条件下的热喀斯特湖发育方向及生态效应进行了探讨。主要结论如下:

(1) 热喀斯特湖持续向多年冻土传递热量,导致多年冻土逐步退化。贯穿融区形成前后,热传导和热对流在多年冻土退化过程中分别占主导作用。在地表温度作用下,形成的活动层厚度为3.35 m。

(2) 贯穿融区形成后,冻土层下水、冻土层上水以及湖水间的水力联系增强,形成新的地下水循环系统。整个模拟时段内,热喀斯特湖表现为负均衡,其排泄量在285～388 a间显著增加。

(3) 地层渗透率大时,地表和湖的排泄量大,热喀斯特湖易于疏干,生态环境将逐渐向荒漠化方向发展；反之,地表和湖的排泄量小,热喀斯特湖呈现扩张趋势,生态环境向沼泽湿地方向发展。

(4) 气候变暖加速了冻土退化和热喀斯特湖的发展,多年冻土逐渐转变为季节冻土,严重影响寒区生态水文过程。