滑坡涌浪诱发冰湖溃决灾害链过程分析与模拟

刘 威，杨宗佶*，游 勇，黄 栋，聂 勇，张广泽，徐正宣，王 栋

(1.中国科学院 水利部 成都山地灾害与环境研究所，四川 成都 610041；2.中铁二院工程集团有限责任公司，四川 成都 610031)

滑坡诱发冰湖溃决灾害链是高山冰川作用区常见的灾害类型之一，具有突发性强、波及范围广、发生频率集中等特点。中国青藏高原区域地处太平洋板块与印度洋板块交接处，地质运动活跃，冰湖众多，为该类灾害链的发生提供了有利条件。例如，2013年7月，西藏那曲地区嘉黎县忠玉乡发生冰湖溃决洪水灾害，导致238户1 160人不同程度受灾，经济损失高达2.7 亿元。近年来，全球气候变暖背景下冰川作用扰动强烈，导致中国高山冰川区滑坡诱发冰湖溃决灾害链的发生概率提高。同时，中国重大工程穿越帕隆藏布流域等冰湖分布密集区域，滑坡诱发冰湖溃决灾害链风险较大，对线路造成长期威胁。因此，开展滑坡诱发冰湖溃决灾害链的机理分析及危险性研究是防灾减灾的迫切需求。

滑坡诱发冰湖溃决灾害链的危害可表现在：一是，滑坡及诱发涌浪的直接危害，受滑坡和涌浪的冲击范围与运移速度等因素影响；二是，冰湖溃决洪水的间接危害，受冰湖容量、冰湖溃决速率、地形等因素影响。由此可见，灾害链的危害具有阶段性且与各阶段动力特征密切相关。因此，要定量评估滑坡诱发冰湖溃决灾害链危险性，就需要深入研究其演化过程，量化分析其动力特征，从而制定合理有效的防灾减灾措施。近年来，有关灾害链的研究正逐步从以模式分类、机理描述等为主的定性阶段向以物理模型、数值模拟等为主的定量阶段发展。目前，针对滑坡诱发冰湖溃决灾害链中各子阶段的演化过程已存在较多研究。在滑坡运动及涌浪传播方面，Zhou等以三峡红岩子滑坡涌浪灾害为例，利用流固耦合方法模拟了其成灾过程，定量给出了滑坡诱发涌浪的最大高度以及衰减特征。Bai等基于深度平均理论，提出非静力双层流模型来描述滑坡入水现象，模型较好地再现了涌浪传播的色散关系，在数值求解方面，利用交错有限差分方法离散模型方程，利用凯勒尔盒子法在垂直方向上重构非静力源项，并通过试验数据验证了模型及算法的可行性。在冰湖溃决及洪水传播方面，黄金辉等通过水槽试验研究了涌浪诱发冰湖溃决的动力过程，探讨冰（雪）崩塌入湖诱发涌浪溃坝的发展过程及溃决流量特征，分析了不同涌浪规模对溃坝形成机制的影响。严佳宏通过振动台试验研究涌浪对冰碛堰塞湖溃决的影响，结果表明，小型终碛湖在峰值加速度高于 0.11

g

的地震波作用下可能因为涌浪发生溃决。Amicarelli等提出一种可考虑液相、固相和固液混合相3种组分的混合物模型，结合光滑粒子流动力学方法模拟了几组溃坝洪水侵蚀试验，较好地反映了下游河床物质随洪水运移特征。值得注意的是，尽管上述方法均可较好反映出滑坡、涌浪、洪水等灾害的动力特征，但灾害链演化过程包括多个阶段，且各阶段物理机理存在显著差异，难以通过单一的物理模型来描述整体灾害链演化过程与阶段间交互影响。针对该问题，Liu等采用耦合灾害链各子阶段物理模型的方法来描述灾害链的整体演化过程，并成功模拟了西藏易贡滑坡—堰塞湖—溃决洪水灾害链成灾过程。然而，由于冰湖所处地区通常地形险峻，现场勘查十分困难，导致该类灾害链的机理研究仍然较为缺乏，且缺乏行之有效的物理模型及定量评估方法。

为此，本文以滑坡—涌浪—冰湖溃决—洪水灾害链为研究对象，通过确定灾害链阶段转化的影响因子，衔接灾害链各阶段物理模型，进而建立可描述灾害链完整演化过程的耦合模型；以四川省甘孜藏族自治州雅拉乡木格措冰湖潜在灾害链为例，在模拟不同复杂情景条件下灾害链演化过程的同时实现危险性定量评价。

1 木格措冰湖

木格措冰湖位于四川省甘孜藏族自治州康定市雅拉乡康定情歌风景区内，图1为所在区域影像图。

图1 木格措冰湖区域平面图Fig. 1 Diagram of Mugecuo glacier lake area

图1中，冰湖面积为1.85 km，长2.9 km，最宽处1.1 km，平均湖面宽约630 m，库容达82.2×10m，平均湖深约为45 m。1986—2018年间的遥感影像监测流域湖泊变化资料分析表明：冰湖总体保持稳定，存在微小的年内和年际变化，导致湖泊面积出现小幅度的波动。湖泊的状态与流域的特征有密切的关系，木格措流域的冰川已消失，存在季节性的降雪，上游的冰湖的水源补给和输出处于平衡状态，冰湖水位总体保持稳定。

根据木格措冰湖流域地貌条件与现场勘查，确定木格措冰湖在东北角出水口为一古滑坡堰塞坝（图2）。由于水流改造，目前堰塞坝中上部物质已被侵蚀，溢流口的河床顺直和平缓（约5%），河床堆积砾石粗化明显。堆积体平面形态呈不规则舌形，堆积体后壁边界清晰，呈圈椅状地貌，边界总长度约0.95 km，顶部高程约4 000 m，高差220 m，最大坡度约50°（平均坡度20°），长530 m，宽437 m，总面积约0.15 km，平均厚度约40 m，体积规模约6×10m。现场调查未发现堰塞坝有整体变形迹象，处于稳定状态。

图2 木格措冰湖滑坡堰塞坝平面图Fig. 2 Diagram of the landslide barrier dam in Mugecuo glacier lake

2 木格措滑坡堆积体

木格措滑坡堆积体位于木格措冰湖北岸（图1），平面形态呈半圆形，堆积体后壁边界为山脊，呈近似圈椅状地貌，边界总长度约1.5 km，后壁海拔高度在4 200～4 600 m范围内，滑坡主滑方向205°，堆积体前缘以木格措冰湖湖面为界，湖水面高程约3 780 m，滑坡堆积体顶部高程约4 466 m，高差686 m，平均坡度30°（最小坡度约5°，最大坡度约50°），纵向长度约1 224 m，最大宽度1 350 m，面积1.42 km。结合现场调查类比临近滑坡体的基岩和滑坡后壁基岩产状，确定滑坡堆积体前缘厚度约80～100 m，中部厚度约40～50 m,平均厚度约50 m，体积约7×10m。

滑体堆积体组成以堆积层和崩积层为主。坡体前缘部分被表层岩质崩积层覆盖，坡体中部以粉质粘土夹碎块石及下部碎石组成，碎石含量在5%～30%，且分布不均，部分位置有大块石出露，粒径8～20 m，结构破碎，在一定范围内形成岩屑堆。碎石土厚度一般为10～50 m，密实度为稍密-密实，土质不均匀，碎石母岩主要为喜山期细粒二长花岗岩(ηγN)，填充粉质黏土、粉土、沙土，其上植物发育良好，根系发达。崩积层主要分布于滑坡体上部及前缘表层，滑坡后壁的断壁下部多松散风化基岩，在长期重力作用下形成高位的岩屑堆，成分以细粒二长花岗岩为主，现场调查及遥感解译发现，表层块石粒径约8～10 m，后缘处高原植被发育不良，偶生有簇状灌木，以高山杜鹃为典型代表。滑坡体前缘部位，厚度0.5 m以上，以大块石、砾石为主，细砾物质较少，表面呈灰色，灰白色，成分也以二长花岗岩为主，现场调查表层块石块径一般0.10～2.00 m，细砾物质填充在块石之间。现场调查并未发现堆积体有整体变形迹象，整体稳定。

3 滑坡诱发冰湖溃决灾害链物理模型

考虑滑坡—涌浪—冰湖溃决—洪水灾害链演化特点，将其划分为4个子阶段：滑坡运动、涌浪传播、冰湖溃决、洪水传播，分别构建滑坡运动物理模型、滑坡入水与涌浪传播物理模型以及溃决洪水演进物理模型，在考虑灾害链动力演化过程中形态改变、时空变化的关键控制要素条件下，搭建滑坡—涌浪—冰湖溃决—洪水灾害链的耦合物理模型。

3.1 滑坡运动物理模型

在滑坡体3维运动方程的基础上，利用深度平均理论对纳维斯托克斯3维方程进行了简化，即:假设滑坡体在运动过程中其

z

方向上的变量保持一致，从而将复杂的3维运动方程简化为2维运动方程，同时进一步考虑了滑坡内部侧压力变化对其动力过程的影响，模型方程如下：式（1）～（3）中：

t

为时间；

h

为滑坡体厚度；

u

、

v

分别为滑坡体沿

x

和

y

方向的速度；µ为基底摩擦系数；

g

、

g

、

g

分别为沿实际地形下坡面、纵坡面和垂直坡面向下的重力加速度分量；ρ = (1-

n

)ρ+

n

ρ为滑坡体密度，其中，ρ为固相颗粒密度，ρ为液相密度，

n

为滑坡体孔隙度；

k

和

k

分别为滑坡体沿

x

和

y

方向的侧向土压力系数，

k

和

k

均为土体内摩擦角

φ

和基底摩擦角

φ

的函数。

3.2 滑坡入水及涌浪传播物理模型

滑坡入水及涌浪传播过程的影响因素有很多，包括滑坡体特征、河流特征、山区地形，边界条件等。基于此，从最基本的质量与动量守恒方程入手，假设滑坡体与堰塞湖水体均为不可压缩，基于深度平均理论，以双层流理论为基本框架，构建针对性的物理动力学模型。在考虑动力学边界条件及深度平均理论简化后，用于描述滑坡入水及涌浪传播过程的双层流模型方程组可表示如下：

式（4）～（9）中，ρ为冰湖水密度，ρ为滑坡体密度，

u

和

v

分别为涌浪沿

x

和

y

方向的速度，

u

和

v

分别为滑坡体沿

x

和

y

方向的速度，

h

和

h

为冰湖水及滑坡体厚度，

g

为重力加速度，τ(

z

)和τ(

z

)为冰湖水与滑坡体的摩擦项，τ(

z

)和τ(

z

)为滑坡体与湖床的摩擦项，

z

为冰湖水与滑坡体的接触面，

z

为滑坡体与湖床的接触面。

3.3 冰湖溃决及洪水传播物理模型

漫顶溃决是堰塞坝主要的溃决方式之一，其溃决过程往往是涌浪在坝顶造成一处缺口，进而产生溃坝洪水。针对涌浪造成的侵蚀现象，相应的计算公式如下：

式中，

E

为溃口侵蚀速率，α和β为侵蚀经验系数，

U

为侵蚀临界速度。洪水传播过程采用浅水波方程，可表示如下：

式（11）～（13）中，

h

为洪水高度，

R

为降雨强度，

u

、

v

分别为洪水为沿

x

和

y

方向上的流速，

S

和

S

分别为洪水沿

x

和

y

方向上所受到的摩擦阻力。

3.4 模型耦合与算法实现

灾害链演化机理的复杂性在于需同时考虑各阶段的条件特征，如滑坡体特性、地形特征、初始及边界条件等，以及不同阶段之间的衔接耦合。针对这一问题，筛选各阶段的输入输出参数，确保模型方程通过参数实现信息交互，具体选取参数及相应计算方法如下：滑坡初始厚度分布、滑坡组成物质特征、地形等参数作为初始条件输入滑坡运动物理模型，得出实时滑坡厚度分布、速度分布及滑动距离等参数，再通过滑动距离判定滑坡是否能进入冰湖；将滑坡临近冰湖的厚度分布与速度分布、冰湖水深分布作为初始条件带入滑坡入水及涌浪传播物理模型，得出滑坡堆积形态、涌浪传播高度、传播速度等参数；将冰湖水深分布、涌浪传播高度、传播速度作为初始条件带入冰湖溃决及洪水传播物理模型，得出实时堰塞湖水量、溃口形态、洪水水深分布、洪水流速等参数。最终，可通过所获的参数信息来实时定量分析灾害链演化过程及动力特征。上述模型方程均利用有限体积法、Well-balanced地形重构技术和?MUSCL界面重构技术进行求解。

3.5 冰湖溃决洪峰流量经验模型

根据堰塞湖应急处置技术导则，采用洪峰流量的经验计算公式对数值模拟结果进行复核，公式如下：式中：

Q

为坝体下游

x

处洪峰流量；

Q

为坝址处洪峰流量；

W

为坝体下泄总水量；

L

为下游断面至坝址距离；

v

为河道洪水期断面最大平均流速，本文中取值2～5 m/s；

K

系数，取值1.10～1.25。

4 木格措滑坡诱发冰湖溃决灾害链过程模拟

假设滑坡堆积体在极端条件下（如地震）局部复活诱发冰湖溃决灾害链，以此模拟灾害链形成过程，采用计算参数见表1。选择鱼司通和三道桥两个断面开展冰湖溃决洪水灾害危害分析。

表1 灾害链过程模拟计算参数
Tab. 1 Parameters used in the simulation of disaster chain

参数 取值滑坡初始厚度/m 20滑坡初始范围/m2 0.27×106滑坡初始体积/m3 5.4×106冰湖初始厚度/m 45冰湖初始范围/m2 1.84×106冰湖初始体积/m3 8.22×107滑坡内摩擦角/(°) 24滑坡底摩擦角/(°) 20曼宁摩擦系数 0.015

滑坡运动与涌浪传播过程的计算结果如图3所示。

图3 不同运动时刻滑坡厚度及涌浪水深分布Fig. 3 Depth distributions of landslide and indced surge at different times

从图3可以看出：滑坡体启动后发生迅速变形，顺滑面下滑，滑坡体运动过程中的最大运动速度为5.9 m/s，且处于滑体前端，至滑入湖底后，在湖水和基底的摩擦阻力作用下逐渐产生堆积，随着后续坡体物质的不断涌入和叠加，滑坡体停止运动后的最大堆积厚度为40.1 m。当滑坡前端入湖面后激发涌浪，涌浪快速向四周传播，诱发的涌浪运动到溃口处的最大运动速度为1.05 m/s，最大高度为5 m。在以上涌浪条件下，计算溃坝下切最大深度为0.82 m，溃口处最大洪峰流量为778 m/s，溃决洪水运动到下游雅拉河鱼司通断面和三道桥断面的峰值流量分别为471.0和385.3 m/s，最大水深分别为1.84 和 0.90 m。同时，根据堰塞湖应急处置技术则取

Q

为778 m/s，

L

为16.6 km，

W

为6.93×10m，

v

为2 m/s，

K

为1.1，可得三道桥断面处洪峰流量为421.2 m/s，经验公式计算结果与数值模拟的结果基本一致。可见，本文所提出的滑坡诱发冰湖溃决灾害链物理模型可较好反映出灾害链各阶段的动力特征。

康定及木格措地区降水年内分配不均，降水主要集中在夏季，5～9月降雨量达621 mm，占年降水量的77%，存在暴雨作用下形成山洪灾害的可能性。如1995年6月15日、7月3日和7月7日,康定城区连续3次遭受以山洪为主的多种山地灾害袭击，城区被淹和遭泥沙淤埋，造成对外交通、通讯中断，直接经济损失约5亿元。因此，在上述研究的基础上，开展考虑暴雨作用下形成山洪条件下的木格措冰湖溃决灾害链形成过程模拟。选取叠加康定地区百年一遇（发生概率为1%，雨强27.4 mm/h）概率降雨下形成的洪水，其余计算参数和表1一致，计算结果如图4所示，在暴雨作用下，坡面产生汇流形成坡面流，坡面流先汇入小支沟处，再由小支沟汇入较大支沟，最终汇入主沟道，增加沟道流量，冰湖溃决发生后，洪水与坡面流存在汇流，导致洪水规模扩大，向下游传播时不断汇集沟道坡面水流。此时，鱼司通断面最大计算流量为1 544.6 m/s，最大水深为4.2 m；三道桥断面最大流量为1 333.6 m/s，最大水深为3.1 m。

图4 暴雨条件下冰湖溃决洪水传播水深分布Fig. 4 Distribution of the lake outbrust depth under the conditon of heavy rainfall

此外，通过现场勘查发现，鱼司通和三道桥上游雅拉河有一定数量的支沟存在泥石流活动，如鱼斯通沟可能存在泥石流堵河风险，因此，考虑在最不利条件下，即“滑坡诱发冰湖溃决灾害链+泥石流堰塞湖级联溃决+暴雨-山洪”组合的可能性，对此情景开展过程模拟。选取叠加康定地区百年一遇（发生概率为1%，雨强27.4 mm/h）概率降雨下形成的洪水，其余计算参数和表1一致。现场勘查表明，鱼司通泥石流堵河形成堰塞坝时，泥石流坝体平均高度为5～6 m。计算结果如图5所示，当上游滑坡诱发冰湖溃决产生的洪水传播至泥石流坝处时，由于坝体拦截作用洪水开始蓄积，随着后续流量的不断汇入，所形成的堰塞湖最大水深8.6 m，面积为1.18×10m，库容为4.2×10m。假设新形成的堰塞湖全部溃决，同时叠加雅拉河流域百年一遇暴雨条件下的洪水，当堰塞湖最大水位超过堰塞坝顶高后，坝体溃决的洪水流量放大效应以及叠加暴雨引起的沟道坡面汇流，共同作用下导致洪水流量陡增（图6）。洪水到鱼司通断面最大流量为2 019.0 m/s，最大水深5.2 m；三道桥断面最大流量为1 871.8 m/s，最大水深为3.5 m。

图5 溃坝洪水传播至泥石流坝前蓄积水深分布Fig. 5 Distribution of the lake outbrust depth at the front of the debris flow dam

图6 考虑暴雨洪水及级联溃决条件下洪水水深分布Fig. 6 Distribution of flood depth under the conditons of heavy rainfall and cascade dam break

5 结 论

滑坡诱发冰湖溃决灾害链的定量分析与风险评价对减灾有重要支撑作用。本文针对该类灾害链的过程特点，采用多物理模型耦合的思路，以四川省甘孜藏族自治州雅拉乡木格措冰湖潜在灾害链为案例，定量预测了木格措滑坡运动、滑坡入水及涌浪传播、冰湖溃决与洪水传播等关键过程，并进一步讨论了暴雨洪水和级联溃决两种不同情景条件对灾害链危害的影响。预测结果表明：考虑地震条件下规模约5.4×10m滑坡失稳，前端入湖面后激发涌浪，涌浪运动到溃口处的最大运动速度为1.05 m/s，最大高度为5 m，鱼司通和三道桥断面峰值流量分别是471.0和385.3 m/s，最大水深分别是1.84和 0.90 m；当叠加暴雨-洪水和上下游级联溃决洪水条件时，下游洪峰流量及水深呈显著增加趋势。以上研究成果可为滑坡诱发冰湖溃决灾害链应急管理与防灾减灾提供关键技术支撑。

需要说明的是，由于缺乏滑坡和堰塞坝勘察资料和湖底地形等详细数据，相关土体物理力学参数只能通过类比法参考工程地质手册确定，滑坡及堰塞坝的剖面只能通过物探结合现场调查和工程地质类比法确定，滑坡规模和冰湖库容主要通过现场调查结合经验公式计算，同时地形图的精度也有限，这些方面的不足，一定程度上制约了计算结果的精度。本研究通过野外现场调查结合全动力过程数值模拟，探索了“地质认识+数值计算”相结合的链生山地灾害的正演预测方法，其预测结果精度还有待未来进一步结合实际进行检验和修正。