强震区沟道堰塞体失稳模式及其动力学特征

Failure Modes and Dynamic Characteristics of the Landslide Dams in Strong Earthquake Area

摘要

强震触发的大量崩塌滑坡所形成的松散固体物质堆积于泥石流沟道，容易形成天然堰塞体，在强降雨和上游流体的冲蚀下极易失稳形成溃决型泥石流.采用自制泥石流试验水槽，通过控制坝体颗粒组成和水动力条件，实施了12组堰塞体失稳模拟试验，获取了堰塞体的破坏过程、溃口流量和相关力学参数的演化特征. 将堰塞体划分为漫顶破坏、滑面破坏、管涌破坏3种不同失稳模式，并结合堰塞体的颗粒组成结构分析了失稳机理及特征，通过动力学过程分析分别建立了不同失稳模式下的堰塞体稳定性判别式. 研究成果对于溃决型泥石流防治工程的规划设计以及提高泥石流防灾减灾水平具有重要意义.

关键词

堰塞体 / 失稳模式 / 动力学特征 / 临界条件 / 滑坡

Key words

landslide dam / failure mode / dynamic characteristic / critical condition / landslides

0 引言

汶川地震导致西南山区地质地貌更加脆弱，大量的崩塌、滑坡及松散固体物质堆积在沟道中，为震后泥石流暴发提供了充足的物源条件，泥石流成为强震区主要的次生地质灾害类型（常鸣等， 2019； Fan et al.， 2019；余斌等， 2020）. 当泥石流沟道中含有一个或多个堰塞体，在强降雨及上游汇流作用下极易失稳形成溃决型泥石流，并引发流量放大效应，将对沟口的人民生命财产安全构成巨大威胁，造成的风险损失远高于一般的泥石流灾害（胡卸文等， 2016； 崔鹏等， 2021）. 而沟道中的堰塞体失稳模式直接决定了溃决型泥石流的灾变启动机制及运动过程. 若能预测泥石流沟道中天然堰塞体的失稳模式及稳定性，就能更合理的计算溃决型泥石流的流速、流量等特征参数. 因此对泥石流沟道中堰塞体的失稳模式和动力学特征进行研究，为溃决型泥石流灾害的防治提供技术支持，具有重要理论价值和现实意义.

近年来，各国学者针对堰塞体溃坝的研究主要集中在堰塞体的基本特征及失稳过程（Costa et al.， 1988； Dong et al.， 2009； 单熠博等， 2020），以及预测堰塞体稳定性（Fan et al.， 2012； Zhao et al.， 2019； 石振明等， 2021）. 在理论研究的基础上，众多学者也进行了溃坝模型实验. Zhou et al.（2015）通过实验发现漫顶破坏是堰塞体失稳的主要模式，而管涌和侵蚀作用一定程度上也会诱发堰塞体失稳. Zhao et al.（2015）通过模拟漫顶溃坝的数值实验后得出斜坡侵蚀是水位漫顶后主要的破坏形式. Zhong et al.（2017）以唐家山滑坡坝失稳为实例开展缩尺试验，认为侵蚀和间歇性作用是堰塞体失稳主要诱因. 刘定竺等（2017）通过物理实验发现堰塞体溃口展宽随溢流的发展而呈现不同形式，以溯源侵蚀的结束为时间节点，溯源侵蚀期间为近等宽阶段，溯源侵蚀结束后为弯曲阶段. 付建康等（2018）通过室内物理实验发现堰塞体漫顶溃坝主要由渗流、漫顶、冲刷和溃决4个阶段组成. 舒安平等（2019）通过对堰塞体溃决过程模拟研究，提出堰塞体溃决存在局部型和整体型两种模式. 除了室内物理试验，很多国内外学者也开展了不同尺度的野外模型试验（Hanson et al.， 2005； Morris et al.， 2007； Zhang et al.，2009）. 然而，目前研究主要是针对单一堰塞体开展失稳机理研究，尚未全面考虑到不同失稳模式下堰塞体所受动力学特征的差异性. 因此综合物理模型试验、动力学原理分析等手段，针对强震区溃决型泥石流沟道中不同失稳模式下的堰塞体分别提出临界判别公式十分必要.

本文通过控制颗粒结构及水动力条件两个因素，实施了12组堰塞体失稳物理模型试验，根据堰塞体的破坏过程，将失稳模式划分为3种模式：漫顶破坏、滑面破坏及管涌破坏；并总结了堰塞体的破坏模式和过程，获取了溃口流量和相关动力学参数的演化特征；结合室内模型试验与动力学特征分析，建立了强震区沟道堰塞体在不同失稳模式下的稳定性判别模型. 研究成果对溃决型泥石流灾害治理工程的规划设计及提高泥石流监测预警水平具有重要指导意义.

1 实验概述

1.1　实验模型及装置

堰塞体失稳物理模型实验平台由实验水槽模型、测试系统、实验物料等构成，见图1.

实验装置：实验装置由水箱、试验水槽及堆积平台构成. 水箱底部为0.7 m×0.7 m的正方形，高1.2 m，总容积0.588 m³. 试验水槽长2.5 m，宽0.4 m，高0.6 m，水槽可根据需要变动坡度，变坡范围5°~25°. 为了在实验过程中便于观测坝体溃决情况，槽体两侧为透明的钢化玻璃，底部为带花纹的钢化板，能提供足够的承重力. 堆积平台长1.5 m，宽1.5 m.

测试系统：由含水率传感器、高速摄像机等构成. 其中含水率传感器用于测量坝体内部含水率变化，传感器分别埋设在坝高8 cm（3处，编号3，4，5）、15 cm处（2处，编号1，2）；高速摄像机主要记录堰塞体失稳实验的全过程影像资料，分别布置在坝体侧面及正面.

堰塞体物料：选取汶川县羊岭沟①号堰塞体作为坝体材料，再通过配制形成试验所需物料级配，图2为羊岭沟堰塞体的分布情况.

1.2　模型相似分析

几何相似方面：选取羊岭沟①号堰塞体为典型案例，按照几何相似比1︰50进行缩尺模型构建，根据该堰塞体实际尺寸和特征，最终确定坝体长度为0.8 m，宽度0.4 m. 坝体迎水面坡度45°，背水面坡度30°等角度参数保持不变.

材料相似方面：考虑容重、级配、孔隙率与原型材料一致，即1︰1；粒径比为1︰50.

边界条件相似方面：水槽两侧为光滑透明钢化玻璃，底部为花纹钢板，在坝体与水槽接触部分均匀涂抹黏土，避免渗流现象，使实验更贴近实际情况.

1.3　模型实验设计

Zhu et al.（2020）通过堰塞体失稳模型实验，发现堰塞体内部渗流是控制堰塞体破坏模式和过程的主要因素之一，而渗流条件主要受堰塞体的粒径分布控制（Jiang et al.，2018； Zhu et al.， 2019）. Casagli et al.（2003）将堰塞体颗粒划分为基质支撑型和颗粒支撑型，两种类型的主要区别是粗颗粒和细颗粒的含量，但是粗颗粒和细颗粒是一个相对概念，和水动力条件有关. 因此，开展的堰塞体失稳物理模型实验的两个控制因子为水槽坡度及堰塞体级配：（1）水槽设计坡度依次为5°、9°、11°、13°、15°、17°；（2）6种实验沙级配的中值粒径分别为1.0 mm、1.5 mm、2.5 mm、3.3 mm、4.3 mm、6.5 mm，实验方案见表1. 实验土颗粒级配如图3所示.

1.4　实验步骤

先将实验土样进行晒干处理，保证初始含水率为0，根据堰塞体的实际颗粒组成情况，用筛析机筛分出不同粒径的试验土，然后按重量比例均匀混合在一起，达到堰塞坝尺寸要求所需用土约110 kg. 将实验物料按预定尺寸堆积，并在固定位置埋设传感器. 打开阀门至恒定流量0.25 L/s，与此同时开启含水率传感器和两台高速摄像机进行观测，直到坝体破坏，剩余坝体稳定为止. 实验结束后，对实验样品进行称重、烘干处理，并完成颗粒分析试验，最终对全部实验资料进行整理分析.

2 堰塞体失稳模式分析

以上所述的12组模型试验中，堰塞体的失稳模式可以划分为3种. 其中A~E、H表现为漫顶破坏模式，F、G表现为滑面破坏模式，I、J、K、L表现为管涌破坏模式，详细的描述如下.

2.1　漫顶破坏模式

漫顶破坏模式可以概括为漫顶溢流-溃口连通-快速下切-坝体稳定4个阶段. 上游水流启动，堰塞体上游水位逐渐升高至漫顶状态，水流漫顶溢流开始堰塞体冲刷堰塞体下游坡面，最开始在坝顶和下游坡面的过渡区冲刷，如图4a；过渡区被冲刷后，向上游方向成溯源侵蚀，向下游形成冲刷下切，形成小型冲沟，并以下切为主，横向扩展较小，如图4b；一旦溯源冲刷至坝后水体处，溃口连通，堰塞体后水体开始下泄，溃口此时向两个方向同时快速扩展，即横向展宽，纵向下切，同时溃口两侧土体发生破环，出现陡坎和坝坡失稳现象，如图4c；随着坝后水位的降低，流速流量逐渐减小，最终流体的剪切应力逐渐减小趋近于坝体颗粒的抗剪强度，坝体重新稳定，如图4d.

在实验过程中，记录了坝体内部含水率和流量的变化情况. 图5a为含水率传感器的变化情况，含水率变化的顺序是3→4→1→2，5号传感器无明显变化. 说明坝体内水位缓慢升高，初始渗流以水平为主（3号和4号传感器率先变化）. 随着水位的升高，坝体上部的含水率传感器开始变化（1和2号）. 而5号传感器没有任何变化，说明整个试验过程中坝体内部渗流都未到达5号传感器所在位置（堰塞体底部前端）. 在t=400 s左右，所有含水率传感器读数均下降至初始状态，此时堰塞体已经发生了失稳.

图5b为流量的变化情况，漫顶溢流阶段，水流溢流冲刷堰塞体坡面，坝后水位无变化，坡面侵蚀带走了部分堰塞体表面颗粒，此时下泄流量略多于入流量；溃口连通后，堰塞体后水体开始下泄，流量陡增；溃口处断面流量到达峰值之后，流量逐渐减小，直至坝体稳定. 整个试验过程中，溃口流量并非连续增加或连续下降的，而是呈现出波动性. 由于水流的侧蚀和下蚀以及坝坡失稳堵塞沟道后再次溃决，导致的溃口形状突变会引起溃口流量的突变，所以往往一次较高流量峰值伴随着一次流量的下降.

2.2　滑面破坏模式

滑面破坏模式可以概括为渗流浸润-前兆滑动-多级滑动-坝体稳定4个阶段. 此类天然堰塞体透水系数较大，且强度极弱，当坝体内水位上升后，使得坡体自重增加，同时在坡体的饱和部分，因渗透水的浮力作用，颗粒间的摩擦阻力降低，如图6a；坡体前半部分沿着某一薄弱滑动面滑动，滑动面呈近似圆弧形，如图6b；当坡体前半部分滑动后，后半部分因失去支撑从而沿着坝顶向下延伸的一个圆弧形滑动面滑移，并堆积于坡脚，总体呈现牵引式滑动，如图6c；多次滑动后降低了堰塞体的整体高度，坝后蓄水倾泻而出，与失稳后的堰塞体混合形成溃决型泥石流，如图6d.

在实验过程中，记录了坝体内部含水量和流量变化. 图7a为含水率传感器的变化情况，含水率变化的顺序是3→4→5→1→2. 说明坝后蓄水期间部分细颗粒沿渗流路径运移，在坝体下部迅速形成了一个完整的渗流通道. 随着水位的升高，坝体上部的含水率传感器开始变化（1号和2号）. 在t=500 s左右，所有传感器位置的含水率同时发生了不同幅度的下降，而后又迅速恢复为正常水平，紧接着在t=515 s左右再次下降后保持稳定. 从实验过程可以看出，这与两次滑坡发生的时间恰好吻合.

图7b为流量的变化情况，由于滑面破坏发生迅速，从开始溃坝时刻开始，流量在漫顶后3 s就到达峰值，紧接着逐渐减小，直至和入流量相等. 流量在衰减过程中出现小幅度的波动，是由于溃坝过程中仍会发生小规模边坡失稳导致堵塞-溃决现象.

2.3　管涌破坏模式

管涌破坏模式可以概括为管涌-小型崩滑-溢流侵蚀-坝体稳定4个阶段. 随着上游来流的不断涌入，坝后水位逐渐升高，导致坝体内孔隙水压力增加，细小的可动颗粒在渗流力的作用下形成管涌通道，如图8a；而后，堰塞体管涌出口处发生小型崩滑，降低了堰塞体的整体稳定性，如图8b；坝后水位上升至坝顶后出现漫顶溢流，堰塞体在坡面水流的冲刷下形成溃口，同时伴随着强烈的侧向展宽和下切侵蚀，溃口流量快速增加并到达峰值，如图8c；随着堰塞体后水量逐渐减小，断面流量到达峰值之后，流量迅速减小直至趋于稳定，如图8d.

在实验过程中，记录了坝体内部含水率和流量的变化情况. 图9a为含水率传感器的变化情况，含水率变化的顺序是3→1→4→5→2. 随着各点含水率相继达到峰值后，含水率依次急剧. t=607 s时，坝体下游表面发生小型崩滑，孔隙水压力和含水率急剧下降. 堰塞坝溃决后，所有含水率传感器读数均下降至初始状态，此时堰塞体已经发生了失稳.

图9b为堰塞体溃决过程中流量的变化情况，当堰塞体坡面发生崩滑时，管涌通道迅速扩大，流量逐渐增加，堰塞体稳定性迅速降低导致堰塞体失稳，此时流量迅速增大，达到峰值流量，而后逐渐降低直至残余坝体稳定. 与漫顶破坏模式相同，整个试验过程中，溃口流量并非连续增加或连续下降的，而是呈现出波动性.

2.4　堰塞体结构对失稳模式和过程的影响

堰塞体的渗流是影响堰塞体破坏模式和过程的重要因素. 为了充分考虑渗流对堰塞体破坏模式和过程的影响，将入库流量保持在较小的情况下，通过对上述12组试验的破坏模式和破坏过程的初步分析，发现改变颗粒粒径分布和水动力条件对堰塞体的破坏模式有显著影响.

依据堰塞坝体材料的级配特点，可将堰塞坝体材料分为基质支撑型和颗粒支撑型两种结构形式. 基质支撑型的坝体材料，粗颗粒之间的空隙由细颗粒填充，粗颗粒相互不接触，密实度较高，渗透性较低；颗粒支撑型的坝体材料，粗颗粒的相互接触形成骨架，密实度较低，渗透性较高，易发生管涌. 一般将2 mm作为划分粗颗粒和细颗粒的界限（Cui et al.， 2017）. 土体内部稳定性取决于其粗颗粒能否防止细颗粒的流失. 当细粒材料的含量在一定范围内时，粗颗粒才可能形成骨架. 如果细粒物质的含量较多，粗粒物质在细粒物质形成的基质中流动，在这种情况下，堰塞体的稳定性主要由细颗粒物质控制.

根据试验结果和堰塞体结构分类，进一步分析堰塞体粒径分布对破坏模式和过程的影响. 实验I~L中粗颗粒的含量均大于50%，为颗粒支撑型，因此均形成了渗流通道，表现为管涌破坏. 实验A~H则属于基质支撑型，除实验F、G为滑面破坏，其余均为漫顶破坏，是由于实验F、G水力梯度更大，容易形成渗流通道，坝体在土体自重力和内部渗流力的影响下产生崩滑.

堰塞体失稳过程中常常伴随着溃口两侧的坝坡失稳，坝坡的抗剪强度可由莫尔-库仑方程表示：

{\tau }_{\mathrm{f}} = c^{\mathrm{\text{'}}} + ( \sigma - \mu ) \mathrm{t} \mathrm{a} \mathrm{n} \phi \mathrm{\text{'}},(1)

式中： τ _f为抗剪强度（kPa）；c´为有效内聚力（kPa）；σ为土体总应力（kPa）；μ为孔隙水压力（kPa）；φ´为有效内摩擦角（°）.

随着细颗粒含量的降低，坝坡强度也会随之降低（Lin and Takahashi， 2012）. 用孔隙内空气压力（μ _a）与孔隙水压力（μ _w）之差表示基质吸力（ψ）. 随着基质吸力的增加，导致的抗剪强度的增加可用角度φ ^b来描述. 将这一效应代入莫尔-库仑方程，可以得到：

{\tau }_{\mathrm{f}} = c^{\mathrm{\text{'}}} + \left(\sigma -{\mu }_{\mathrm{a}}\right) \mathrm{t} \mathrm{a} \mathrm{n} {\phi }^{\mathrm{\text{'}}} + ( {\mu }_{\mathrm{a}} - {\mu }_{\mathrm{w}} ) \mathrm{t} \mathrm{a} \mathrm{n} {\phi }^{\mathrm{b}},(2)

因此，细粒颗粒含量越高，坝体强度越高. 考虑到坝体局部区域在溃决过程中存在非饱和物质，在溃口两侧更容易形成较大规模的坝坡失稳. 根据实验结果，也可以得出随着细粒颗粒含量的逐渐增加，坝坡失稳次数也随之增加. 溃口两侧坝坡失稳堆积于溃口中并堵塞沟道，导致溃坝过程中出现流量多峰现象.

3 堰塞体失稳判别式构建

3.1　漫顶破坏模式

通过试验观测及野外调查，堰塞体发生漫顶溢流之后，由于强烈的溯源侵蚀及冲刷下切，冲沟迅速扩宽，最后导致堰塞体失稳. 这一阶段的定量分析表明，发生漫顶溢流后，部分堰塞体材料在坡面水流作用下开始移动，启动颗粒物质受力如下：水下的重力W'，坡面的摩擦力f，坡面的支持力N，水流的拖拽力F _D及上扬力F_y，如图10所示.

颗粒所受水下重力为：

W \mathrm{’} = \frac{\pi }{\mathrm{6}} D^{\mathrm{3}} \left({\gamma }_{\mathrm{s}}-{\gamma }_{\mathrm{w}}\right),(3)

式中：D为颗粒物的直径（m），γ _S为堰塞体的重度（kN/m³）；γ_w 为水的重度（kN/m³）.

坡面物质对启动颗粒的支持力为：

N = \frac{\pi }{\mathrm{6}} D^{\mathrm{3}} \left({\gamma }_{\mathrm{s}}-{\gamma }_{\mathrm{w}}\right) \mathrm{c} \mathrm{o} \mathrm{s} \theta - F_{\mathrm{y}}.(4)

堰塞体表面颗粒一般比较松散，内聚力很低，因此不考虑颗粒间的内聚力. 坡面摩擦力为：

f = \left[\frac{\pi }{\mathrm{6}}{D}^{\mathrm{3}}\left({\gamma }_{\mathrm{s}}-{\gamma }_{\mathrm{w}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\theta -F_{\mathrm{y}}\right] \mathrm{t} \mathrm{a} \mathrm{n} \phi，(5)

式中： θ为堰塞体坡面倾角（°）；φ为堰塞体有效内摩擦角（°）.

拖拽力为：

F_{\mathrm{D}} = C_{\mathrm{D}} \pi D^{\mathrm{2}} \frac{{\gamma }_{\mathrm{w}}{v}_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}g},(6)

式中：C _D为阻力系数，取决于物体形状和雷诺数，参考杰缅季夫的实验，取0.712，无量纲；g为重力加速度，取9.8 m/s²；v_c 为漫顶流体流速（m/s），采用强震区泥石流临界启动流速公式（张惠惠， 2011）：

v_{\mathrm{c}} = F {\left(\mathrm{1}+\epsilon \right)}^{\frac{\mathrm{1}}{\mathrm{2}}} {\left(f\mathrm{c}\mathrm{o}\mathrm{s}\beta -\mathrm{s}\mathrm{i}\mathrm{n}\beta \right)}^{\frac{\mathrm{1}}{\mathrm{2}}} {\left(\frac{d}{k}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}} h^{\frac{\mathrm{1}}{\mathrm{6}}},(7)

F = \left\{\begin{array}{c}\mathrm{10.01},\mathrm{个别}\mathrm{启动}\\ \mathrm{10.7},\mathrm{少量}\mathrm{启动}\\ \mathrm{11.7},\mathrm{泥石}\mathrm{流形}\mathrm{成}\end{array}\right.

式中：v_c 为启动流速（m/s）；ε为床面泥沙颗粒相对荫暴度系数，取值为0~1，无量纲；f 为床面砂粒间的摩擦系数，取值为0.63，无量纲；β为床面与水平面的夹角（°）；d为颗粒粒径（m）；k为糙率尺寸，应考虑床沙的组成，根据床沙级配中的代表粒径或其倍数来确定（m）；h为启动水深（m）.

堰塞体扬压力为底部水的浮托及静水压力，可根据扬压力计算公式求解：

F_y = \frac{\mathrm{1}}{\mathrm{2}} {\gamma }_{\mathrm{w}} \left(\frac{\mathrm{2}}{\mathrm{3}}D+D\mathrm{s}\mathrm{i}\mathrm{n}\theta \right) \mathrm{c} \mathrm{o} \mathrm{s} \theta.(8)

将沿坡面向上的阻抗力与沿坡面向下的启动力比值定义为稳定性系数K _S，可求得稳定性系数：

\begin{array}{l}{K}_{\mathrm{S}}=\frac{{\tau }_{\mathrm{f}}}{{\tau }_{\mathrm{d}}}=\{[\frac{\pi }{\mathrm{6}}{D}^{\mathrm{3}}({\gamma }_{\mathrm{s}}-{\gamma }_{\mathrm{w}})\mathrm{c}\mathrm{o}\mathrm{s}\theta -\frac{\mathrm{1}}{\mathrm{2}}{\gamma }_{\mathrm{w}}(\frac{\mathrm{2}}{\mathrm{3}}D+\mathrm{s}\mathrm{i}\mathrm{n}\theta )\mathrm{c}\mathrm{o}\mathrm{s}\theta ]\mathrm{t}\mathrm{a}\mathrm{n}\phi +\frac{\mathrm{1}}{\mathrm{2}}{\gamma }_{\mathrm{w}}(\frac{\mathrm{2}}{\mathrm{3}}D+D\mathrm{s}\mathrm{i}\mathrm{n}\theta )\mathrm{c}\mathrm{o}\mathrm{s}\theta \mathrm{s}\mathrm{i}\mathrm{n}\theta \}\\ /[C_{\mathrm{D}}\mathrm{\pi }{D}^{\mathrm{2}}\frac{{\gamma }_{\mathrm{w}}{v}_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}\mathrm{g}}+\frac{\mathrm{\pi }}{\mathrm{6}}{D}^{\mathrm{3}}({\gamma }_{\mathrm{s}}-{\gamma }_{\mathrm{w}})\mathrm{s}\mathrm{i}\mathrm{n}\theta ]\end{array}.（9)

3.2　滑面破坏模式

主要发生在堰塞体坡度较陡，渗透系数较大时. 往往是由于强度降低而产生坝体崩滑，并逐渐向上游逐级扩展造成的. 坝体崩滑后会在坡脚下堆积，在来流作用下将会形成溃决型泥石流. 因此可以借鉴边坡问题的研究方法来分析堰塞体稳定性.

此类坝体细颗粒含量低，颗粒间粘聚力可忽略不计，将其视为无粘性土坡，假设堰塞体为各向同性的、完全浸水的，且不存在渗流作用，在坡面上取一土单元体，受力分析如图11所示.

因此堰塞体单元的下滑剪切力为：

{\tau }_{\mathrm{d}} = \gamma \text{'} \mathrm{s} \mathrm{i} \mathrm{n} f \alpha + {\gamma }_{\mathrm{w}} \mathrm{s} \mathrm{i} \mathrm{n} \alpha = {\gamma }_{\mathrm{s}\mathrm{a}\mathrm{t}} \mathrm{s} \mathrm{i} \mathrm{n} \alpha,(10)

式中： τ _d为堰塞体单元的下滑剪切力（kN/m³）；γ＇为单元堰塞体自重（kN/m³）；α为堰塞体背水面坡度（°）；γ _w为水的重度（kN/m³）；γ _sat为堰塞体饱和重度（kN/m³）.

阻止堰塞体下滑的力则为单元体与下部土体的抗剪力，等于单元体自重在坡面法向方向分力引起的摩擦力，即：

{\tau }_{\mathrm{f}} = \gamma \text{'} \mathrm{c} \mathrm{o} \mathrm{s} \alpha \mathrm{t} \mathrm{a} \mathrm{n} \phi,(11)

式中： τ _f为堰塞体单元的抗剪力（kN/m³）；φ为堰塞体内摩擦角（°）.

故堰塞体的稳定性系数K _S为：

K_{\mathrm{S}} = \frac{{\tau }_{\mathrm{f}}}{{\tau }_{\mathrm{d}}} = \frac{\gamma \text{'}\mathrm{t}\mathrm{a}\mathrm{n}\phi }{{\gamma }_{\mathrm{s}\mathrm{a}\mathrm{t}}\mathrm{t}\mathrm{a}\mathrm{n}\alpha },(12)

3.3　管涌破坏模式

管涌破坏一般发生在颗粒支撑型的堰塞体中. 骨架孔隙中的可动颗粒受到的力有：渗流水流的拖曳力、水下重力、与孔隙壁之间的摩擦力和颗粒之间的相互作用力. 对于堰塞体内部发生的管涌而言，只有渗透水对可动颗粒的拖曳力为起动力，其它力都为阻力. 且可动颗粒在孔隙壁中移动时，摩擦力是最大的阻力，而颗粒之间的碰撞力、电场力等都非常小且十分复杂. 假设渗流通道水平，用重力沿渗流通道向下的启动力抵消颗粒之间的碰撞力、电场力等阻力，阻止颗粒起动的力以管道之间的摩擦力为主.

孔隙中单个可动颗粒所受到的摩擦力为：

f = \frac{\pi }{\mathrm{6}} D^{\mathrm{3}} ( {\gamma }_{\mathrm{s}} - {\gamma }_{\mathrm{w}} ) \mathrm{t} \mathrm{a} \mathrm{n} \phi,(13)

式中： γ _s为堰塞体的重度（kN/m³）；γ _w为水的重度（kN/m³）；φ为堰塞体有效内摩擦角（°）；D为可动颗粒的粒径（m），取太沙基控制粒径d ₁₅为可动颗粒的最大粒径（吴梦喜，2019）.

发生渗流部位的单位体积堰塞体颗粒的渗流压力为：

P = {\gamma }_{\mathrm{w}} i,(14)

i = \frac{\mathrm{d}h}{\mathrm{d}l} = \frac{h_{\mathrm{w}}+l\mathrm{s}\mathrm{i}\mathrm{n}\eta }{l},(15)

式中： i为水力梯度，无量纲：h _w为坝后堰塞湖水深（m）；l为堰塞体长度（m）；η为堰塞体沟床坡度（°）.

仅由可动颗粒承受动水压力，单个颗粒所受到的渗透力为：

F_{\mathrm{P}} = {\gamma }_{\mathrm{w}} i \frac{\mathrm{\pi }{D}^{\mathrm{2}}}{S},(16)

S = \frac{\mathrm{6}\alpha (\mathrm{1}-n)}{d_{\theta k}},(17)

式中： S为单位体积堰塞体通道中可动颗粒的表面积（m^-1）；α为颗粒形状系数，无量纲，对于各种砂粒α=1.16~1.17，对于锐角颗粒α=1.50~1.67；n为孔隙率，无量纲；d_θk 为土体的等效粒径，一般计算方法为：

d_{\theta k} = \sqrt[]{\frac{\mathrm{2.13}k{(\mathrm{1}-n)}^{\mathrm{2}}}{\mathrm{g}{n}^{\mathrm{3}}}},(18)

式中： k为堰塞体的渗透系数（m/s）；n为孔隙率，无量纲；g为重力加速度，取9.81 m/s².

由于可动颗粒实际移动路线并非直线，为了更准确地描述可动颗粒在通道中的运动情况，引入弯曲度模型（陆高明，2016）. 堰塞体中弯曲度为可动颗粒实际移动路线长度与堰塞体长度的比值，弯曲度的计算公式为：

\begin{array}{l}\tau =\mathrm{0.75}+\frac{\mathrm{0.25}}{\sqrt[]{\mathrm{1}-\frac{\mathrm{4}}{\pi }(\mathrm{1}-n)}}+\\ \frac{\mathrm{\pi }-\mathrm{2}}{\mathrm{4}}\sqrt[]{\frac{\mathrm{1}-n}{\mathrm{\pi }}}\end{array}.(19)

引入弯曲度模型后，单位体积堰塞体渗流通道中可动颗粒的表面积为：

S = \frac{\mathrm{6}\alpha (\mathrm{1}-n)\tau }{d_{\theta k}}.(20)

联立上式得到渗透水流作用在单个颗粒上的渗透力为：

\begin{array}{l}{F}_{\mathrm{P}}=[{\gamma }_{\mathrm{w}}i\mathrm{\pi }{D}^{\mathrm{2}}\sqrt[]{\frac{\mathrm{2.13}k{(\mathrm{1}-n)}^{\mathrm{2}}}{\mathrm{g}{n}^{\mathrm{3}}}}]/\\ \{\mathrm{6}\alpha (\mathrm{1}-n)[\mathrm{0.75}+\\ \frac{\mathrm{0.25}}{\sqrt[]{\mathrm{1}-\frac{\mathrm{4}}{\mathrm{\pi }}(\mathrm{1}-n)}}+\frac{\mathrm{\pi }-\mathrm{2}}{\mathrm{4}}\sqrt[]{\frac{\mathrm{1}-n}{\mathrm{\pi }}}]\}\end{array} .(21)

结合式（13）、式（21），得到管涌破坏模式下的堰塞体稳定性系数K _S计算公式为：

\begin{array}{l}{K}_{\mathrm{S}}=\frac{f}{F_{\mathrm{P}}}=[{\gamma }_{\mathrm{w}}\frac{h_{\mathrm{w}}+l\mathrm{s}\mathrm{i}\mathrm{n}\eta }{l}\sqrt[]{\frac{\mathrm{2.13}k{(\mathrm{1}-n)}^{\mathrm{2}}}{\mathrm{g}{n}^{\mathrm{3}}}}]/\{\alpha \\ \left(\mathrm{1}-n\right)[\mathrm{0.75}+\frac{\mathrm{0.25}}{\sqrt[]{\mathrm{1}-\frac{\mathrm{4}}{\mathrm{\pi }}\left(\mathrm{1}-n\right)}}+\frac{\mathrm{\pi }-\mathrm{2}}{\mathrm{4}}\sqrt[]{\frac{\mathrm{1}-n}{\mathrm{\pi }}}]D\left({\gamma }_{\mathrm{s}}-{\gamma }_{\mathrm{w}}\right)\mathrm{t}\mathrm{a}\mathrm{n}\phi \}\end{array}.（22）

4 结论

本文通过控制坝体颗粒结构组成及水动力条件，进行了12组堰塞体失稳模拟实验，将堰塞体的破坏模式划分为3种，并提出了不同破坏模式下的堰塞体的稳定性系数计算公式，主要得到以下结论：

（1）通过实施12组溃坝模拟实验和成果分析，提出了堰塞体的3种失稳模式，即漫顶破坏模式、滑面破坏模式、管涌破坏模式. 其中，漫顶破坏模式概括为漫顶溢流-溃口连通-快速下切-坝体稳定4个阶段；滑面破坏模式概括为渗流浸润-前兆滑动-多级滑动-坝体稳定4个阶段；管涌破坏模式概括为管涌-小型崩滑-溢流侵蚀-坝体稳定4个阶段.

（2）根据堰塞体的颗粒级配组成将堰塞体结构分为基质支撑型和颗粒支撑型，不同结构类型决定了堰塞体失稳的启动机理和溃决过程，由基质支撑型材料组成的堰塞体多为漫顶溢流破坏，而颗粒支撑型材料组成的堰塞体均表现为管涌破坏；且细粒颗粒含量越高，坝体强度越高，坝坡失稳次数随之减小.

（3）基于堰塞体的3种破坏模式，总结了堰塞体的溃决机理及特征，通过动力学分析分别建立了不同破坏模式下的堰塞体稳定性判别式. 研究成果对于溃决型泥石流防治工程的规划设计以及提高泥石流防灾减灾水平具有重要意义.

参考文献

[1] Casagli, N., Ermini, L., Rosati, G.,2003. Determining Grain Size Distribution of Material Composing Landslide Dams in the Northern Apennine: Sampling and Processing Methods. Engineering Geology, 69 (1): 83-97.

[2] Chang, M.,Dou, X.Y., Tang, C., et al., 2019. Hazard Assessment of Typical Debris Flow Induced by Rainfall Intensity. Earth Science, 44(8):2794-2802 (in Chinese with English abstract).

[3] Costa, J. E., Schuster, R. L., 1988. The Formation and Failure of Natural Dams. Geological Society of America Bulletin, 100(7): 1054-1068. https://doi.org/10.1130/0016-7606(1988)100<1054:tfafon>2.3.co;2

[4] Cui P, Guo, J, 2021. Evolution Models, Risk Prevention and Control Countermeasures of the Valley Disaster Chain. Advanced Engineering Sciences, 53(3):5-18 (in Chinese with English abstract).

[5] Cui, Y. F., Zhou, X. J., Guo, C. X., 2017. Experimental Study on the Moving Characteristics of Fine Grains in Wide Grading Unconsolidated Soil under Heavy Rainfall. Journal of Mountain Science, 14(3): 417-431. https://doi.org/10.1007/s11629-016-4303-x

[6] Dong, J. J., Tung, Y. H., Chen, C. C., et al., 2009. Discriminant Analysis of the Geomorphic Characteristics and Stability of Landslide Dams. Geomorphology, 110(3/4): 162-171. https://doi.org/10.1016/j.geomorph. 2009. 04.004

[7] Fan, X.M., Scaringi, G., Korup, O.,et al., 2019. Earthquake‐Induced Chains of Geologic Hazards: Patterns, Mechanisms, and Impacts. Reviews of Geophysics, 57(2): 421-503. https://doi.org/10.1029/2018rg000626

[8] Fan, X., Tang, C. X., van Westen, C. J., et al., 2012. Simulating Dam-Breach Flood Scenarios of the Tangjiashan Landslide Dam Induced by the Wenchuan Earthquake. Natural Hazards and Earth System Sciences, 12(10): 3031-3044. https://doi.org/10.5194/nhess-12-3031-2012

[9] Fu, J.K., Luo, G., Hu, X.W., et al., 2018. Physical Model Experiment on Overtopping Overflow Failure of Landslide Dam. Journal of Jilin University(Earth Science Edition), 48(1): 203-212 (in Chinese with English abstract).

[10] Hanson, G.J., Cook, K.R., Hunt, S.L., et al.,2005. Physical Modeling of Overtopping Erosion and Breach Formation of Cohesive Embankments. Transactions of ASAE, 48(5): 1783-1794.

[11] Hu, X.W., Han, M., Liang, J.X., et al., 2016. Some Key Problems on Debris Flow in Wenchuan Earthquake Area. Journal of Southwest Jiaotong University, 51(2): 331-340 (in Chinese with English abstract).

[12] Jiang, X.G., Wei, Y.W., Wu, L.,et al., 2018. Experimental Investigation of Failure Modes and Breaching Characteristics of Natural Dams. Geomatics, Natural Hazards and Risk, 9(1): 33-48. https://doi.org/10.1080/19475705. 2017.1407367

[13] Lin, K., Takahashi, A., 2012. Strength Reduction of Cohesionless Soil Due to Internal Erosion Induced by One-Dimensional Upward Seepage Flow. Soils and Foundations, 52(4): 698-711. https://doi.org/10.1016/j.sandf. 2012. 07.010

[14] Liu, D.Z., Cui, P., Jiang, D.W., et al., 2017. Experimental Study on Breach Broadening Process of Landslide Dam. Science of Soil and Water Conservation, 15(6): 19-26 (in Chinese with English abstract).

[15] Liu, J.F., You, Y., Chen, X.C., et al., 2010. The Characteristics and Countermeasures of Dam-Breaking Debris Flow after Wenchuan Earthquake: A Case Study of the Tangfang Gully in Pingwu County,Sichuan Province. Advanced Engineering Sciences, 42(5): 68-75 (in Chinese with English abstract).

[16] Lu, G.M., Xia, Y.Y., Rui, R., 2016. Improvement of Piping's Critical Slope Ratio Equation Based on the Tortuosity and Chords Model. Journal of Wuhan University of Technology, 38(3): 41-47 (in Chinese with English abstract).

[17] Morris, M. W., Hassan, M. A. A. M., Vaskinn, K. A., 2007. Breach Formation: Field Test and Laboratory Experiments. Journal of Hydraulic Research, 45(sup1): 9-17. https://doi.org/10.1080/00221686.2007.9521828

[18] Shan, Y.B., Chen, S.S., Zhong, Q.M., 2020. A Rapid Evaluation Method of Landslide Dam Stability. Chinese Journal of Rock Mechanics and Engineering, 39(9):1847-1859 (in Chinese with English abstract).

[19] Shi, Z.M., Zhou, M.J., Peng, M., et al., 2021. Research Progress on the Mechanisms and Breaching Flood of Overtopping Failure of Landslide Dams Caused by Landslides and Avalanches. Chinese Journal of Rock Mechanics and Engineering,1-16(in Chinese with English abstract).

[20] Shu, A.P., Zhu, F.Y., Wang, S., et al., 2019. Starting Processes and Dynamic Characteristics of Dam-Break Debris Flow. Journal of Hydraulic Engineering, 50(6): 661-669 (in Chinese with English abstract).

[21] Wu, M.X., Gao, G.Y., Yang, J.X., et al., 2019. A method of Predicting Critical Gradient for Piping of Sand and Gravel Soils. Rock and Soil Mechanics, 40(3):861-870 (in Chinese with English abstract).

[22] Yu, B., Yang, L.Y., Liu, Q.H., et al., 2020. A Precise Prediction Model on Debris Flows Caused by Runoff Mechanism Based on Channel Width and Particle Size. Earth Science, 45(4):1447-1456 (in Chinese with English abstract).

[23] Zhang, H.H., 2011. An Study On Debris Flow Disaster Caused by Flood-Triggering Channel Accumulation in Post-Earthquake Meizoseismal Areas(Dissertation). Chengdu University of Technology, Chengdu: 24-35 (in Chinese with English abstract).

[24] Zhang, J. Y., Li, Y., Xuan, G. X., et al., 2009. Overtopping Breaching of Cohesive Homogeneous Earth Dam with Different Cohesive Strength. Science in China Series E: Technological Sciences, 52(10): 3024-3029. https://doi.org/10.1007/s11431-009-0275-1

[25] Zhao, T. L., Chen, S. S., Fu, C. J., et al., 2019. Centrifugal Model Tests and Numerical Simulations for Barrier Dam Break Due to Overtopping. Journal of Mountain Science, 16(3): 630-640. https://doi.org/10.1007/s11629-018-5024-0

[26] Zhao, W. Y., Chen, X. Q., You, Y., et al., 2015. Dam-Break Characteristics of Landslide Dams with Different Types of Open Channel Discharge Sections. Environmental Earth Sciences, 74(6): 5331-5340. https://doi.org/10.1007/s12665-015-4543-z

[27] Zhong, Q. M., Chen, S. S., Mei, S. A., et al., 2017. Numerical Simulation of Landslide Dam Breaching Due to Overtopping. Landslides, 15(6): 1183-1192. https://doi.org/10.1007/s10346-017-0935-3

[28] Zhou, G. G. D., Cui, P., Zhu, X. H., et al., 2015. A Preliminary Study of the Failure Mechanisms of Cascading Landslide Dams. International Journal of Sediment Research, 30(3): 223-234. https://doi.org/10.1016/j.ijsrc.2014.09.003

[29] Zhu, X. H., Peng, J. B., Jiang, C., et al., 2019. A Preliminary Study of the Failure Modes and Process of Landslide Dams Due to Upstream Flow. Water, 11(6): 1115. https://doi.org/10.3390/w11061115

[30] Zhu, X. H., Peng, J. B., Liu, B. X., et al., 2020. Influence of Textural Properties on the Failure Mode and Process of Landslide Dams. Engineering Geology, 271(2): 105613. https://doi.org/10.1016/j.enggeo.2020.105613

[31] 常鸣,窦向阳,唐川, 等, 2019.降雨驱动泥石流危险性评价.地球科学, 44(8):2794-2802.

[32] 崔鹏,郭剑, 2021.沟谷灾害链演化模式与风险防控对策.工程科学与技术,53(3):5-18.

[33] 付建康, 罗刚,胡卸文, 等, 2018. 滑坡堰塞坝越顶溢流破坏的物理模型实验.吉林大学学报(地球科学版), 48(1): 203-212.

[34] 胡卸文, 韩玫,梁敬轩, 等, 2016. 汶川地震灾区泥石流若干关键问题.西南交通大学学报, 51(2): 331-340.

[35] 刘定竺, 崔鹏,蒋德旺, 等, 2017. 堰塞坝溃口展宽过程实验研究.国水土保持科学, 15(6): 19-26.

[36] 柳金峰, 游勇,陈兴长, 等, 2010. 震后堵溃泥石流的特征及防治对策研究——以四川省平武县唐房沟为例.四川大学学报(工程科学版), 42(5): 68-75.

[37] 陆高明, 夏元友, 芮瑞, 2016. 基于弯曲度与链索模型的管涌临界坡降比公式改进. 武汉理工大学学报, 38(3): 41-47.

[38] 单熠博,陈生水,钟启明, 2020.堰塞体稳定性快速评价方法研究.岩石力学与工程学报,39(9):1847-1859.

[39] 石振明,周明俊,彭铭,等, 2021.崩滑型堰塞坝漫顶溃决机制及溃坝洪水研究进展.岩石力学与工程学报, 1-16.

[40] 舒安平, 朱福杨,王澍, 等, 2019. 溃坝泥石流起动过程及其动力学特征. 水利学报, 50(6): 661-669.

[41] 吴梦喜, 高桂云, 杨家修,等, 2019. 砂砾石土的管涌临界渗透坡降预测方法. 岩土力学, 40(3):861-870.

[42] 余斌,杨凌崴,刘清华,等, 2020.基于沟床宽度与颗粒粒径的泥石流精细化预报模型.地球科学,45(4):1447-1456.

[43] 张惠惠, 2011. 震后强震区洪水启动沟道堆积物引发泥石流灾害研究(硕士学位论文).成都:成都理工大学,24-35.

基金资助

国家重点研发计划(2018YFC1505402)

国家自然科学基金项目(42077245;41521002)