工程科学与技术   2022, Vol. 54 Issue (2): 41-48
滑坡涌浪诱发冰湖溃决灾害链过程分析与模拟
刘威1, 杨宗佶1, 游勇1, 黄栋1, 聂勇1, 张广泽2, 徐正宣2, 王栋2     
1. 中国科学院 水利部 成都山地灾害与环境研究所,四川 成都 610041;
2. 中铁二院工程集团有限责任公司,四川 成都 610031
基金项目: 中国科学院STS项目(KFJ–STS–QYZD–172);中国科学院重点部署项目(KFZD–SW–425);第二次青藏高原综合科学考察研究(2019QZKK0904);国家自然科学基金青年项目(41907241)
摘要: 近年来,滑坡诱发冰湖溃决灾害链发生频率日趋增加,对下游居民及工程设施安全造成严重威胁。为定量描述滑坡诱发冰湖溃决灾害链危害,通过系统分析该灾害链的演化过程及阶段特征,将其划分为滑坡运动、涌浪传播、冰湖溃决和洪水传播4个子阶段;分别研究各子阶段演化的物理机理及关键因子,建立相应的数学模型与计算方法,如滑坡运动采用Savage–Hutter模型、涌浪及洪水传播采用浅水波方程等;通过确定各子阶段之间的衔接因子,如滑坡运动特征、涌浪传播特征、坝体溃决特征等,实现各子阶段的过程耦合与数据传递,建立滑坡—涌浪—冰湖溃决—洪水灾害链的全过程物理模型,最终完成该灾害链的过程模拟与危险分析。为验证所建立模型与计算方法的可行性,以四川省甘孜藏族自治州雅拉乡木格措冰湖潜在灾害链为案例,结合现场勘查与卫星影像数据,开展不同情境,如暴雨或上下游级联溃决等条件下灾害链的过程模拟与危险分析,且将模拟结果与经验公式计算数据进行对比。结果表明,本文所提方法可较好的完成滑坡诱发冰湖溃决灾害链过程模拟,模型模拟数据与经验公式计算数据较为吻合,且考虑暴雨或上下游级联溃决条件时,下游洪水流量呈显著增加趋势。本文研究成果提供了一种新的滑坡诱发冰湖溃决灾害链研究思路,对该灾害链的防灾减灾技术具有支撑作用。
关键词: 灾害链    滑坡涌浪    冰湖溃决    物理模型    数值模拟    过程分析    
Dynamic Analysis and Numerical Simulation of Disaster Chain from Landslide,Surge Wave to Glacier Lake Outburst
LIU Wei1, YANG Zongji1, YOU Yong1, HUANG Dong1, NIE Yong1, ZHANG Guangze2, XU Zhengxuan2, WANG Dong2     
1. Inst. of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu 610041, China;
2. China Railway Eryuan Eng. Group Co. Ltd., Chengdu 610031, China
Abstract: Recently the frequency of disaster chain of landslide, surge wave and glacier lake outburst increases, which poses a serious threat to the safety of downstream residents and engineering facilities. In order to describe quantitatively the hazards of this disaster chain, its evolution process and corresponding characteristics are analyzed and further divided into four stages including landslide motion, surge wave propagation, glacier lake outburst and flood propagation. The physical mechanisms and key factors of these stages are investigated and based on this the corresponding model and numerical scheme are also established. For example, the Savage–Hutter model and shallow water model are applied to describe landslide motion and surge/flood propagation, respectively. By determining the factors that connect each stage, the coupling and data transmission of these stages are achieved and therefore a comprehensive model of this disaster chain is presented for process simulation and hazard analysis. In order to verify the feasibility of the presented model and numeric scheme, a potential disaster chain in Mugecuo glacier lake near Yala village, Ganzi state of Sichuan province is studied. Based on the data of field investigation and satellite imagery, the process of disaster chain under several scenarios such as heavy rainfall and cascade dam break are performed and the simulated results are compared with that of empirical formula. These results show that the presented method is quite effective to simulate the disaster chain of landslide, surge wave and glacier lake outburst and the calculated data agree with that of empirical formula. It also indicates that the flood discharge can be increased significantly due to the effects of heavy rainfall or cascade dam break. This research provides a new idea for studying the disaster chain of landslide, surge and glacier lake outburst and has some consultative value to prevent and mitigate this type of disaster chain.
Key words: disaster chain    landslide-induced surge    glacier lake outburst    physical model    numerical simulation    process analysis    

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

滑坡诱发冰湖溃决灾害链的危害可表现在:一是,滑坡及诱发涌浪的直接危害,受滑坡和涌浪的冲击范围与运移速度等因素影响[11-12];二是,冰湖溃决洪水的间接危害,受冰湖容量、冰湖溃决速率、地形等因素影响[4,13]。由此可见,灾害链的危害具有阶段性且与各阶段动力特征密切相关。因此,要定量评估滑坡诱发冰湖溃决灾害链危险性,就需要深入研究其演化过程,量化分析其动力特征,从而制定合理有效的防灾减灾措施。近年来,有关灾害链的研究正逐步从以模式分类[14]、机理描述[15]等为主的定性阶段向以物理模型[16]、数值模拟[17]等为主的定量阶段发展。目前,针对滑坡诱发冰湖溃决灾害链中各子阶段的演化过程已存在较多研究。在滑坡运动及涌浪传播方面,Zhou等 [18]以三峡红岩子滑坡涌浪灾害为例,利用流固耦合方法模拟了其成灾过程,定量给出了滑坡诱发涌浪的最大高度以及衰减特征。Bai等[19]基于深度平均理论,提出非静力双层流模型来描述滑坡入水现象,模型较好地再现了涌浪传播的色散关系,在数值求解方面,利用交错有限差分方法离散模型方程,利用凯勒尔盒子法在垂直方向上重构非静力源项,并通过试验数据验证了模型及算法的可行性。在冰湖溃决及洪水传播方面,黄金辉等[20]通过水槽试验研究了涌浪诱发冰湖溃决的动力过程,探讨冰(雪)崩塌入湖诱发涌浪溃坝的发展过程及溃决流量特征,分析了不同涌浪规模对溃坝形成机制的影响。严佳宏[21]通过振动台试验研究涌浪对冰碛堰塞湖溃决的影响,结果表明,小型终碛湖在峰值加速度高于 0.11g 的地震波作用下可能因为涌浪发生溃决。Amicarelli等[22]提出一种可考虑液相、固相和固液混合相3种组分的混合物模型,结合光滑粒子流动力学方法模拟了几组溃坝洪水侵蚀试验,较好地反映了下游河床物质随洪水运移特征。值得注意的是,尽管上述方法均可较好反映出滑坡、涌浪、洪水等灾害的动力特征,但灾害链演化过程包括多个阶段,且各阶段物理机理存在显著差异,难以通过单一的物理模型来描述整体灾害链演化过程与阶段间交互影响。针对该问题,Liu等[23]采用耦合灾害链各子阶段物理模型的方法来描述灾害链的整体演化过程,并成功模拟了西藏易贡滑坡—堰塞湖—溃决洪水灾害链成灾过程。然而,由于冰湖所处地区通常地形险峻,现场勘查十分困难,导致该类灾害链的机理研究仍然较为缺乏,且缺乏行之有效的物理模型及定量评估方法。

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

1 木格措冰湖

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

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

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

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

图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 km2。结合现场调查类比临近滑坡体的基岩和滑坡后壁基岩产状,确定滑坡堆积体前缘厚度约80~100 m,中部厚度约40~50 m,平均厚度约50 m,体积约7×107 m3

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

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

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

3.1 滑坡运动物理模型

在滑坡体3维运动方程的基础上,利用深度平均理论对纳维斯托克斯3维方程进行了简化[23-24],即:假设滑坡体在运动过程中其z方向上的变量保持一致,从而将复杂的3维运动方程简化为2维运动方程,同时进一步考虑了滑坡内部侧压力变化对其动力过程的影响[25],模型方程如下:

$ \frac{{\partial h}}{{\partial t}} + \frac{{\partial \left( {hu} \right)}}{{\partial x}} + \frac{{\partial \left( {hv} \right)}}{{\partial y}} = 0 $ (1)
$ \begin{aligned}[b] \frac{{\partial \left( {hu} \right)}}{{\partial t}} &+ \frac{{\partial \left( {h{u^2}} \right)}}{{\partial x}} + \frac{{\partial \left( {huv} \right)}}{{\partial y}} =\hfill \\& {g_x}h - \frac{1}{2}{k_x}{g_\textit{z}}h\frac{{\partial h}}{{\partial x}} - \frac{u}{{\sqrt {{u^2} + {v^2}} }}\mu h\left[ {{g_\textit{z}}\left( {1 - \frac{{{\rho _{\rm{f}}}}}{\rho }} \right)} \right] \hfill \\ \end{aligned} $ (2)
$ \begin{aligned}[b] \dfrac{{\partial \left( {hv} \right)}}{{\partial t}} &+ \dfrac{{\partial \left( {huv} \right)}}{{\partial x}} + \dfrac{{\partial \left( {h{v^2}} \right)}}{{\partial y}} =\hfill \\& {g_y}h - \dfrac{1}{2}{k_y}{g_\textit{z}}h\dfrac{{\partial h}}{{\partial y}} - \dfrac{v}{{\sqrt {{u^2} + {v^2}} }}\mu h\left[ {{g_\textit{z}}\left( {1 - \dfrac{{{\rho _{\rm{f}}}}}{\rho }} \right)} \right] \hfill \\[-9pt] \end{aligned} $ (3)

式(1)~(3)中:t为时间;h为滑坡体厚度;u、v分别为滑坡体沿xy方向的速度;μ为基底摩擦系数;gxgygz分别为沿实际地形下坡面、纵坡面和垂直坡面向下的重力加速度分量;ρ = (1–n)ρf+s为滑坡体密度,其中,ρs为固相颗粒密度,ρf为液相密度,n为滑坡体孔隙度;kxky分别为滑坡体沿xy方向的侧向土压力系数,kxky均为土体内摩擦角 $\varphi $ int和基底摩擦角 $\varphi $ bed的函数。

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

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

$ \frac{{\partial {h_1}}}{{\partial t}} + \frac{{\partial {h_1}{u_1}}}{{\partial x}} + \frac{{\partial {h_1}{v_1}}}{{\partial y}} = 0 $ (4)
$ \frac{{\partial {h_2}}}{{\partial t}} + \frac{{\partial {h_2}{u_2}}}{{\partial x}} + \frac{{\partial {h_2}{v_2}}}{{\partial y}} = 0 $ (5)
$ \begin{aligned}[b] \frac{{\partial \left( {{h_1}{u_1}} \right)}}{{\partial t}} + &\frac{{\partial \left( {{h_1}u_1^2 + \dfrac{1}{2}gh_1^2} \right)}}{{\partial x}} + \frac{{\partial \left( {{h_1}u_1{v_1}} \right)}}{{\partial y}} =\hfill \\& - g{h_1}\frac{{\partial \left( {{\textit{z}_{\rm{b}}} + {h_2}} \right)}}{{\partial x}} + \frac{{{\tau _{x\textit{z}}}\left( {{\textit{z}_{\rm{m}}}} \right)}}{{{\rho _1}}} \hfill \\[-9pt] \end{aligned} $ (6)
$ \begin{aligned}[b] \frac{{\partial \left( {{h_1}{v_1}} \right)}}{{\partial t}} &+ \frac{{\partial \left( {{h_1}{u_1}{v_1}} \right)}}{{\partial x}} + \dfrac{{\partial \left( {{h_1}v_1^2 + \dfrac{1}{2}gh_1^2} \right)}}{{\partial y}} = \hfill \\& - g{h_1}\frac{{\partial \left( {{\textit{z}_{\rm{b}}} + {h_2}} \right)}}{{\partial y}} + \frac{{{\tau _{y\textit{z}}}\left( {{\textit{z}_{\rm{m}}}} \right)}}{{{\rho _1}}} \hfill \\[-9pt] \end{aligned} $ (7)
$ \begin{aligned}[b] \frac{{\partial \left( {{h_2}{u_2}} \right)}}{{\partial t}} +& \frac{{\partial \left( {{h_2}u_2^2 + \dfrac{1}{2}{k_{x}}gh_2^2} \right)}}{{\partial x}} + \frac{{\partial \left( {{h_2}{u_2}{v_2}} \right)}}{{\partial y}} = \hfill \\& - {k_{x}}rg{h_2}\frac{{\partial {h_1}}}{{\partial x}} - {k_{x}}{h_2}g\frac{{\partial {{\text{z}}_{\rm{b}}}}}{{\partial x}} + \frac{{{\tau _{x\textit{z}}}\left( {{\textit{z}_{\text{b}}}} \right) - {\tau _{x{\textit{z}}}}\left( {{\textit{z}_{\rm{m}}}} \right)}}{{{\rho _2}}} -\hfill \\& {k_{x}}r\sin \;{\varphi _{{{\rm{int}}} }}g{h_2}\frac{{\partial {h_1}}}{{\partial y}} - {k_{x}}\sin\; {\varphi _{{{\rm{int}}} }}g{h_2}\frac{{\partial {h_2}}}{{\partial y}} -\hfill \\& {k_x}{\text{sin}}\;{\varphi _{{{\rm{int}}} }}g{h_2}\frac{{\partial {\textit{z}_{\rm{b}}}}}{{\partial y}} \hfill \\[-17pt] \end{aligned} $ (8)
$ \begin{aligned}[b] \frac{{\partial \left( {{h_2}{v_2}} \right)}}{{\partial t}} &+ \frac{{\partial \left( {{h_2}{u_2}{v_2}} \right)}}{{\partial x}} + \frac{{\partial \left( {{h_2}v_2^2 + \dfrac{1}{2}{k_{y}}gh_2^2} \right)}}{{\partial y}} =\hfill \\& - {k_{{y}}}rg{h_2}\frac{{\partial {h_1}}}{{\partial y}} + {k_{{y}}}g{h_2}\frac{{\partial {{\text{z}}_{\rm{b}}}}}{{\partial y}} + \frac{{{\tau _{y\textit{z}}}\left( {{\textit{z}_{\text{b}}}} \right) - {\tau _{y\textit{z}}}\left( {{\textit{z}_{\rm{m}}}} \right)}}{{{\rho _2}}} -\hfill \\& {k_{y}}r\sin\; {\varphi _{{{\rm{int}}} }}g{h_2}\frac{{\partial {h_1}}}{{\partial x}} - {k_{y}}\sin\; {\varphi _{{{\rm{int}}} }}g{h_2}\frac{{\partial {h_2}}}{{\partial x}} -\hfill \\& {k_{y}}{\text{sin}}\;{\varphi _{{{\rm{int}}} }}g{h_2}\frac{{\partial {{\textit{z}}_{\rm{b}}}}}{{\partial x}} \hfill \\[-15pt] \end{aligned} $ (9)

式(4)~(9)中,ρ1为冰湖水密度,ρ2为滑坡体密度,u1v1分别为涌浪沿xy方向的速度,u2v2分别为滑坡体沿xy方向的速度,h1h2为冰湖水及滑坡体厚度,g为重力加速度,τxz(zm)和τyz(zm)为冰湖水与滑坡体的摩擦项,τxz(zb)和τyz(zb)为滑坡体与湖床的摩擦项,zm为冰湖水与滑坡体的接触面,zb为滑坡体与湖床的接触面。

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

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

$ E = \left\{ \begin{gathered} 0, | {{{({u}_1^2 + {v}_1^2)^{\frac{1}{2}}}}} | \lt {U_{\rm{l}}}; \hfill \\ \beta {h_1}{( ({{{{u}_1^2 + {v}_1^2}) \mathord{/ {\vphantom {{{\mathbf{u}}_1^2} {U_{\rm{l}}^2 - 1}}} } {U_{\rm{l}}^2 - 1}}} )^\alpha }, | {{{({u}_1^2 + {v}_1^2)^{\frac{1}{2}}}}} | \ge {U_{\rm{l}}} \hfill \\ \end{gathered} \right. $ (10)

式中,E为溃口侵蚀速率,αβ为侵蚀经验系数,Ul为侵蚀临界速度。洪水传播过程采用浅水波方程,可表示如下:

$ \frac{{\partial {h_{\rm{f}}}}}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {{h_{\rm{f}}}{u_{\rm{f}}}} \right) + \frac{\partial }{{\partial y}}\left( {{h_{\rm{f}}}{v_{\rm{f}}}} \right) = R $ (11)
$ \frac{\partial }{{\partial t}}\left( {{h_{\rm{f}}}{u_{\rm{f}}}} \right) + \frac{\partial }{{\partial x}}\left( {{h_{\rm{f}}}u_{\rm{f}}^2 + \frac{1}{2}gh_{\rm{f}}^2} \right) + \frac{\partial }{{\partial y}}\left( {{h_{\rm{f}}}{u_{\rm{f}}}{v_{\rm{f}}}} \right) = - g{h_{\rm{f}}}\left( {\frac{{\partial {\textit{z}_{\rm{b}}}}}{{\partial x}} + {S_{{\rm{f}}x}}} \right) $ (12)
$ \frac{\partial }{{\partial t}}\left( {{h_{\rm{f}}}{v_{\rm{f}}}} \right) + \frac{\partial }{{\partial x}}\left( {{h_{\rm{f}}}{u_{\rm{f}}}{v_{\rm{f}}}} \right) + \frac{\partial }{{\partial y}}\left( {{h_{\rm{f}}}v_{\rm{f}}^2 + \frac{1}{2}gh_{\rm{f}}^2} \right) = - g{h_{\rm{f}}}\left( {\frac{{\partial {\textit{z}_{\rm{b}}}}}{{\partial y}} + {S_{{\rm{f}}y}}} \right) $ (13)

式(11)~(13)中,hf为洪水高度,R为降雨强度,ufvf分别为洪水为沿xy方向上的流速,SfxSfy分别为洪水沿xy方向上所受到的摩擦阻力。

3.4 模型耦合与算法实现

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

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

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

$ {Q_{x{\rm{m}}}} = \frac{{W{Q_{\rm{m}}}}}{{W + \dfrac{{{Q_{\rm{m}}}L}}{{vK}}}} $ (14)

式中:Qxm为坝体下游x处洪峰流量;Qm为坝址处洪峰流量;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

滑坡运动与涌浪传播过程的计算结果如图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 m3/s,溃决洪水运动到下游雅拉河鱼司通断面和三道桥断面的峰值流量分别为471.0和385.3 m3/s,最大水深分别为1.84 和 0.90 m。同时,根据堰塞湖应急处置技术则取Qm为778 m3/s,L为16.6 km,W为6.93×106 m3v为2 m/s,K为1.1,可得三道桥断面处洪峰流量为421.2 m3/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 m3/s,最大水深为4.2 m;三道桥断面最大流量为1 333.6 m3/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×105 m2,库容为4.2×105 m3。假设新形成的堰塞湖全部溃决,同时叠加雅拉河流域百年一遇暴雨条件下的洪水,当堰塞湖最大水位超过堰塞坝顶高后,坝体溃决的洪水流量放大效应以及叠加暴雨引起的沟道坡面汇流,共同作用下导致洪水流量陡增(图6)。洪水到鱼司通断面最大流量为2 019.0 m3/s,最大水深5.2 m;三道桥断面最大流量为1 871.8 m3/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×106 m3滑坡失稳,前端入湖面后激发涌浪,涌浪运动到溃口处的最大运动速度为1.05 m/s,最大高度为5 m,鱼司通和三道桥断面峰值流量分别是471.0和385.3 m3/s,最大水深分别是1.84和 0.90 m;当叠加暴雨–洪水和上下游级联溃决洪水条件时,下游洪峰流量及水深呈显著增加趋势。以上研究成果可为滑坡诱发冰湖溃决灾害链应急管理与防灾减灾提供关键技术支撑。

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

参考文献
[1]
Cui Peng,Dang Chao,Cheng Zunlan,et al. Debris flows resulting from glacial-lake outburst floods in Tibet,China[J]. Physical Geography, 2010, 31(6): 508-527. DOI:10.2747/0272-3646.31.6.508
[2]
Benn D I,Benn T,Hands K,et al. Response of debris-covered glaciers in the mount everest region to recent warming,and implications for outburst flood hazards[J]. Earth-Science Reviews, 2012, 114(1/2): 156-174. DOI:10.1016/j.earscirev.2012.03.008
[3]
Xu Daoming,Feng Qinghua. Dangerous glacial lake and outburst features in Xizang Himalayas[J]. Acta Geographica Sinica, 1989, 44(3): 343-352. [徐道明,冯清华. 西藏喜马拉雅山区危险冰湖及其溃决特征[J]. 地理学报, 1989, 44(3): 343-352. DOI:10.3321/j.issn:0375-5444.1989.03.011]
[4]
Liu Jingjing,Cheng Zunlan,Li Yong,et al. Characteristics of glacier-lake breaks in Tibet[J]. Journal of Catastrophology, 2008, 23(1): 55-60. [刘晶晶,程尊兰,李泳,等. 西藏冰湖溃决主要特征[J]. 灾害学, 2008, 23(1): 55-60. DOI:10.3969/j.issn.1000-811X.2008.01.013]
[5]
Cheng Zunlan,Zhu Pingyi,Dang Chao,et al. Hazards of debris flow due to glacier-lake outburst in southeastern Tibet[J]. Journal of Glaciology and Geocryology, 2008, 30(6): 954-959. [程尊兰,朱平一,党超,等. 藏东南冰湖溃决泥石流灾害及其发展趋势[J]. 冰川冻土, 2008, 30(6): 954-959.]
[6]
Wang Xin,Liu Shiyin,Guo Wanqin,et al. Hazard assessment of moraine-dammed lake outburst floods in the Himalayas,China[J]. Acta Geographica Sinica, 2009, 64(7): 782-790. [王欣,刘时银,郭万钦,等. 我国喜马拉雅山区冰碛湖溃决危险性评价[J]. 地理学报, 2009, 64(7): 782-790. DOI:10.3321/j.issn:0375-5444.2009.07.002]
[7]
Sun Meiping,Liu Shiyin,Yao Xiaojun,et al. The cause and potential hazard of glacial lake outburst flood occurred on July 5,2013 in Jiali County,Tibet[J]. Journal of Glaciology and Geocryology, 2014, 36(1): 158-165. [孙美平,刘时银,姚晓军,等. 2013年西藏嘉黎县“7.5”冰湖溃决洪水成因及潜在危害[J]. 冰川冻土, 2014, 36(1): 158-165.]
[8]
Liu Jingjing,Tang Chuan,Cheng Zunlan,et al. Imapct of temperature on glacier-lake outbursts in Tibet[J]. Journal of Jilin University (Earth Science Edition), 2011, 41(4): 1121-1129. [刘晶晶,唐川,程尊兰,等. 气温对西藏冰湖溃决事件的影响[J]. 吉林大学学报(地球科学版), 2011, 41(4): 1121-1129.]
[9]
Cheng Zunlan,Tian Jinchang,Zhang Zhengbo,et al. Debris flow induced by glacial-lake break in Southeast Tibet[J]. Earth Science Frontiers, 2009, 16(6): 207-214. [程尊兰,田金昌,张正波,等. 藏东南冰湖溃决泥石流形成的气候因素与发展趋势[J]. 地学前缘, 2009, 16(6): 207-214.]
[10]
Yao Lingkan,Qiu Yanling,Wei Yongxing. Challenges in construction of railway and highway from Sichuan to Tibet through eastern margin of Tibetan platea[J]. Journal of Southwest Jiaotong University, 2012, 47(5): 710-734. [姚令侃,邱燕玲,魏永幸. 青藏高原东缘进藏高等级道路面临的挑战[J]. 西南交通大学学报, 2012, 47(5): 710-734.]
[11]
Huang Bolin,Yin Yueping,Chen Xiaoting,et al. Experimental modeling of tsunamis generated by subaerial landslides:two case studies of the Three Gorges Reservoir,China[J]. Environmental earth Sciences, 2014, 71(9): 3813-3825. DOI:10.1007/s12665-013-2765-5
[12]
Ai Hongzhou,Yao Lingkan,Zhou Yiliang. Laboratory investigations of earthquake-and landslide-induced composite surges[J]. Journal of Mountain Science, 2017, 14(8): 1537-1549. DOI:10.1007/s11629-016-4339-y
[13]
Sadeghiamirshahidi M,Vitton S J. Tropical storm-induced landslide-dammed lakes and debris flow hazards at Ocotepeque,Western Honduras[J]. Landslides, 2019, 16(1): 55-64. DOI:10.1007/s10346-018-1067-0
[14]
Cui Yun,Kong Jiming,Tian Shujun,et al.The critical role for heavy rainfall in the evolution of the mountain hazards chains[J].Journal of Mountain Science,2011,29(1):87–94.
崔云,孔纪名,田述军,等.强降雨在山地灾害链成灾演化中的关键控制作用[J].山地学报.2011,29(1):87–94.
[15]
Huang Yidan.Then nonlinear response characteristics of transmission system in mountain disaster chain[D].Chengdu:Soutwest Jiaotong University,2014.
黄艺丹.山地灾害链承传系统非线性响应特性研究[D].成都:西南交通大学,2014.
[16]
Liu Wei,He Siming. Dynamic simulation of a mountain disaster chain:landslides,barrier lakes,and outburst floods[J]. Natural Hazards, 2018, 90(2): 757-775. DOI:10.1007/s11069-017-3073-2
[17]
Xu Wenjie,Chen Zuyu,He Bingshun,et al. Research on river-blocking mechanism of Xiaojiaqiao landslide disasters of chain effects[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(5): 933-942. [徐文杰,陈祖煜,何秉顺,等. 肖家桥滑坡堵江机制及灾害链效应研究[J]. 岩石力学与工程学报, 2010, 29(5): 933-942.]
[18]
Zhou Jiawen,Xu Fuguang,Yang Xingguo,et al. Comprehensive analyses of the initiation and landslide-generated wave processes of the 24 June 2015 Hongyanzi landslide at the Three Gorges Reservoir,China[J]. Landslides, 2016, 13(3): 589-601. DOI:10.1007/s10346-016-0704-8
[19]
Bai Yefei,Cheung K F. Depth-integrated free-surface flow with a two-layer non-hydrostatic formulation[J]. International Journal for Numerical Methods in Fluids, 2012, 69(2): 411-429. DOI:10.1002/fld.2566
[20]
Huang Jinhui,Liu Jiankang,Cheng Zunlan,et al. An experiment of the effects of waves on glacial lake outburst induced by waves overtopping[J]. Mountain Research, 2014, 32(2): 241-248. [黄金辉,刘建康,程尊兰,等. 涌浪规模对冰碛湖溃决的影响实验[J]. 山地学报, 2014, 32(2): 241-248.]
[21]
Yan Jiahong.Rish analysis of moraine lake outburst under seismic source waves[D].Chengdu:Soutwest Jiaotong University,2016.
严佳宏.地震涌浪作用下冰碛堰塞湖溃决风险分析[D].成都:西南交通大学,2016.
[22]
Amicarelli A,Kocak B,Sibilla S,et al. A 3D smoothed particle hydrodynamics model for erosional dam-break floods[J]. International Journal of Computational Fluid Dynamics, 2017, 31(10): 413-434. DOI:10.1080/10618562.2017.1422731
[23]
Liu Wei,He Siming. A two-layer model for simulating landslide dam over mobile river beds[J]. Landslides, 2016, 13(3): 565-576.
[24]
Liang Qiuhua,Borthwick A G. Adaptive quadtree simulation of shallow flows with wet-dry fronts over complex topography[J]. Computers & Fluids, 2009, 38(2): 221-234. DOI:10.1016/j.compfluid.2008.02.008
[25]
Liu Wei,He Siming,Li Xinpo,et al. Two-dimensional landslide dynamic simulation based on a velocity-weakening friction law[J]. Landslides, 2016, 13(5): 957-965. DOI:10.1007/s10346-015-0632-z
[26]
Zou Q,Su Z M,Zhu X H.Mechanism of landslide-debris flow-barrier lake disaster chain after the Wenchuan earthquake[C]//Earthquake-Induced Landslides.Berlin:Springer,2013:917–924.
[27]
Cao Zhixian,Yue Zhiyuan,Pender G.Landslide dam failure and flood hydraulics. Part Ⅱ:Coupled mathematical modeling[J]. Natural Hazards, 2011, 59(2): 1021-1045. DOI:10.1007/s11069-011-9815-7