工程科学与技术   2020, Vol. 52 Issue (2): 29-37
金沙江白格堰塞湖溃决过程数值模拟
钟启明1,2, 陈生水1,2, 单熠博1     
1. 南京水利科学研究院,江苏 南京 210029;
2. 水利部土石坝破坏机理与防控技术重点实验室,江苏 南京 210029
基金项目: 国家重点研发计划项目(2018YFC1508604);国家自然科学基金项目(51779153;51539006);中央级公益性科研院所基本科研业务费专项资金(Y319003)
摘要: 2018年10月10日和11月3日,在中国四川省与西藏自治区交界处的白格村同一位置连续发生两次滑坡,完全堵塞了金沙江形成堰塞湖。由于“10·10”滑坡形成的堰塞湖水位抬升迅速,堰塞湖于10月12日自然漫顶溃决。“11·3”滑坡堵塞了“10·10”堰塞体溃决形成的流道,形成了更大的堰塞湖,鉴于客观条件允许,采取了开挖泄流槽降低堰塞湖溃决水位的措施,至11月12日,堰塞湖发生漫顶溃决,溃口洪水峰值流量为31 000 m3/s。由于“11·3”白格堰塞湖溃决案例拥有较为完整的实测资料,为堰塞湖溃决过程的研究提供了宝贵的基础数据。基于堰塞体的溃决机理,建立了可考虑堰塞湖的水动力条件、堰塞体的形态和材料特征的堰塞湖溃决过程数学模型。该模型采用宽顶堰公式模拟溃口洪水流量,并根据堰塞湖入湖和溃口流量以及堰塞湖的水位–湖面面积关系曲线确定堰塞湖水位的变化;采用基于水流剪应力和堰塞体材料临界剪应力、可考虑宽级配堰塞体材料特性的冲蚀公式模拟材料的冲蚀过程;假设溃口在纵向下切和横向展宽过程中坡角保持不变,采用极限平衡法分析溃口在发展过程中发生的边坡失稳;采用按时间步长迭代的数值计算方法模拟堰塞湖溃决时的水土耦合过程。采用作者建立的模型对“11·3”白格堰塞湖溃决案例进行反演分析后发现,本模型计算获得的溃口流量过程、堰塞湖水位变化过程、溃口发展过程与实测值基本吻合。参数敏感性分析表明:冲蚀系数对溃决过程具有重要的影响;残留坝高通过影响下泄库容也对溃决过程产生作用;开挖泄流槽可降低堰塞湖溃决时的库容,从而对溃口流量过程产生影响,是降低灾害损失的有效手段。
关键词: 白格    堰塞湖    溃决    数值模拟    参数敏感性分析    
Numerical Modeling of Breaching Process of Baige Dammed Lake on Jinsha River
ZHONG Qiming1,2, CHEN Shengshui1,2, SHAN Yibo1     
1. Nanjing Hydraulic Research Inst., Nanjing 210029, China;
2. Key Lab. of Earth–Rock Dam Failure Mechanism and Safety Control Techniques, Ministry of Water Resources, Nanjing 210029, China
Abstract: On October 10th and November 3rd, 2018, two successive landslides occurred at the same place of Baige village, the border of Sichuan Province and Tibet Autonomous Region, China, which blocked Jinsha River twice. Owing to the rapid rising of water level of the “10·10” dammed lake, on October 12th, the landslide dam breached naturally. Then, the flow channel formed after the breaching of “10·10” dammed lake was blocked by the “11·3” landslide, resulting in an even larger dammed lake. The measures of excavation of drainage channel were taken to decrease the water level of the dammed lake. On November 12th, the dammed lake breached with the peak flow rate of 31 000 m3/s. For the “11·3” Baige dammed lake, the detailed hydrological data were documented, which provided valuable basic data for the study of the breaching process of dammed lake. Based on the breaching mechanism of landslide-formed dam, a numerical model which can consider the hydrodynamic conditions of the lake and the morphological characteristics of dam materials of the landslide dam, was developed. The broad-crested weir equation was adopted to simulate the breach flow discharge, and the variation of water level of the dammed lake was determined according to inflow and breaching flow, as well as the rating curve of water level and surface area of the dammed lake. Based on the shear stress of water flow and the critical shear stress of dam materials, the erosion process was simulated by using the erosion formula for wide-graded dam materials. Under the assumption that the breach slope angles remain unchanged during the longitudinal cutting and transverse broadening processes, the limit equilibrium method is used to analyze the slope instability during the breach evolution process. The coupling process of soil and water during dam breaching was simulated by the time step iterative numerical method. Then, the back analysis of the “11·3” Baige dammed lake was conducted, and the comparison showed that the breach hydrograph, water level variation of dammed lake, and breach morphology evolution calculated by the proposed model are consistent with the measured data. The sensitivity analysis showed that soil erodibility had significant influence on the breaching process, and residual dam height also played a role in the breaching process by affecting the discharged water storage. In addition, the excavation of spillway can reduce the reservoir capacity when the dammed lake breached, thus affecting the breach hydrograph, which is an effective measure to reduce disaster losses.
Key words: Baige    dammed lake    dam breach    numerical modeling    sensitivity analysis of parameters    

堰塞湖一般是由于降雨或地震导致山体滑坡堵塞河道形成,特殊的形成条件决定了堰塞湖在世界各国山区均有分布,而中国独特的自然地理环境更容易形成堰塞湖[1]。其中,仅2008年中国汶川地震时就形成了257个堰塞湖[2];2018年10至11月,中国金沙江和雅鲁藏布江连续发生滑坡堵江事件,引起了国内外广泛关注。众所周知,堰塞湖溃决过程模拟对致灾后果评价和应急预案的制定具有重要的意义,目前,堰塞湖溃决过程数学模型主要可分为参数模型和基于物理机理的数学模型[3-4]。与人工填筑的土石坝相比,堰塞体拥有独特的结构形态和材料特性,一般而言,堰塞体具有更大的长宽比(顺河向/横河向)[5-7],堰塞体材料分布不均且具有宽级配特征[8-10];由于堰塞体未经人工碾压,故材料的冲蚀特性沿深度方向存在着变化[11]。因此,直接采用人工水库大坝的溃坝模型模拟堰塞湖的溃决过程可能会出现较大的误差。

参数模型一般基于历史溃坝案例,通过逻辑回归方程模拟溃决特征参数(主要包括溃口峰值流量、溃口最终尺寸与溃决历时),目前的参数模型多用于描述土石坝的溃决[12],专门针对堰塞体溃决的参数模型较少。Costa[5,13]最早提出了预测堰塞体溃口峰值流量的参数模型;Peng等[8]基于52座有实测资料的已溃堰塞体,也构建了专门用于预测堰塞体溃口峰值流量、溃口最终尺寸与溃决历时的参数模型。值得一提的是,近年来,国内学者基于易贡、唐家山、小岗剑等堰塞湖溃决案例,开发了一些基于物理机理的数学模型[9-11],为堰塞体溃决过程的模拟提供了新的有效手段。

通过堰塞体溃决案例的实测资料和物理模型试验发现[11,14-16],堰塞体的冲蚀模式主要表现为溯源冲蚀和沿程侵蚀,且溃决过程中堰塞体下游坝坡坡角呈逐渐减小的趋势;由于堰塞体材料的宽级配特性,随着溃决后期堰塞湖水动力条件的减弱,水流无法进一步冲蚀堰塞体,溃决过程随之结束。另外,堰塞体地勘资料表明[8,17],堰塞体沿高程自上而下存在着材料逐渐粗化的现象,因此在模拟时应根据堰塞体材料的分布特征将堰塞体分层处理,这样才能合理模拟堰塞体材料的冲蚀特性。

本文基于堰塞体的溃决机理,考虑堰塞体形态结构特征、材料特性和堰塞湖的水动力条件,建立了一个可描述堰塞体溃决过程的数学模型,对拥有完整实测资料的“11·3”白格堰塞湖的溃决过程进行反演分析,验证数值模型的合理性,并对影响堰塞湖溃决洪水流量过程的主要因素进行分析。

1 白格堰塞湖溃决过程概述

2018年10月10日22:06,中国四川省与西藏自治区交界的白格村(北纬31°4′56.41″,东经98°42′17.98″)金沙江右岸发生滑坡[1,18-20],形成堰塞体,完全堵塞了金沙江,据实测分析,“10·10”白格堰塞体的体积约为2 750×104 m3[21]。由于水位上涨迅速,且道路不畅,无法及时采取人工干预措施,堰塞湖于10月12日自然漫顶溃决[21-22]。2018年11月3日17:21,同一地点发生了二次滑坡[23-24]。“11·3”滑坡体的体积约为160×104 m3,滑坡体在下滑过程中沿途裹挟约140×104 m3的坡面残留碎屑和风化岩石高速下滑,堵塞“10·10”堰塞湖溃决形成的新流道,造成二次堵江[21]图1)。两次白格堰塞湖的形态特征参数见表1图2给出了两次堰塞体的断面形态示意图[25]

图1 两次滑坡形成的白格堰塞体 Fig. 1 Baige landslide dams formed by two subsequential landslides

表1 白格堰塞湖形态特征 Tab. 1 Morphological characteristics of Baige dammed lakes

图2 两次堰塞体断面形态示意图 Fig. 2 Schematic diagrams of cross sections of the two landslide dams

针对第2次白格堰塞体高度及堰塞湖水位上涨速度,采取了开挖泄流槽降低堰塞湖溃决时水位的人工干预措施。2018年11月8日,泄流槽开始施工,至11月11日11:00泄流槽开挖完成,泄流槽顺河向长220 m,底宽3 m,最大深度15 m,泄流槽边坡坡比1∶1.3[21]。11月12日04:45,堰塞湖发生溃决。

“11·3”白格堰塞湖形成后,利用干流已建水文站(岗拖、波罗、叶巴滩、巴塘、奔子栏、塔城、石鼓、上虎跳峡),同时,在堰塞体上、下游3 km处分别设立了水位观测站[26],构成了应急监测站网,用于测报堰塞湖入湖流量、溃决洪水流量和水位变化。对于入湖流量,采用预报模型参数移植和上下游水文站倍比放缩与耦合气象数值预报的分布式水文模型相结合的方法进行预报[26];对于溃口流量,采用地理信息系统(geographic information system,GIS)和数字高程模型(digital elevation model,DEM)网格数据相结合的快速空间信息处理技术,获取了水位–库容关系曲线[26],并通过实测堰塞湖的水位变化推求。

“11·3”白格堰塞湖形成时的水位为2 892.84 m,此时金沙江的流量约为700 m3/s[21]。利用构建的应急监测网站,全程监测了“11·3”白格堰塞湖的溃决过程,为相关研究提供了宝贵的实测基础数据。“11·3”白格堰塞湖的溃决过程可简要描述如下:2018年11月12日04:45,堰塞湖水位上升至2 952.52 m,达到泄流槽进口处底高程,湖水进入泄流槽;10:50,泄流槽全程进水,流量约为1~3 m3/s,泄流槽入口处水流流速约为1.0~1.5 m/s;入夜,水流流速增加至1.5~2.4 m/s;11月13日10:00流速增加至3.0 m/s;13:45,堰塞湖达到最高水位2 956.4 m,相应库容为5.78×108 m3;14:30,泄流槽入口处水流流速上升至6.5 m/s;18:00,溃口出现峰值流量,为31 000 m3/s,此时的水流流速为10 m/s;20:00,溃口流量下降至7 700 m3/s;截止于2018年11月14日08:00,溃口流量与基流基本一致,堰塞湖水位下降至2 905.75 m,溃决过程停止。“11·3”白格堰塞湖溃决过程如图3所示。

图3 “11·3”白格堰塞湖溃决过程实拍图 Fig. 3 Photos of breaching process of “11·3” Baige dammed lake

2 堰塞湖溃决过程数学模型

堰塞湖溃决过程数学模型主要包括3个模块:水动力模块、材料冲蚀模块和溃口发展模块。模型采用按时间步长迭代的数值计算方法模拟堰塞体溃决过程中的水土耦合作用。

2.1 水动力模块

堰塞体在溃决过程中堰塞湖需满足水量平衡关系,可用式(1)表示:

${A_{\rm{s}}}\left( {{{\textit{z}}_{\rm{s}}}} \right)\frac{{{\rm{d}}{{\textit{z}}_{\rm{s}}}}}{{{\rm{d}}t}} = {Q_{{\rm{in}}}} - {Q_{\rm{b}}}$ (1)

式中,As(zs)为不同水位对应的堰塞湖湖面面积,zs为堰塞湖水位,t为时间,Qin为上游来流量,Qb为溃口流量。

模型采用宽顶堰流量公式模拟溃口流量过程[27]

${Q_{\rm{b}}} = 1.7b{H^{1.5}} + 1.3m{H^{2.5}}$ (2)

式中:b为溃口底部宽度;H为溃口水深,可表示为zszb,其中,zb为溃口底部高程;m为溃口边坡坡比的倒数(水平向/垂直向)。

2.2 材料冲蚀模块

模型采用基于剪应力的冲蚀公式模拟每层堰塞体材料的冲蚀过程[28]

$\frac{{{\rm{d}}{{\textit{z}}_{\rm{b}}}}}{{{\rm{d}}t}} = {k_{\rm{d}}}\left( {{\tau _{\rm{b}}} - {\tau _{\rm{c}}}} \right)$ (3)

式中,kd为冲蚀系数,τb为溃口底床处的水流剪应力,τc为堰塞体材料的临界剪应力。

冲蚀系数与堰塞体材料的临界剪应力可通过试验[17,29-30]或采用经验公式获取。针对堰塞体材料的宽级配特性,模型分别采用式(4)、(5)计算kd[17]τc[31]

${k_{\rm{d}}} = 20\;075{e^{4.77}}C_{\rm{u}}^{ - 0.76}$ (4)

式中:e为堰塞体材料的孔隙比;Cu为堰塞体材料的不均匀系数,可表示为Cu=d60/d10,其中,d60d10分别表示小于该粒径的试样质量占试样总质量的60%和10%对应的材料粒径。

${\tau _{\rm{c}}} = \frac{2}{3}g{d_{50}}\left( {{\rho _{\rm{s}}} - {\rho _{\rm{w}}}} \right)\tan\; \varphi $ (5)

式中,g为重力加速度,d50为堰塞体材料的平均粒径, $\;\rho_{\rm{s}} $ 堰塞体材料的密度, $\;\rho_{\rm{w}} $ 为水的密度, $\varphi $ 为堰塞体材料的内摩擦角。

溃口底床处的水流剪应力可采用曼宁公式计算[32]

${\tau _{\rm{b}}} = \frac{{{\rho _{\rm{w}}}g{n^2}Q_{\rm{b}}^2}}{{{A^2}{R^{{1 / 3}}}}}$ (6)

式中:n为曼宁糙率系数,可表示为 $n = d_{50}^{1/6}/12 $ [32] $A $ 为过流断面面积;R为过流断面水力半径。

2.3 溃口发展模块

通过白格堰塞体溃决过程录像发现,溃口坡角在冲蚀过程中基本上保持一致直至发生边坡失稳。为了便于计算分析,模型假设在溃决水流作用下,对于每层材料,溃口的纵向下切与横向扩展速度一致,并且溃口边坡坡角不变直至失稳,溃口顶宽和底宽横向扩展速度可表示为:

$\frac{{{\rm{d}}B}}{{{\rm{d}}t}} = \frac{{{n_{{\rm{loc}}}} \cdot \left( {{{{\rm{d}}{{\textit{z}}_{\rm{b}}}} / {{\rm{d}}t}}} \right)}}{{\sin\; \beta }}$ (7)

式中:B为溃口顶宽;nloc为溃口位置参数,其中,nloc=1表示单侧冲蚀,nloc=2表示两侧冲蚀; $\;\beta $ 为溃口边坡坡角。

$\frac{{{\rm{d}}b}}{{{\rm{d}}t}} = {n_{{\rm{loc}}}}\frac{{{\rm{d}}{{\textit{z}}_{\rm{b}}}}}{{{\rm{d}}t}}\left( {\frac{1}{{\sin\; \beta }} - \frac{1}{{\tan\; \beta }}} \right)$ (8)

针对堰塞体在溃决过程中顺河向坡角逐渐变化的情况,模型假设堰塞体下游坡坡脚处的冲蚀深度为0,并沿下游坡向上线性增加直至堰塞体顶部,由于顶部的冲蚀深度可通过式(3)计算,因此下游坡坡角的变化过程可通过计算获取。

溃口的持续下切与横向扩展会导致边坡的失稳,基于白格堰塞体溃决过程的现场视频,采用极限平衡法模拟溃口边坡的失稳。另外,由于模型对于堰塞体材料分层考虑,每层材料的物理力学指标不尽相同,为了便于分析,模型假设滑动面为平面,边坡失稳后的坡角以溃口所在层的材料确定。当满足式(9)时,溃口边坡发生失稳:

${F_{\rm{d}}} > {F_{\rm{r}}}$ (9)

式中,Fd为边坡失稳的驱动力,Fr为边坡失稳的抵抗力。可分别表示如下:

${F_{\rm{d}}} = W\sin\; \alpha = \frac{1}{2}{\gamma _{\rm{s}}}H_{\rm{s}}^2\left( {\frac{1}{{\tan \alpha }} - \frac{1}{{\tan\; \beta }}} \right)\sin\; \alpha $ (10)
$\begin{aligned}[b]{F_{\rm{r}}} =& W\cos\; \alpha \tan\; \varphi + \frac{{C{H_{\rm{s}}}}}{{\sin\; \alpha }}= \\ & \frac{1}{2}{\gamma _{\rm{s}}}H_{\rm{s}}^2\left( {\frac{1}{{\tan\; \alpha }} - \frac{1}{{\tan\; \beta }}} \right)\cos\; \alpha \tan\; \varphi + \frac{{C{H_{\rm{s}}}}}{{\sin\; \alpha }} \end{aligned}\!\!\!$ (11)

式中, $W $ 为失稳块体的重量, $\alpha $ 为失稳后新溃口边坡的坡角, $\gamma _{\rm{s}} $ 为堰塞体材料的容重,Hs为溃口边坡高度,C为堰塞体材料的黏聚力。

由于堰塞体具有宽级配的特性,因此在溃决后可能存在残留坝高,模型有2种计算方法可供选用:1)通过计算水流在溃口底床的切应力和每层堰塞体材料的临界剪应力确定溃口能否继续下切;2)通过设定溃口最终底高程对残留坝高进行设置,当溃口下切至设置高程后无法继续向下发展,只能横向扩展。

对于“11·3”白格堰塞体,地勘资料表明[25],堰塞体颗粒总体较细,但下部颗粒粒径较大。由于堰塞体下部粗颗粒抗冲蚀能力较强,因此溃决过程中主要是上部细颗粒发生冲蚀;另外,对于本次反演分析,由于已知溃口最终底高程,采用第2)种方法确定残留坝高。

2.4 数值计算方法

模型采用按时间步长迭代的数值计算方法模拟堰塞体溃决过程中的水土耦合作用,输入初始参数,设置计算时长tc和时间步长Δt,采用如下步骤计算堰塞体的溃决过程:1)采用式(1)和(2)计算Qb;2)采用式(1)计算dzs/dt;3)采用式(3)~(6)计算dzb/dt;4)分别采用式(7)和(8)计算dB/dt和db/dt;5)在每个时间步结束时分别采用zs–(dzs/dt)Δtzb–(dzb/dt)ΔtB–(dB/dt)Δtb–(db/dt)Δt计算下一个步长输入的zszbBb;6)在每个时间步结束时采用式(9)~(11)校核溃口边坡的稳定性;7)开始下一个循环的迭代计算,直至计算时长结束。模型的计算流程见图4

图4 模型计算流程图 Fig. 4 Flow chart of the model calculation

3 “11·3”白格堰塞湖溃决过程反演分析 3.1 输入参数

选择“11·3”白格堰塞体垭口高程为堰塞体顶高程,以原河床底高程2 870 m为堰塞体底高程,通过地勘资料,确定堰塞体的高度为96 m,堰塞体顶部对应的高程为2 966 m,堰塞体顺河向顶宽和横河向长度分别为200 m和700 m。由于堰塞体的形态不规则,为了便于计算分析,选取初始溃口处的断面形态作为代表,上、下游坡坡比分别设定为1∶2.7和1∶5.5。根据堰塞湖的水位–库容关系曲线[22],推导出堰塞湖水位–湖面面积关系曲线如图5所示。为了模拟堰塞湖溃决过程,初始湖水位设置为泄流槽入口处底高程2 952.52 m,初始溃口(泄流槽)底宽为3 m;通过比较堰塞体顶高程和初始水位,可得初始溃口深度为13.48 m,溃口边坡坡比为1∶1.3(垂直向/水平向)。综合考虑泄流槽的底高程和堰塞体材料的特性,选择文献[25]提供的L4–2(高程为2 960 m)位置处的材料物理力学特性作为模型的输入参数。并通过式(4)和(5),得出冲蚀系数和堰塞体材料的临界剪应力分别为368 mm3/(N·s)和13.6 Pa。实测资料显示,堰塞湖溃决后的最终湖水位为2 905.75 m,假设堰塞体的残留坝高为35.75 m。模型的计算参数见表2

表2 数学模型输入参数 Tab. 2 Input parameters of numerical model

图5 白格堰塞湖水位–湖面面积关系曲线 Fig. 5 Water level and water surface area curve of Baige dammed lake

3.2 计算结果与实测结果比较分析

对于堰塞湖的溃决模拟,溃口峰值流量(Qp)、溃口最终顶宽(Bf)、溃口最终底宽(bf)、溃口最终深度(Df)以及溃口峰值流量出现时间(Tp)是重要的输出参数,并且,溃口洪水流量过程线是合理分析堰塞湖溃决洪水下游演进情况和致灾后果评价的重要基础,因此上述参量的准确性是模型合理性的重要依据。下面主要通过对比计算结果和实测结果,验证建立的数学模型的合理性。表3分别给出QpBfbfTp等4个参数的比较情况,由于模型设置了残留坝高,因此不对参数Df进行比较。另外,图67对溃口流量过程和堰塞湖水位变化过程的计算值和实测值进行了比较。由于堰塞湖溃决过程中水下的溃口未能测量,故只提供计算获取的溃口发展过程,并与堰塞湖溃决后的实测溃口最终尺寸进行了比较(图8)。值得一提的是,通过联解水量平衡方程与蓄泄关系,搭建静库容调洪演算模型,根据10、15、30 min等不同时间步长的水位变幅,推算多时间尺度的流量过程,长江水利委员会水文局综合获取了最终溃口处的流量过程,并通过搭建动库容调洪演算模型复核了计算结果[26]

表3 白格堰塞湖溃决输出参数实测值与计算值比较 Tab. 3 Comparison of calculated and measured data of output parameters for Baige dammed lake

图6 溃口流量过程实测值与计算值比较 Fig. 6 Comparison of measured and calculated breach hydrographs

图7 堰塞湖水位变化过程实测值与计算值比较 Fig. 7 Comparison of measured and calculated variations of water level in dammed lake

图8 溃口发展过程实测值与计算值比较 Fig. 8 Comparison of measured and calculated breach evolution processes

表3图68可以看出:所建立的数学模型可较好地反演“11·3”白格堰塞湖溃决的全过程,各重要输出参数的最大误差控制在±10%以内;另外,溃口流量过程及堰塞湖水位变化过程的计算值与实测值基本吻合,验证了模型的合理性。

4 参数敏感性分析

为了更进一步研究模型的核心输入参数对堰塞湖溃决过程的影响,开展了参数敏感性分析。大量研究表明,冲蚀系数对溃坝过程有着重要的影响,为了检验冲蚀系数对“11·3”白格堰塞湖溃决过程的影响,对冲蚀系数开展参数敏感性分析。堰塞湖溃决后不同的残留坝高对应不同的下泄库容,进而对堰塞体溃决过程产生影响,也对不同残留坝高情况下的堰塞湖溃决过程进行模拟,分析残留坝高对溃决过程的影响。另外,为了验证泄流槽开挖对堰塞湖减灾的效果,模拟了不开挖泄流槽情况下堰塞体的溃决过程,分析泄流槽对溃决过程的影响。其他的参数,如平均粒径、堰塞体材料的黏聚力和内摩擦角等对堰塞湖的溃决过程影响较小[33],在此不进行分析。

4.1 冲蚀系数敏感性分析

为了研究冲蚀系数的敏感性,分别将kd缩小和扩大1倍,其他参数保持不变,然后采用建立的数学模型对“11·3”白格堰塞湖溃决过程进行模拟。冲蚀系数敏感性分析结果见表4,其中相对误差为kd变化后计算结果与原始结果(kd=1.0时)的比较;图9同时给出了不同冲蚀系数情况下的溃口洪水流量过程。通过表4图9可以看出,堰塞湖溃决各输出参数对冲蚀系数均较为敏感,其中溃口峰值流量最为敏感,说明了冲蚀系数的合理选择对堰塞湖溃决过程模拟的重要性。

表4 不同冲蚀系数对白格堰塞湖溃决过程的影响 Tab. 4 Influence of different soil erosion coefficients for the breach process of Baige dammed lake

图9 不同冲蚀系数条件下的溃口流量过程 Fig. 9 Breach hydrographs for different soil erodibility coefficients

4.2 残留坝高敏感性分析

为了研究堰塞体溃决后残留坝高对溃决过程的影响,将堰塞体溃决后最终底部高程由2 905.75m变为2 934.00 m(1/3溃坝)或2 870.00 m(全溃坝),其他参数保持不变。表5给出了残留坝高敏感性的分析结果,图10同时给出了不同残留坝高情况下的溃口洪水流量过程。由表5可知,随着残留坝高的降低,溃口峰值流量和溃口最终顶宽呈上升趋势,而溃口底宽逐渐减小,峰值流量出现时间逐渐滞后。究其原因,应与堰塞湖溃决时的水动力条件相关。由图10可以看出:随着高程的降低,白格堰塞湖库面面积迅速减小,因此当溃口最终底高程由2 905.75 m降低至2 870.00 m时,溃口峰值流量仅增加12.4%,相较于溃口最终底高程由2 934.00 m降低至2 905.75 m时的溃口峰值流量明显减小。究其原因,应为随着堰塞湖水位的快速下降,水动力条件明显减弱,溃口顶宽增大趋势减缓,而由于溃口下切深度逐渐增大,底宽反而呈逐渐减小的趋势;由于溃口下切过程中水动力条件逐渐减弱,峰值流量出现的时间也逐渐滞后。

表5 不同残留坝高对堰塞湖溃决过程的影响 Tab. 5 Influence of different residual dam heights for the breach process of Baige dammed lake

图10 不同残留坝高条件下的溃口流量过程 Fig. 10 Breach hydrographs for different residual dam heights

4.3 泄流槽开挖敏感性分析

对于大多数堰塞湖,受制于事件突发和交通闭塞,常无法采取有效的人工干预措施。实践证明,在客观条件允许的情况下,通过开挖泄流槽的方式可有效降低堰塞湖溃决时的水位,进而起到减小下泄库容和降低溃口峰值流量的作用,以减轻堰塞湖溃决的致灾后果。为了研究泄流槽的开挖效果,采用建立的数学模型模拟未开挖泄流槽的溃决过程,假设上游来流一直保持700 m3/s,模型的其他输入参数保持不变。表6给出了开挖泄流槽与否的对比分析结果,图11同时给出了泄流槽是否开挖情况下的溃口洪水流量过程。根据堰塞湖入流情况,不开挖泄流槽时,白格堰塞湖约在11月15日15:30发生自由漫溢,在11月17日22:37 溃口出现峰值流量,下泄水量约为8.22×108 m3;而开挖泄流槽时的下泄水量为5.84×108 m3表6图11的泄流槽开挖敏感性分析结果表明,对于库容较大的堰塞湖,开挖泄流槽可明显降低溃口峰值流量,减小溃口尺寸,在条件允许的情况下,可作为人工干涉的首选方式。

表6 开挖泄流槽与否对堰塞湖溃决过程的影响 Tab. 6 Influence of excavating drainage channel or not for the breach process of Baige dammed lake

图11 开挖泄流槽与否情况下的溃口流量过程 Fig. 11 Breach hydrographs for excavating drainage channel or not

5 结 论

本文基于堰塞体的溃决机理,建立了包含水动力模块、材料冲蚀模块和溃口发展模块的堰塞体溃决过程数学模型,并采用按时间步长迭代的数值计算方法对堰塞体溃决过程中的水土耦合过程进行模拟。采用建立的模型对拥有完整水文实测资料的“11·3”白格堰塞湖的溃决过程进行了反演分析,通过与实测溃口流量过程及溃口最终尺寸的对比,验证了数学模型的合理性。

参数敏感性分析表明:冲蚀系数对于堰塞体溃决过程的模拟具有至关重要的作用,对该参数的选择应给予充分重视,而目前大多采用经验公式或小尺度试验来获取该参数,故宽级配堰塞体材料的冲蚀特性应进一步深入研究。残留坝高对堰塞湖溃口后的下泄库容有着重要影响,本文模型提供了两种残留坝高的确定方法,可为堰塞体溃决过程的模拟提供参考。另外,对于库容较大的堰塞湖,在客观条件允许的情况下,泄流槽的开挖被证明是一种行之有效的减轻堰塞湖溃决洪水灾害的人工干预措施,可作为堰塞湖减灾措施的首选方案。

参考文献
[1]
刘宁,杨启贵,陈祖煜.堰塞湖风险处置[M].武汉:长江出版社,2016.
[2]
Cui P,Zhu Y Y,Han Y S,et al. The 12 May Wenchuan earthquake-induced landslide-dammed lakes:Distribution and preliminary risk evaluation[J]. Landslides, 2009, 7(6): 209-223. DOI:10.1007/s10346-009-0160-9
[3]
Wu W M. Earthen embankment breaching[J]. Journal of Hydraulic Engineering, 2011, 137(12): 1549-1564. DOI:10.1061/(ASCE)HY.1943-7900.0000498
[4]
Zhong Q M,Wu W M,Chen S S,et al. Comparison of simplified physically based dam breach models[J]. Natural Hazards, 2016, 84(2): 1385-1418. DOI:10.1007/s11069-016-2492-9
[5]
Costa J E.Floods from dam failures[R].Washington:U.S.Geological Survey,1985.
[6]
Korup O. Recent research on landslide dams—A literature review with special attention to New Zealand[J]. Progress in Physical Geography, 2002, 26(2): 206-235. DOI:10.1191/0309133302pp333ra
[7]
Ermini L,Casagli N. Prediction of the behaviour of landslide dams using a geomorphological dimensionless index[J]. Earth Surface Processes and Landforms, 2003, 28(1): 31-47. DOI:10.1002/esp.424
[8]
Peng M,Zhang L M. Breaching parameters of landslide dams[J]. Landslides, 2012, 9(1): 13-31. DOI:10.1007/s10346-011-0271-y
[9]
Chen Z Y,Ma L Q,Yu S,et al. Back analysis of the draining process of the Tangjiashan barrier lake[J]. Journal of Hydraulic Engineering, 2015, 141(4): 05014011. DOI:10.1061/(ASCE)HY.1943-7900.0000965
[10]
Zhong Q M,Chen S S,Mei S A,et al. Numerical simulation of landslide dam breaching due to overtopping[J]. Landslides, 2018, 16(6): 1183-1192. DOI:10.1007/s10346-017-0935-3
[11]
Chang D S,Zhang L M. Simulation of the erosion process of landslide dams due to overtopping considering variations in soil erodibility along depth[J]. Natural Hazards and Earth System Sciences, 2010, 10(4): 933-946. DOI:10.5194/nhess-10-933-2010
[12]
Mei Shiang,Chen Shengshui,Zhong Qiming,et al. Parametric model for breaching analysis of earth–rock dam[J]. Advances Engineering Sciences, 2018, 50(2): 60-66. [梅世昂,陈生水,钟启明,等. 土石坝溃坝参数模型研究[J]. 工程科学与技术, 2018, 50(2): 60-66. DOI:10.15961/j.jsuese.201700828]
[13]
Costa J E,Schuster R L.The formation and failure of natural dams[J].Geological Society of America Bulletin,1988,100(7):1054–1068.
[14]
Liu N,Chen Z Y,Zhang J X,et al. Draining the Tangjiashan barrier lake[J]. Journal of Hydraulic Engineering, 2010, 136(11): 914-923. DOI:10.1061/(ASCE)HY.1943-7900.0000241
[15]
Niu Z P,Xu W L,Li N W,et al. Experimental investigation of the failure of cascade landslide dams[J]. Journal of Hydrodynamics, 2012, 24(3): 430-441. DOI:10.1016/S1001-6058(11)60264-3
[16]
Zhou G G D,Zhou M J,Shrestha M S,et al. Experimental investigation on the longitudinal evolution of landslide dam breaching and outburst floods[J]. Geomorphology, 2019, 334: 29-43. DOI:10.1016/j.geomorph.2019.02.035
[17]
Chang D S,Zhang L M,Xu Y,et al. Field testing of erodibility of two landslide dams trigged by the 12 May Wenchuan earthquake[J]. Landslides, 2011, 8(3): 321-332. DOI:10.1007/s10346-011-0256-x
[18]
Deng Jianhui,Gao Yunjian,Yu Zhiqiu,et al. Analysis on the formation mechanism and process of Baige landslides damming the upper reach of Jinsha River,China[J]. Advances Engineering Sciences, 2019, 51(1): 9-16. [邓建辉,高云建,余志球,等. 堰塞金沙江上游的白格滑坡形成机制与过程分析[J]. 工程科学与技术, 2019, 51(1): 9-16. DOI:10.15961/j.jsuese.201801438]
[19]
Ouyang C J,An H C,Zhou S,et al. Insights from the failure and dynamic characteristics of two sequential landslides at Baige village along the Jinsha River,China[J]. Landslides, 2019, 16(7): 1397-1414. DOI:10.1007/s10346-019-01177-9
[20]
Zhang Z,He S M,Liu W,et al. Source characteristics and dynamics of the October 2018 Baige landslide revealed by broadband seismograms[J]. Landslides, 2019, 16(4): 777-785. DOI:10.1007/s10346-019-01145-3
[21]
Cai Yaojun,Luan Yuesheng,Yang Qigui,et al. Study on structural morphology and dam-break characteristics of Baige barrier dam on Jinsha River[J]. Yangtze River, 2019, 50(3): 15-22. [蔡耀军,栾约生,杨启贵,等. 金沙江白格堰塞体结构形态与溃决特征研究[J]. 人民长江, 2019, 50(3): 15-22. DOI:10.16232/j.cnki.1001-4179.2019.03.004]
[22]
Chen Zuyu,Zhang Qiang,Hou Jingming,et al. Back analysis on dam-breach flood of “10·10” Baige barrier lake on Jinsha River[J]. Yangtze River, 2019, 50(5): 1-4. [陈祖煜,张强,侯精明,等. 金沙江“10·10”白格堰塞湖溃坝洪水反演分析[J]. 人民长江, 2019, 50(5): 1-4. DOI:10.16232/j.cnki.1001-4179.2019.05.001]
[23]
Fan X M,Xu Q,Alonso–Rodriguez A,et al. Successive landsliding and damming of the Jinsha River in eastern Tibet,China:Prime investigation,early warning,and emergency response[J]. Landslides, 2019, 16(5): 1003-1020. DOI:10.1007/s10346-019-01159-x
[24]
Liang G L,Wang Z,Zhang G W,et al. Two huge landslides that took place in quick succession within a month at the same location of Jinsha River[J]. Landslides, 2019, 16(5): 1059-1062. DOI:10.1007/s10346-019-01165-z
[25]
Zhang L M,Xiao T,He J,et al. Erosion-based analysis of breaching of Baige landslide dams on the Jinsha River,China,in 2018[J]. Landslides, 2019, 16(10): 1965-1979. DOI:10.1007/s10346-019-01247-y
[26]
Cheng Haiyun. Hydrology emergency monitoring and forecast on “11·3” Baige barrier lake,Jinsha River[J]. Yangtze River, 2019, 50(3): 23-27. [程海云. “11·3”金沙江白格堰塞湖水文应急监测预报[J]. 人民长江, 2019, 50(3): 23-27. DOI:10.16232/j.cnki.1001-4179.2019.03.005]
[27]
Singh V P.Dam breach modeling technology[M].Dordrecht:Kluwer Academic Publishes,1996.
[28]
Graf W H.Hydraulics of sediment transport[M].Littleton:Water Resources Publications,1984.
[29]
Briaud J L,Ting F C K,Chen H C,et al. Erosion function apparatus for scour rate prediction[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2001, 127(2): 105-113. DOI:10.1061/(ASCE)1090-0241(2001)127:2(105)
[30]
Hanson G J,Cook K R. Apparatus,test procedures,and analytical methods to measure soil erodibility in situ[J]. Applied Engineering in Agriculture, 2004, 20(4): 455-462. DOI:10.13031/2013.16492
[31]
Annandale G W.Scour technology mechanics and engineering practice[M].New York:Mcgraw–Hill,2006.
[32]
Wu W M. Simplified physically based model of earthen embankment breaching[J]. Journal of Hydraulic Engineering, 2013, 139(8): 837-851. DOI:10.1061/(ASCE)HY.1943-7900.0000741
[33]
Wang L,Chen Z Y,Wang N X,et al. Modeling lateral enlargement in dam breaches using slope stability analysis based on circular slip mode[J]. Engineering Geology, 2016, 209: 70-81. DOI:10.1016/j.enggeo.2016.04.027