崩滑堰塞湖的形成-孕灾-致灾机理与模拟方法

钟 启 明，钱 亚 俊，单 熠 博

(1.南京水利科学研究院，江苏 南京 210029； 2.水利部土石坝破坏机理与防控技术重点实验室，江苏 南京 210029)

1 研究背景

堰塞湖是一定量固体物质阻塞山区河谷或河道导致上游壅水形成的具有一定库容的水体，而堵塞河谷或河道的固体物质被称为堰塞体[1-2]。堰塞湖在世界范围内广泛分布，依据成因，堰塞湖可分为熔岩、滑坡、崩塌、泥石流和冰碛堰塞湖等5类[3]，全球350个大型堰塞湖分布如图1所示[4-5]。世界范围内1 393个堰塞湖案例的统计数据表明[6]，形成堰塞湖的诱因依次是地震(50.5%)、降雨(39.3%)、融雪(2.4%)、人为原因(2.2%)、火山喷发(0.9%)，其他未知原因的占4.7%。由此可以看出，地震和降雨是形成堰塞湖的主导因素，两种成因的堰塞湖约占总数的90%。

图1 全球350个大型堰塞湖分布图[4-5]Fig.1 The distribution map of 350 large dammed lakes around the world[4-5]

美国地质调查局Costa和Schuster对全球73个堰塞湖的寿命统计发现[1]，85%的堰塞湖寿命小于1 a。我国学者Peng和Zhang[7]、石振明等[8]、Shen等[6]也分别通过对全球204，276个和352个堰塞湖的寿命统计得出了类似结论。堰塞湖的破坏模式包括水流漫顶冲刷、渗透破坏或边坡失稳，其中89%为水流漫顶冲刷溃决，10%为渗透破坏溃决[1]。

近年来，受地形地貌、地质构造及气象水文等条件综合作用，我国堰塞湖事件呈多发、频发态势。据作者统计，21世纪以来(2000年1月至2019年12月)，我国发生有文献记载的堰塞湖案例362个。2000年4月9日，西藏波密县易贡乡发生巨型山体滑坡，形成易贡堰塞湖[2]，并于同年6月10日溃决，21亿m3洪水下泄，溃口峰值流量达到12.4万m3/s，导致我国墨脱、波密、林芝3县90余乡近万人受灾，印度布拉马普特拉河沿岸7个邦94人死亡，250万人无家可归[5，9]；2008年，“5·12”汶川地震形成了257处崩滑堰塞湖[10]，其中唐家山是集雨面积最广、蓄水量最大、威胁最严重的堰塞湖，在人工干预下于2008年6月7日应急泄流，由于泄流前紧急转移下游风险人口约25万人，所幸未造成人员伤亡[11]；2018年10～11月，我国金沙江和雅鲁藏布江各接连发生两次滑坡事件，形成了白格和加拉堰塞湖，并在短期内发生溃决，对西藏、四川和云南的人民群众生命财产安全构成巨大威胁[12-13]。

与人工填筑的土石坝不同，堰塞体一般由崩滑土石料快速堆积而成，没有经过充分压实，结构较为松垮、组成物质杂乱，局部存在由大颗粒骨架架空形成的高渗透区域，渗流和力学稳定性较差[14]，且缺乏必要的洪水溢流设施，容易发生溃决造成严重的洪水灾害，对下游公众生命财产和基础设施安全构成巨大威胁。

因此，快速地评估堰塞体的稳定性、准确地预测溃口流量过程对堰塞湖的科学应急处置至关重要。本文针对崩塌型和滑坡型这两类发生频率最高的堰塞湖，重点围绕崩滑堰塞湖形成-孕灾-致灾机理与模拟方法开展研究，总结其形成机理，剖析影响其稳定性的关键特征参数，并介绍作者提出的堰塞体稳定性快速评价方法和堰塞湖溃决过程数值模拟方法，选择我国21世纪形成的3个典型崩滑堰塞湖案例，利用作者的研究成果对其灾害链进行模拟分析。

2 崩滑堰塞湖形成机理与堰塞体颗粒组成

2.1 崩滑堰塞湖形成机理

研究表明，崩滑堰塞体的几何形态、粒径分布和材料力学性质等重要特征取决于其形成过程，并极大地影响了堰塞湖的孕灾和致灾过程[15-16]。根据国内外常用的分类方法[17-19]，并考虑崩滑体的运动特征，可将崩滑堰塞湖的形成原因细分为3类：崩塌、滑坡和碎屑流。

崩塌形成堰塞湖是指岸坡上部岩土体被裂缝切割、拉裂后，在外荷载作用下失去稳定，脱离母岩急剧向下翻滚、跳跃，从而堵塞河道[19-20]。滑坡形成堰塞湖是指岸坡在外荷载作用下，当坡内软弱结构面上的剪应力超过了抗剪强度，滑坡体沿滑裂面发生整体滑动。一般可将滑坡形成堰塞湖的过程概括为岸坡上部拉裂、下部滑移隆起、中段快速剪断及整体的高速滑动[21]。碎屑流形成堰塞湖是指在滑坡或崩塌过程中，坡体物质在高速差异滑移速度和相互碰撞作用下产生平移剪切运动、跳跃及滚动等综合流动形式，将高速滑坡碎屑流冲到河对岸并阻挡停积于河床上，从而堵塞河道[18-19]。

2.2 堰塞体颗粒分布特征

目前，国内外对崩滑堰塞湖基本特性的研究主要集中在堰塞湖事件的统计与类型识别、堰塞体的形态特征、堰塞湖形成的影响因素等方面[16]。调查研究发现[2，5]，堰塞体的颗粒组成对其稳定性和溃决过程有重要影响。一般而言[16]，岩质整体滑坡形成的堰塞体稳定性好，溃决时溃口发展速度慢、溃决程度低；少量崩塌块石堆积形成的堰塞体可达成入流和渗流平衡而避免溃决；崩滑碎屑流形成的堰塞体易于快速溃决，危害性更大。

我国学者Fan等[15，22]基于崩滑发生区域的地质条件、堰塞体的地貌特征和材料物理力学特性，根据崩滑堰塞体的颗粒分布特征将2008年“5·12”汶川地震形成的堰塞体分为3类：第1类，堰塞体由于深部岩体的滑塌导致，内部结构可分为2层或3层，底部由较为完整的岩层构成，顶部为岩石碎屑和细颗粒土体，如图2(a)所示；第2类，堰塞体由边坡崩塌滑落的石块形成，内部结构可分为两层，底部为高度破碎的岩石碎屑和细颗粒土体，顶部由颗粒较大的碎石组成，如图2(b)所示；第3类，堰塞体由远程高速滑坡的碎屑流构成，由于滑坡体经过较长的移动，堰塞体颗粒较为松散且含有较多细颗粒，如图2(c)所示。

图2 基于颗粒分布特征的堰塞体分类Fig.2 Classification of landslide dams based on particle size distribution characteristics

3 崩滑堰塞体稳定性快速评价方法

崩滑堰塞体形成后，其在外荷载作用下的稳定性是揭示堰塞湖孕灾机理的核心。1999年，Casagli和Ermini[23]最早提出了堆积指标法(BI)来判断堰塞体的稳定性，并改进提出了无量纲堆积体指标法(DBI)[4]。随后，国内外学者基于已溃和未溃的堰塞湖资料，采用统计学的方法提出了一系列崩滑堰塞体稳定性的快速评价方法[24]，如Korup[25]、Dong等[26]和Stefanelli等[27]提出的评价方法，这些评价方法大多基于堰塞体的地貌学指标和堰塞湖的水动力条件，未考虑堰塞体的颗粒组成。目前国内外常用的堰塞体稳定性快速评价方法如表1所列。

表1 国内外常用堰塞体稳定性快速评价方法对比Tab.1 Comparison of rapid evaluation methods for stability of landslide dams at home and abroad

由于堰塞湖的成因各异，堰塞体的颗粒分布特征存在巨大差异，例如主要由大块石组成的堰塞体较土质堰塞体具有更好的稳定性，但目前尚缺乏能合理考虑堰塞体颗粒组成的评价方法。作者通过考虑堰塞体的形态特征、颗粒组成及被阻塞河道的水动力条件，采用逻辑回归的数值计算方法，建立了一套新的堰塞体稳定性快速评价方法[28]。根据可获取的堰塞体颗粒组成信息的多寡，该快速评价方法又可分为精细化和简化方法。

3.1 考虑颗粒组成的堰塞体稳定性精细化快速评价方法

若堰塞体材料的颗分情况已知，则选择高度、宽度(即顺河向距离)和体积来表征堰塞体的形态特征，选择堰塞湖流域面积来表征被阻塞河道的水动力特征，选择d90，d60，d30和d5等代表粒径表征颗粒组成特征，Shan等建立了堰塞体稳定性精细化快速评价方法，数学表达式如下[28]：

LsIVAS=-0.264lgI+1.166lgVd

-1.551lgAb-0.168lgSd-4.847

(1)

式中：Ls(IVAS)为堰塞体稳定性精细化评价指标，当Ls(IVAS)>0时，堰塞体是稳定的，当Ls(IVAS)<0时，堰塞体是不稳定的；I为堰塞体高度与宽度的比值，即I=Hd/W；Sd为精细化颗粒组成指标，Sd=(d90-d60)/(d30-d5)。

3.2 考虑颗粒组成的堰塞体稳定性简化快速评价方法

若堰塞体材料的颗粒特征仅为定性描述，则同样选择高度、宽度和体积来表征堰塞体的形态特征，选择堰塞湖流域面积来表征被阻塞河道的水动力特征，选择颗粒特征参数来表征颗粒组成特征，作者建立了堰塞体稳定性简化快速评价方法，数学表达式如下[28]：

LsIVAM=-0.198lgI+1.387lgVd

-1.432lgAb+4.169Mi-8.674

(2)

式中：Ls(IVAM)为堰塞体稳定性简化评价指标，当Ls(IVAM)>0时，堰塞体是稳定的；当Ls(IVAM)<0时，堰塞体是不稳定的。Mi为简化颗粒组成指标，当堰塞体材料以块石为主时(粒径>200 mm的颗粒重量超过50%，且粒径>2 mm的颗粒重量超过80%)，Mi取0.75～1.00；当堰塞体材料为块石加粗粒土时(粒径>200 mm的颗粒重量不超过50%，粒径>2 mm的颗粒重量超过80%)，Mi取0.50～0.75；当堰塞体材料以粗粒土为主时(粒径>2 mm的颗粒重量占20%～80%)，Mi取0.25～0.50；当堰塞体材料以细粒土为主时(粒径≤2 mm的颗粒重量超过80%)，Mi取 0～0.25；当颗粒含量位于某一取值范围之内，Mi采用线性插值法进行计算。

4 崩滑堰塞湖溃决机理

近年来，由于崩滑堰塞湖案例频发，国内外学者们围绕堰塞体的材料冲蚀特性、溃口形态和流量的演化规律开展了一系列的研究[29-30]。主要的研究手段包括：现场测试、小尺度模型试验和离心模型试验。

2008年“5·12”汶川地震形成的唐家山堰塞湖和2018年“11·03”白格堰塞湖都是在人工开挖泄流槽的情况下发生了溃决，科技人员通过各种手段获取了现场实测资料[11，31]，为崩滑堰塞湖溃决机理的研究提供了第一手的宝贵资料。以“11·03”白格堰塞湖泄流过程为例，可将溃决过程分为3个阶段[32]：

(1) 均匀冲蚀阶段。在冲蚀的初期阶段，由于溃口(初始泄流槽)水位较低、流速较小，仅堰塞体表层细颗粒被带走，溃口未发生明显下切现象，漫溢水流表现出均匀冲蚀的特征，此阶段溃口出流量小于来流量，堰塞湖水位继续抬升，如图3(a)所示。

(2) 溯源冲蚀阶段。随着上游水位的不断抬升，流速持续增大，由于堰塞体下游坡脚处的流速最大，因而率先发生冲蚀而导致结构性破坏，随后堰塞体发生溯源冲蚀直至达到上游坡顶部，如图3(b)所示。

(3) 沿程冲蚀阶段，当溯源冲蚀发展到溃口顶部后，由于溃口底高程的下降，水头突然增大，导致流量猛增直至峰值流量，溃口不断下切和展宽，并伴随溃口边坡的失稳，由于溃口流量大于来流量，堰塞湖水位不断下降直至漫溢水流无法冲蚀堰塞体材料为止，如图3(c)所示。

图3 2018年“11·03”白格堰塞湖溃决过程Fig.3 Breach process of the Baige landslide dammed lake at November 3，2018

另外，国内外学者开展了一系列针对堰塞湖溃决机理的小尺度模型试验和离心模型试验[33-37]，通过对溃口在横断面和纵断面的演化过程，可以总结出与现场观测资料类似的堰塞体溃决机理如下：堰塞体顶部溃口的向下冲蚀和下游坡面的溯源冲蚀→堰塞体下游坡角减小→溯源冲蚀发展至上游坡顶部→溃口快速下切和展宽→溃口边坡失稳(破坏面为平面)→溃口冲蚀结束。

值得一提的是，与均质土坝溃决机理不同，堰塞体溃决过程中下游坝坡逐渐变缓，且溃口的最终深度受堰塞体颗粒分布特征的影响。

5 崩滑堰塞湖溃决过程数值模拟方法

5.1 模型基本假设

基于崩滑堰塞湖的溃决机理，考虑堰塞体颗粒分布特征，在堰塞湖溃决过程数值模拟中采用以下假设：① 为简便起见，将堰塞体的横断面和纵断面分别视为倒梯形和正梯形；② 基于崩滑类型，结合地质调查获得的堰塞体颗粒分布特征，对崩滑堰塞体进行分层，且各层呈水平状分布；③ 各层顶部和底部的冲蚀系数由原位试验或经验公式确定，且各层的冲蚀系数由上至下呈线性变化；④ 在堰塞体溃决过程中，溃口边坡坡角在失稳前一直保持不变；⑤ 由于各层土体力学性质不同，溃口边坡失稳后的坡角由溃口所在位置处的土体特性确定；⑥ 溃口边坡失稳时沿平面滑动破坏，且滑落块体瞬间被冲走；⑦ 下游坡的溯源冲蚀速率等于溃口底部的向下冲蚀速率除以该时刻下游坡的坡比(垂直/水平)；⑧ 下游坡坡脚处的冲蚀深度为0，冲蚀深度沿下游坡脚向上至下游坡顶部呈线性增大，根据溃口底部和溯源冲蚀的深度确定下游坡坡角的减小量。

5.2 数值模拟方法

基于数值模拟的假设，建立了一个可描述崩滑堰塞湖溃决过程的数值模拟方法。该模拟方法主要包括3个组成部分[32，38]：水动力模块、材料冲蚀模块和溃口发展模块。采用按时间步长迭代的计算方法模拟崩滑堰塞湖溃决过程中的水土耦合作用，输入初始参数，设置计算时长tc和时间步长Δt，计算流程如图4所示。图4中，zb为溃口底部高程，zs为堰塞湖水位，t为时间，kd为冲蚀系数，τb为顶部溃口底床水流剪应力，τc为土颗粒临界剪应力，m为溃口边坡系数(水平/垂直)，H为堰塞体顶部溃口处水深(H=zs-zb)，Qb为溃口流量，Qin为入流量，As为堰塞湖湖面面积，B为溃口顶宽，b为溃口底宽，nloc为溃口位置参数(nloc=1表示溃口只能向一侧发展，nloc=2表示溃口可向两侧发展)，Fd为驱动力，Fr为抵抗力，α为失稳后溃口边坡的坡角，β为溃口边坡坡角，φ为土体的内摩擦角，C为土体的黏聚力，γs土体的容重，Hs为溃口边坡高度。

图4 崩滑堰塞湖溃决过程数值模拟计算流程Fig.4 Flow chart of numerical simulation of landslide dammed lake breach process

6 我国21世纪典型崩滑堰塞湖案例分析

为了对本文提出的崩滑堰塞湖形成-孕灾-致灾模拟方法的可行性进行验证，针对不同的崩滑方式，分别选择2008年“5·12”汶川地震形成的唐家山堰塞湖[5，11]、小岗剑堰塞湖[5，39-40]和2018年“11·03”白格堰塞湖[12，41]等3个典型崩滑堰塞湖，开展堰塞体稳定性快速评价和堰塞湖溃决过程模拟。

按照崩滑堰塞湖的形成机理和地质调查资料获取的堰塞体颗粒分布特征，根据前文2.2中的分类标准，唐家山属于第1类堰塞体，小岗剑属于第2类堰塞体，“11·03”白格属于第3类堰塞体。

6.1 典型堰塞体稳定性快速评价

分别采用作者提出的考虑颗粒组成的堰塞体稳定性精细化(式(1) )和简化快速评价方法(式(2) )对3个崩滑堰塞体的稳定性进行评价，并选择国内外常用的堰塞体稳定性评价方法，如Casagli和Ermini[23]提出的BI法、Ermini和Casagli[4]提出的DBI法、Korup[25]提出的Is法和Ia法、Dong等[26]提出的Ls(AHWL)和Ls(AHV)法、Stefanelli等[27]提出的HDSI法，与作者提出的方法进行比较。堰塞体稳定性快速评价计算参数如表2所列，各种评价方法的对比情况如表3所列。

表2 各评价方法计算参数Tab.2 The calculation parameters for each evaluation method

表3 各评价方法计算结果对比Tab.3 Comparison of calculation results of each evaluation method

调研结果表明，唐家山、小岗剑和“11·03”白格堰塞湖均为高危堰塞湖，若不开挖泄流槽也会自然溃决，因此其堰塞体均为不稳定堰塞体。计算结果表明：本文提出的Ls(IVAS)法和Ls(IVAM)对3个堰塞体的稳定性给出了准确的评价结论，同样，DBI法、Ia法、Ls(AHWL)法和Ls(AHV)法也给出了准确的评价结论；但是，BI法对唐家山和“11·03”白格堰塞体给出了错误的评价结论，Is法无法对3个堰塞体的稳定性给出明确的评价结论，HDSI无法对唐家山堰塞体的稳定性给出明确的评价结论。

6.2 典型堰塞湖溃决过程分析

唐家山堰塞湖通过开挖泄流槽的方式，于2008年6月10日泄流；小岗剑堰塞湖通过人工爆破的方式，于2008年6月12日泄流。“11·03”白格堰塞湖亦采用人工开挖泄流槽的方式，于2018年11月12日泄流。地质调查资料显示，唐家山堰塞体自上而下大致可分为3层，顶部以碎石土为主，第2层为块碎砾石土，第3层为似层状碎裂岩体；小岗剑堰塞体主要以孤石和块石为主，含少量碎石，细粒土填充于块石骨架之间，且顶部含有较多大粒径石块，自上而下大体可分为两层；白格堰塞体是由高位滑坡形成，滑坡土体经过长距离的运移颗粒发生了崩解，以砂砾石夹碎石土为主，整体颗粒较小，自上而下大体可分为两层。唐家山、小岗剑和“11·03”白格堰塞湖的溃决过程计算参数如表4和表5所列。

表4 我国3个典型堰塞湖物理力学参数Tab.4 Physical and mechanical parameters of three typical landslide dammed lakes in China

表5 我国3个典型堰塞体各层计算参数Tab.5 Calculation parameters of each layer of three typical landslide dams in China

为了验证模型的合理性，提取计算结果中的溃口峰值流量(Qp)、溃口最终顶宽(Bf)、溃口最终底宽(bf)、溃口最终深度(Df)以及溃口峰值流量出现时刻(Tp)等堰塞湖溃决特征参数与实测值进行比对(见表6)。

由表6可以看出：唐家山堰塞湖溃决时的溃口峰值流量、溃口最终顶宽、底宽和深度的计算误差控制在±5%以内，峰值流量出现时刻的计算误差控制在±10%以内；小岗剑堰塞湖溃决时的溃口峰值流量的计算误差控制在±5%以内，溃口最终底宽和深度、峰值流量出现时刻的计算误差控制在±10%以内；白格“11·03”堰塞湖溃决时的溃口峰值流量、溃口最终顶宽和峰值流量出现时刻的计算误差控制在±5%以内，溃口最终深度的计算误差控制在±15%以内，仅溃口最终底宽的误差较大，超过±30%。综上可知，崩滑堰塞湖溃决过程数值模拟方法可较好地反馈分析3个堰塞湖的溃决过程，比对结果验证了数值模拟方法的合理性。

表6 我国3个典型堰塞湖溃决特征参数计算值与实测值比较Tab.6 Comparison of the calculated and measured breach characteristic parameters of three typical landslide dammed lakes in China

7 结论与建议

本文围绕崩滑堰塞湖形成机理、堰塞体稳定性快速评价方法、堰塞湖溃决机理和溃决过程数值模拟方法进行了系统的阐述，介绍了作者提出的崩滑堰塞湖形成-孕灾-致灾过程数值模拟方法并验证了其合理性。主要结论和建议如下。

(1) 崩滑堰塞湖风险评估对于应急抢险具有重要的指导意义，可围绕其形成-孕灾-致灾的灾害链开展，其中的重要环节就是对堰塞体稳定性的快速评价和对堰塞湖溃决过程的模拟。研究表明，崩滑堰塞体的形态特征、颗粒组成及被阻塞河道的水动力条件是其稳定性的决定因素；另外，颗粒分布特征决定了堰塞体在溃决洪水作用下的材料冲蚀特性和溃口演化规律。

(2) 本文介绍的崩滑堰塞湖形成-孕灾-致灾模拟方法以崩滑堰塞湖形成机理为切入点，重点分析堰塞体材料颗粒组成和分布特征，在此基础上通过考虑堰塞体的形态特征、颗粒组成及被阻塞河道的水动力条件，建立了其稳定性快速评价方法；通过颗粒分布特征将堰塞体分层，建立了基于水动力模块、材料冲蚀模块和溃口发展模块的堰塞湖溃决过程模型，并通过典型案例验证了其合理性。

(3) 崩滑堰塞湖的灾害链牵涉到材料和荷载不确定性问题，以及复杂的水土耦合问题，应在材料测试设备与方法、冲蚀试验系统与技术和溃决模拟理论与算法等方面开展深入研究，更加科学地评估其形成-孕灾-致灾过程。