非均匀结构堰塞坝溃决过程的3维数值模拟

Three‒Dimensional Numerical Simulation of the Breaching Process of Landslide Dams with Heterogeneous Structures

摘要

堰塞坝是斜坡失稳体快速堆积而成的天然坝体，在斜坡运移分选和河谷约束堆积作用下普遍呈现3维非均匀结构，这导致目前基于均匀结构假设的溃坝模拟误差巨大，严重影响峰值流量预测和风险防控决策。为深入研究非均匀结构对堰塞坝溃坝过程的影响，充分考虑泥沙的不定向输移，建立3维堰塞坝溃坝模型，并采用大涡模拟（LES）与考虑泥沙状态转化的质量守恒方程，准确模拟水沙输运特征与复杂的溃口发育演变。使用该模型分别模拟均匀结构、垂直非均匀结构和水平非均匀结构物理模型试验，计算得到的溃口发展趋势和流量过程线与物理模型试验基本接近。在均匀坝体中，随着材料抗侵蚀能力的增强，溃坝过程中的侵蚀特征由层状冲刷逐渐演变为单级陡坎和多级陡坎。在竖向非均匀结构中，层间材料的组合影响明显，“上细下粗”的顺序结构会导致陡坎侵蚀特征明显，且溃坝流量曲线上没有明显的峰值；相反，“上粗下细”的反序结构在顶部形成“盔甲层”，导致坝体容易在高水位下突然迅速溃决，产生较大的峰值流量。在水平非均匀结构中，过流侧上部是核心区域，其下部和对岸上部分别影响下切速率和扩展速率，对岸下部对溃决影响有限。此外，粗颗粒沉积会引发溃口纵剖坡度变缓，削弱溃口的过流能力，导致溃口的溯源侵蚀受阻，溃决流量峰值较低。

Abstract

Objective Landslide dams are natural blockages formed by the rapid accumulation of slope failure materials, typically resulting from mass movements such as rockslides or debris avalanches. Due to gravitational sorting during downslope transport and complex topographic constraints in narrow valleys, these dams generally exhibit inherently three-dimensional (3D), spatially heterogeneous internal structures. However, most existing breach simulation models assume material homogeneity for computational convenience, neglecting the influence of real-world structural non-uniformity. This simplification results in significant deviations in breach morphology and peak discharge predictions, ultimately impairing the reliability of hazard assessments and early warning systems for downstream communities. This study clarifies the fundamental influence of 3D material heterogeneity on the breach process of landslide dams, enhancing the predictive capabilities of numerical models in geohazard risk ma-nagement. Methods A novel 3D landslide dam breach model was developed by coupling Large Eddy Simulation (LES) with sediment mass conservation equations that incorporated phase transitions between solid and suspended states, accurately capturing the erosion dynamics in heterogeneous dams. The model considered the non-directional transport and deposition behavior of sediment under complex flow conditions. A Volume-of-Fluid (VOF) method was employed to simulate the evolution of the free surface and interface, enabling the precise tracking of water-sediment interactions. The model was validated through physical model experiments simulating three structural types: homogeneous, vertically heterogeneous, and laterally heterogeneous dams. Key breach parameters, including incision rate, breach geometry, and outflow hydrographs, were compared to experimental data to verify the model's reliability and accuracy. Results For homogeneous dams, breach behavior varied significantly with the material grain size. Fine-grained dams exhibited rapid, layered erosion, forming triangular longitudinal profiles, with breach durations under 60 seconds and peak discharges reaching 3.32 L/s. In contrast, coarse-grained dams underwent multistage headcut erosion, requiring higher flow shear stress for sediment entrainment. This delayed breach development resulted in longer times to peak (up to 120 s) and reduced peak discharge (2.1 L/s). Reverse vortices formed at the headcut bases enhanced local scouring but exerted limited influence on overall erosion rates. Medium-grained dams exhibited intermediate characteristics, characterized by single-stage headcuts and corresponding breach metrics that fell between those of fine- and coarse-grained dams. As the median grain size increased, breach morphology transitioned from uniform scouring to progressively complex headcut erosion patterns, with delayed peak times and reduced peak flows. For vertically heterogeneous dams, the breach process was susceptible to the configuration of layered materials. The upper layer (V1) influenced breach initiation: in tests where it consisted of fine particles, breach formation was accelerated; in contrast, coarse-grained V1 layers delayed erosion onset and increased upstream impoundment volumes by up to 40%. The middle layer (V2) governed vertical incision rates: fine-grained layers accelerated downward erosion, whereas coarse layers formed headcuts that restricted further deepening. The lower layer (V3) controlled basal stability: coarse-grained foundations enhanced dam resistance to scouring, while fine-grained ones raised undercutting and subsequent collapse of the overlying mass. In addition, an interactive effect was observed between the various layer combinations. For example, "inverse grading" structures (coarse-over-fine) tended to form armor layers on the surface, delaying breach initiation and leading to sudden failure modes with elevated peak discharges. In contrast, "normal grading" (fine-over-coarse) favored headcut erosion. These structural patterns critically influenced both the erosion mechanisms and the hydraulic response. In laterally heterogeneous dams, spatial variability in material properties across dam zones significantly altered breach morphology and flow dynamics. The breach core zone (C1) determined incision and discharge characteristics, while adjacent zones (C2 and D1) influenced breach asymmetry and propagation direction. The D1 zone played a major role in controlling breach widening, whereas D2 exhibited minimal influence. Comparative simulations (Tests 7~9) revealed differences: Test 7 yielded a wide breach, Test 8 produced a narrow and deep breach, and Test 9 exhibited slow incision but extensive lateral expansion. In addition, the deposition of coarse particles within the breach resulted in a reduction of the slope angle and flow velocity, which suppressed the initiation of further coarse particle entrainment and impeded upstream headward erosion. These feedback mechanisms resulted in a lower peak discharge and a more gradual evolution of the breach. Conclusion This study demonstrates that the spatial heterogeneity of dam materials is a primary factor governing the evolution of breaches in landslide dams. Vertically inverse grading structures (coarse-over-fine) contribute to delayed breach initiation and energy storage, often resulting in abrupt failure events with peak discharges exceeding those of homogeneous dams by more than 30%. In laterally heterogeneous dams, the sorting of materials on the overtopping side directly influences the breach depth-to-width ratio (ranging from 0.8 to 2.5), controlling the discharge capacity. A positive feedback mechanism associated with coarse sediment deposition is identified: the accumulation of sediment reduces the breach slope, which in turn decreases flow velocity and further enhances deposition, transforming the breach mode from rapid incision to a slow-release erosion regime. The developed 3D VOF-LES-based hydro-sediment coupled model overcomes the limitations of the traditional homogeneous assumption and, for the first time, enables high-resolution simulation of breach morphology evolution under realistic heterogeneous conditions. With prediction errors maintained within 10% under complex experimental conditions, this model provides a robust tool for enhancing risk assessments and emergency planning in regions prone to landslide dam breaches.

Graphical abstract

关键词

堰塞坝 / 非均质结构 / 3维数值模拟 / 溃决侵蚀 / 溃决流量

Key words

landslide dam / heterogeneous structure / three‒dimensional numerical simulation / breaching erosion / outflow discharge

本刊网刊

04 ‒ 30

网络出版日期:2024 ‒ 06 ‒ 03

堰塞坝是由崩塌、滑坡等斜坡失稳体快速堆积并堵塞河道而形成的天然坝体^[1‒3]，受失稳地质体类型、材料、运动、河道边界等因素影响，其内部材料分布普遍呈现明显的3维空间非均匀性特征，并极大影响溃坝过程和溃坝程度。以唐家山堰塞坝为例，其坝体沿深度方向主要分为3层，水平方向也呈现明显的颗粒分选现象^[4]。唐家山堰塞坝剖面图如图1所示。通过人工干预，唐家山堰塞坝仅发生部分溃决，峰值流量为6 500 m³/s。深层岩质滑坡造成的小林村堰塞坝，经历了碰撞破碎和长距离颗粒分选，在水平方向形成了非均匀结构^[5]；该堰塞坝以泥沙为主体结构，发生完全溃决，峰值流量高达70 700 m³/s。尽管两者均为基岩滑坡成坝，但溃坝峰值流量差异巨大且难以预测。为更好地评估溃坝风险和进行科学减灾决策，需进一步研究非均匀结构对溃坝的影响机制。

近年来，国内外学者针对堰塞坝的溃决过程与机理开展了大量研究，模型试验是重要技术手段之一。大量学者开展各因素对堰塞坝溃决的影响研究，包括坝体几何形态^[6‒7]、材料组成^[8‒12]和孔隙水压力等外在因素^[13‒17]，但目前的溃决模型试验大多将模型坝体简化为均质坝，对结构特征的关注较少。石振明等^[18]针对非均匀结构对堰塞坝溃决机理的影响开展了相关模型试验，但由于试验的监测数据有限，未能系统揭示非均匀结构对溃坝过程的组合影响。

随着计算机运算和存储技术的蓬勃发展，数值模拟已经被越来越广泛地应用于堰塞坝溃决机理的研究。根据模型计算原理不同，可分为参数模型、简化模型和精细化模型^[19]。参数模型通过回归分析方法拟合历史数据得到经验方程来预测溃决过程的关键参数，包括溃决峰值流量、峰现时间和最终溃口形态^[4,20‒21]。但由于其不能模拟溃坝过程，只能快速且粗糙地对堰塞坝溃决风险进行评估。简化模型对坝体几何形状演变模式进行简化，以模拟滑坡坝的溃决过程，如DABA模型^[22‒23]、DB‒IWHR模型^[24]和DB‒NHRI模型^[25]。虽然简化模型计算效率高，并考虑必要的物理过程，但由于严重依赖预定义的溃坝发展模式，其模拟结果与实际过程仍有较大出入。精细化模型主要通过假定溃口横、纵断面演变模式来逐步计算溃坝洪水过程；基于深度平均理论，采用有限体积法^[26]、有限差分法^[27]或SPH方法^[28‒29]对堰塞坝溃决的真实水土耦合过程进行模拟。

总体而言，此前大多数数值模拟采用1维或2维模型^[30]，3维精细化数值模型已成为当前溃坝模拟研究的热门方向^[31‒32]，但在溃口发展与坝体纵剖面侵蚀特征方面仍缺乏系统研究。同时，当前的3维溃坝模型均假设堰塞坝为均匀结构，未考虑非均匀结构对溃决过程，尤其是溃口形状的影响，故对非均匀堰塞坝的溃决机理仍需要深入研究。

因此，本文基于水沙运动方程与流体体积法（VOF）建立3维精细数值模型，对均匀结构、竖向非均匀与水平非均匀结构开展9个材料空间分布不同的工况模拟，并用水槽试验结果对各工况的侵蚀特征与溃决流量进行对比验证。同时，对竖向与水平两种非均匀结构各增加3个补充工况，系统讨论坝体内不同材料分配组合对溃决的影响，为堰塞体溃决的风险评估提供科学参考。

1 数值方法

本文数值模型由水动力模型和泥沙输运模型两个模块组成。水动力模型采用雷诺平均Naiver‒Stokes（RANS）方程和大涡模拟（LES）计算溃决过程中复杂的水力特性。泥沙输运模型考虑了泥沙颗粒在推移质与悬移质两种状态间的转化，使用质量守恒方程计算溃口形状。此外，通过流体体积法（VOF）实现相邻网格单元的泥沙输运，从而实现水沙液面的定位。

1.1 水动力模型

由于溃坝水流问题具有较强的湍流特征，传统的雷诺平均Navier‒Stokes方程（RANS）在描述其湍流特性时存在着很大的不确定性，而直接数值模拟（DNS）又需要极高的计算资源，因此采用LES方法计算湍流特征。LES采用滤波函数，将涡旋结构中比滤波函数尺度小的涡旋过滤掉，以进一步分解出描写大涡流的方程^[33]：

\frac{\partial }{\partial t}\left(\rho {\overline{u}}_i\right)+\frac{\partial }{\partial x_i}\left(\rho \overline{u}{\overline{u}}_j\right)=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial {\overline{u}}_i}{\partial x_j}\right)-\frac{\partial {\tau }_{ij}}{\partial x_j}（1）

\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}\left(\rho {\overline{u}}_i\right)=\mathrm{0}（2）

式（1）～（2）中：\rho为流体密度；x_i 、x_j 为两个不同方向的位置坐标；t为时间；{\overline{u}}_i、{\overline{u}}_j为不同方向的流速，i、j表示x或y或z方向；μ为湍流黏性系数；{\tau }_{ij}为亚格子尺度应力，它体现了被滤掉的小涡对大涡的影响，可表示为：

{\tau }_{ij}=\rho {\overline{u_{i}u}}_j-{\overline{u_i\overline{u}}}_j（3）

由LES计算得到的结果具有高精度、较为真实的流动结构和湍流特征，适用于解决溃坝水流等需要获得较为精确流动信息的问题。自由液面信息则采用VOF法获得。

1.2 泥沙输运模型

堰塞体溃决过程中，泥沙颗粒可以分为推移质和悬移质两种状态，包含4种运动方式：悬移质运动、沉降运动、推移质运动及冲刷运动（图2（a））。因此，溃口形态的发展采用包含相关运动方式的质量守恒公式进行模拟：

\varphi \frac{\partial z}{\partial t}=\left(\frac{\partial q_{\mathrm{b}x}}{\partial x}+\frac{\partial q_{\mathrm{b}y}}{\partial y}+E+D\right)（4）

式中，z为底床高程，q_bx 和q_by 分别为x和y方向单宽体积推移质输移率，\varphi为最大堆积分数，D为沉降通量，E为冲刷通量。

由于单个颗粒周围的水动力和界面处的边界层很难计算，因此采用经验模型进行坝体的冲蚀模拟。根据坝体材料的抗冲刷性能（临界希尔兹数）判断泥沙颗粒是否达到冲刷启动条件。泥沙启动后可进行冲刷运动，使用经验公式计算出相应水流条件下的推移质向悬移质的转化通量E（即冲刷通量）^[34]，公式如下：

E={\alpha }_{\mathrm{i}}{\mathit{n}}_{\mathrm{s}}{d_{\mathrm{*}}}^{\mathrm{0.3}}{\left({\theta }_{\mathrm{i}}-{\theta }_{\mathrm{c}\mathrm{i}}\right)}^{\mathrm{1.5}}\sqrt[]{\frac{‖g‖d_{\mathrm{i}}\left({\rho }_{\mathrm{i}}-{\rho }_{\mathrm{w}}\right)}{{\rho }_{\mathrm{w}}}}（5）

式中，{\alpha }_{\mathrm{i}}为坝体材料的挟带率系数， n_s为充填溃口底床的法向矢量，d_{\mathrm{*}}为无量纲坝料粒径参数，d_{\mathrm{i}}为不同材料的粒径，{\rho }_{\mathrm{i}}为坝料密度，{\rho }_{\mathrm{w}}为水密度，θ_i为临界希尔兹数，θ_ci为无量纲临界希尔兹数。同时，泥沙颗粒也可以进行推移质运动，结合泥沙粒径与水流条件，选取Meyer‒Peter和Muller的经验公式表示推移质单宽输沙^[35]：

q_{\mathrm{b}}=K{\left({\theta }_{\mathrm{i}}-{\theta }_{\mathrm{c}\mathrm{i}}\right)}^{\mathrm{1.5}}{\left[g\frac{\left({\rho }_{\mathrm{i}}-{\rho }_{\mathrm{w}}\right)}{{\rho }_{\mathrm{w}}}{d}_{\mathrm{50}}^{\mathrm{3}}\right]}^{\mathrm{0.5}}（6）

式中，K为推移质系数，d₅₀为坝料平均粒径。在推移质的输移过程中，悬移质颗粒因重量而发生沉降运动，沉降通量D可使用Soulsby^[33]提出的沉降速度进行计算：

D=\frac{v_{\mathrm{f}}}{d_{\mathrm{i}}}\left[{\left(\mathrm{10}.{\mathrm{36}}^{\mathrm{2}}+\mathrm{1.049}{d}_{\mathrm{*}}^{\mathrm{3}}\right)}^{\mathrm{0.5}}-\mathrm{10.36}\right]c_{\mathrm{i}}（7）

式中，v_f为流体的运动黏度，c_i为不同坝料的悬移质浓度。改变后的溃口形态将影响水流的运动，导致水力特性发生相应变化，进而进入下一步冲刷和沉积判断。由于悬移质随水运动且不影响坝体形态，故无需单独求解。

采用基于欧拉模型的VOF法描述复杂的3维溃口形态变化。在使用VOF法计算水沙交界面时，每个网格都被赋予一个值表示泥沙的充满程度，用于确定界面的位置和形状（图2（b））。网格里泥沙的最大充满程度为65%，即当网格里的泥沙体积达到65%时，该网格就被判定为淤积泥沙。网格单元中的泥沙体积减少时，意味着该网格存在泥沙向悬移质转化运动或向相邻网格单元的输移运动。相邻单元间的泥沙输移量与主要输移方向均由网格处的水力条件，尤其是沙床剪切应力决定。泥沙在网格间的不定向移动意味着不仅可以计算坝体材料沿主流方向的下切侵蚀，也可以模拟溃口侧坡处发生的横向冲刷（图2（c）），进而准确描述溃口的3维形态。此外，数值模型不考虑渗流作用。上述3维溃坝模型基于FLOW‒3D软件平台耦合实现：在每个时间步中计算河床高度，得到局部希尔兹数和相关参数，并根据溃口形状公式计算河床高程变化。根据新的地形通过求解水动力方程，更新流场的流体体积分数和水力要素，实现水沙两相模型耦合计算。

1.3 模型设置

坝体模型的几何尺寸如图3（a）所示，坝体纵剖面（顺河向）为梯形，坝高为24 cm，坝顶宽为24 cm，坝底宽为108 cm，上游坝坡坡比为1.0∶2.0（26.6°），下游坝坡坡比为1.0∶1.5（33.7°）。坝体横剖面（横河向）为矩形，坝长为40 cm，等于水槽宽度。根据堰塞坝成坝机制和抢险实例，在靠近水槽边壁坝顶一侧开挖泄流槽，其横截面为三角形，截面积为20 cm²（高5 cm、宽8 cm）。在泄流槽与上游坡面交界处设置监测面1，在坝体下游坡脚处设置监测面2（图3（b）），分别用来监测溃口断面形态与溃决流量。

根据堰塞坝的堵江成坝机制，首先选取均匀、竖向非均匀、水平非均匀3种坝体结构类型共9个工况，见表1～3。竖向非均匀结构的内部材料分布情况如图3（c）所示，坝体沿深度方向分为3层，分别命名为V1、V2和V3；水平非均匀结构按照坝体横断面上的4个象限划分将坝体分为4个区域，分别命名为C1、C2、D1和D2，并对坝体处的网格进行加密处理（图3（d））。其中，C1和C2为坝体过流侧（即开挖泄流槽一侧），D1和D2为坝体对岸侧，每个区域的材料颗粒级配不同。根据2008年汶川地震诱发形成的东河口、唐家山和小岗剑堰塞坝的现场土样数据^[14,22]，选取3种典型的堰塞坝体材料级配类型，级配曲线如图4所示。不同材料的d₅₀分别为1.6、3.6和6.0 mm，并根据粒径相对大小，依次定义为细粒材料（F）、中粒材料（M）和粗粒材料（G），模型输入参数见表4。

2 计算结果

为便于后续的定量分析，根据侵蚀起点的位置将溃决过程分为两个阶段，如图5所示。溃口形成阶段（阶段Ⅰ）：从泄流槽全线过流开始，到侵蚀点移动至上游坡面时结束。溃口发展阶段（阶段Ⅱ）：从侵蚀发展至上游坡面开始（溃决流量迅速增加并很快达到峰值）到出入流平衡结束（侵蚀起点沿上游坝坡向下发展至稳定）。

2.1 均匀坝溃决过程

工况1在溃决过程中主要发生层状冲刷，纵剖面逐渐发展成三角形，溃决历时较短，出现较大流量峰值。为方便观测，图6从3个视角展示工况1的溃决过程。由图6可知，溃决开始后，下游坡面首先发生冲刷，溃口向上游移动的过程中纵坡率基本不变。当坝体纵剖面由梯形变为近似三角形时，侵蚀起点发育至上游坝坡。溃决水流在通过斜坡后流速提高，在下游坡脚处形成纵向高流速区域，加速了溃口的侵蚀速率。同时，溃口上游横向流速增加，开始侧向侵蚀。随着溃口快速下切和展宽，溃决流量迅速增加至峰值3.32 L/s（图7），仅低于实测值5.6%，达峰时间较试验结果提前6.7%。峰值过后，溃决流量逐渐减小至与入流平衡，坝前水位缓慢下降，溃口的纵坡率也逐渐减小，最终发展至上下游宽、中间窄的喇叭状，直至溃决结束。

图8为工况2、3的溃决过程。由图8可知：随着坝体材料中值粒径增大，纵剖面侵蚀特征由层状侵蚀（工况1）转变为陡坎侵蚀（工况2），进一步发展为多级陡坎侵蚀（工况3）；同时峰值流量减小，峰现时间推迟。工况2中粒材料启动需更大水流剪切力，表现出更强的抗侵蚀能力；溃决初期水流流速较慢，下切侵蚀能力较弱，工况2发生单级陡坎侵蚀。工况3坝体材料换为最难启动的粗粒材料，下游坝坡更稳定，初期只能形成高度较小的多级陡坎；流量峰值较工况2降低是由于漫顶水流流经陡坎后具备较大势能，产生反向旋流，加速陡坎底部掏蚀。陡坎高度越小，陡坎侵蚀越慢，溃决历时越长，峰值流量也越小。3个工况的溃决曲线与模型试验结果趋势一致，拟合较好（图7）。

2.2 竖向非均匀坝溃决过程

工况4的冲蚀特征受不同土层组合的影响，如图9（a）所示，上层为层状侵蚀，中层出现单级陡坎，下层基本稳定。

由于泄流槽的存在，V1层细粒材料的厚度较小，与中层交界处很快出现坡度变化点，并且在V2层发展为单级陡坎；V3层的粗粒材料抗侵蚀能力比较强，溃决水流无法破坏坝底，该层在溃决过程中基本保持稳定。陡坎高度受限于V2层的厚度，陡坎溯源速度比较慢，溃决流量缓慢增加，并且不出现峰值，该结果与试验实测值吻合。

工况5的溃口坡度基本不变（图9（a）），V3层冲刷明显，坝体溃决迅速，峰值流量较小。

溃决开始以后，工况5上层的细粒材料很快启动，并且被溃决水流带走，溃口很快下切发展至中下层。而溃决水流越过V2层后掏蚀坝底的中粒材料，下层的不稳定加速了上部坝体的侵蚀失稳，导致溃口下切侵蚀速率快。流量峰值于58 s达到峰值，为2.87 L/s，将数值结果与试验对比可以发现，工况5的峰值流量较实测值大5.9%，达峰时间较实测提前了9.2%，如图9（b）所示。

工况6的V1层侵蚀缓慢，纵剖面在阶段Ⅰ出现坡度变化，溃决流量相对最大（图9（a））。

由于V1层为抗侵蚀能力较强的中粒材料，溃口的下切侵蚀速率较低，延缓了溃坝过程。同时，泄流槽过流能力较弱，导致坝前水位不断雍高。当溃口缓慢发展至上游坝坡时，大量蓄水短时间内下泄，工况6的溃决流量于92 s达到峰值4.41 L/s，明显大于工况5。工况6流量峰值比实测值小了6.8%，达峰时间比模型试验提前了9.5%（图9（b））。两个关键溃决参数的相对误差均在10%之内，表明数值模拟结果与模型试验较吻合。

2.3 水平非均匀坝溃决过程

工况7在阶段Ⅰ有明显的陡坎侵蚀特征，阶段Ⅱ后溃口侧坡开始横向展宽，最终溃口宽度较大（图10（a））。C1与C2交界处出现高度较大的单级陡坎，加速了溃口的溯源侵蚀，75 s时达峰值2.80 L/s，较实测值大8.5%，且达峰时间较模型试验早11.6%。阶段Ⅰ中，溃口展宽不明显，但当溃口发展至上游坝坡后，部分水流横向汇入泄流槽，对溃口的侧坡产生冲蚀。展宽至D1区的中粒材料后，展宽速率降低，溃口趋于稳定。

工况8在溃决过程中以下切侵蚀为主，溃口基本不产生横向发育，导致其溃口更深且窄（图10（b））。过流侧为细粒材料，主要发生层状冲刷，溃口下切速率较大，导致溃决流量快速达到峰值3.73 L/s。如图11（a）所示，该工况与试验吻合度最高，流量峰值相对误差为2.2%，达峰时间相对误差为6%。在阶段Ⅱ中，上游溃口处横向流速较低，侧坡侵蚀程度小，最终溃口比工况7更深、更窄。

工况9在溃决过程中的剖面纵坡度基本一致，溃口下切缓慢，展宽程度最大（图10（c））。由于中粒材料的阻碍作用，C1区冲刷速度较慢，从而导致达峰时间晚，峰值流量仅为2.83 L/s。同时，由于溃口下切受阻，大量水流只能横向流动。在上游溃口处出现大范围的横向高流速区域，导致溃口大幅横向展宽，很快发育至D1区，最终溃口形状浅且宽（图11（b））。

3 非均匀结构响应分析

3.1 竖向非均匀结构对溃口下切的影响

为深入探究非均匀结构对溃决过程的影响，针对竖向非均匀坝和水平非均匀坝分别增加了3个补充工况。竖向非均匀坝补充工况见表5，各工况可根据初期的溃决特征分为无陡坎、单级陡坎与多级陡坎3类，如图12所示。综合侵蚀特征与溃决流量，讨论各层对溃坝的不同作用及组合影响。

竖向非均匀结构中，不同结构层对溃决特征，尤其对溃口下切发展模式的影响不同。V1层主要影响溃口的形成阶段。如图13所示，工况4与5的V1层为细粒材料，导致溃决开始后溃口迅速形成，并很快冲刷至中层。V2层影响溃口下切速率，进而决定了阶段Ⅰ的历时。工况11的中层细粒材料迅速掏蚀，溃口很快发展至上游；工况12的V2层为中粒材料，导致溃口的下切发育速率较慢。V3层决定了坝底稳定性，影响纵剖面的溃决特征。工况12底部为细粒材料，对中上部缺乏足够的支撑，导致V3层掏蚀严重，溃口以层状冲刷方式溯源发展。

不同层间的材料组合对溃决的影响更加明显，其中，上中层的材料组合主要影响溃决历时与溃决流量。当中上层为“上细下粗”的顺序结构时，溃决时间较短，溃决水流快速增加至峰值或稳定（工况4、5、10）。当出现“上粗下细”的反序结构时，受上层的盔甲层的影响，可能会引起较大的流量峰值。例如：工况6的上层为中粒材料，溃决初期缓慢的溃决水流较难使其启动，导致陡坎发育速率缓慢，其下软弱层未发生冲蚀；在阻碍侵蚀进程、达峰时间延后的同时，泄流能力受阻，坝前水位不断提高，坝前雍水量增加。所以当溃口发展至上游坝坡时，坝体突然溃决，流量峰值最大。

中下层的材料组合决定溃决过程中的侵蚀特征。当这两层的材料表现为顺序结构时，界面交界处的坡度产生变化，在V2层处出现陡坎侵蚀特征（工况4、11）。溃决水流流经中下层交界面时受阻增大，流速减小，侵蚀能力也随之减弱。同时，由于在V2层发生侵蚀，其纵坡度会迅速增大直至形成陡坎，陡坎的存在也使得V2层的水流流速进一步增大，形成明显的高流速区，进而促进陡坎的溯源发育。当中下层为反序结构时，溃口纵剖面未出现明显坡度变化（工况5、12）。这是由反序结构导致溃决水流更易掏蚀坝底，而下层的不稳定加速了上部坝体的侵蚀失稳，因此溃决过程以层状冲刷为主；同时，溃口溯源速率较快，坝体容易突然出现溃决。

结构层的组合影响坝体的稳定性。工况4的上中下3层的材料粒径逐级增加，为顺序结构，此结构的坝体稳定性高。该工况的溃口发展模式为单级陡坎侵蚀，在陡坎向前发育中，由于溃决水流持续下泄，流量缓慢增至入流量后保持稳定，无峰值出现（图（13））。工况12则为反序结构，细粒材料的坝底加剧了坝体的不稳定。同时，中上层共同形成V3层的厚盔甲层，导致该工况中流量峰值出现较晚。

3.2 水平非均匀结构对溃口展宽的影响

水平非均匀结构的补充工况见表6。

C1区是水平非均匀坝溃决的核心区域，下切侵蚀特征主要由过流侧材料性质差异决定。水平非均匀坝最终溃口横断面对比如图14所示。工况7过流侧材料分布与工况4类似，因此也出现相似的单级陡坎；但由于工况7底部缺乏稳定层，水流充分向下冲蚀，导致其陡坎高度明显大于工况4，溃决方式也由工况4的均匀泄流转变为突然溃决。根据C1的不同材料特性，水平非均匀工况可分为粗粒和细粒两类。其中，细粒工况最终溃口形态通常较宽（图14（a）），这是由于细粒材料的休止角较小，侧坡侵蚀较快，溃口很快展宽至D1区域。

C2区影响溃口下切速率，并与D1区共同影响溃口发展的优势方向。以工况8为例，C2区材料的抗侵蚀能力明显低于D1区的中粒材料，所以溃口以下切侵蚀为主，最终溃口相对工况7深且窄。

D1区主要影响溃口的展宽（图14），工况7与13的C2区虽均为中粒材料，但后者D1区换为粗粒后，最终溃口明显变窄。这种溃口展宽的差异也影响了溃决流量，工况13流量峰值稍低于工况7，如图15所示，原因是溃口的缩窄导致过流能力减弱。工况9与15溃口形状基本一致（图14），这是由于粗颗粒的阻碍作用，溃口未侧向发育至D1区就达到平衡阶段，并不再变化。

D2区对溃决的影响有限。所有工况溃口发育基本都没有发育到D2区，其材料性质对溃坝过程的影响较小。

受材料性质、渗流和饱和度等因素影响，溃口侧坡的发展规律非常复杂，其数值模拟具有较大难度。采用休止角计算溃口形态虽然能够反映材料分区的影响，但由于未考虑土颗粒的黏性特性，无法很好地模拟溃坝过程中常见的侧坡失稳崩塌现象。因此，为更精确地模拟溃口发展，尤其是其展宽模式，考虑在数值模型中引入土颗粒黏性特性，这将有助于更准确地模拟溃口的演变过程。

3.3 颗粒沉积对溃决流量的影响

粗颗粒沉积会减缓溃口纵剖面坡度，溃口泄流能力降低，在加剧沉积现象的同时阻碍溃口溯源侵蚀过程，导致溃决流量峰值降低。图16为工况11的粗颗粒沉积分布。由图16可知：溃决初期（t=50 s），溃口发生溯源侵蚀，坝底基本稳定，纵坡接近直线且坡度较大。上层粗粒材料启动后大多被水流挟带流走，仅有小部分沉积在溃口下游，此特征与纵坡坡度较大相关。

斜坡上的粗颗粒受力情况如图17所示，颗粒主要受重力F_w、支撑力F_l、水流剪切力F_d与阻力F_r作用。其中，F_d为溃决水流对泥沙颗粒产生的前驱力，方向与水流方向相同。由于数值模型不考虑渗流过程，所以忽略渗流力的影响。

因此，斜坡上粗颗粒的平衡状态为：

F_{\mathrm{d}}+F_{\mathrm{w}}\mathrm{s}\mathrm{i}\mathrm{n}\beta =F_{\mathrm{w}}\mathrm{c}\mathrm{o}\mathrm{s}\beta（8）

式中，\beta为斜坡坡度。

由式（8）可知，斜坡坡度越大，粗颗粒自重的前驱分量越大，导致其启动所需的水流剪切力越小，沉积颗粒更容易冲刷输运。同时纵坡坡率一致的溃口过流能力较大，溃决流量迅速增加，加速溃口的下切侵蚀。

当t=75 s时，溃口逐渐发展为坡度沿程逐渐减小的变坡度纵坡，同时，伴随着溃口侧坡的不断失稳，溃口持续横向展宽，大量上层粗颗粒进入溃决水流，并在溃口纵坡中层形成局部沉积。粗颗粒的沉积则加剧纵坡度变缓，导致水流流速降低，溃口的过流能力减弱，溃决流量增长梯度降低。随着粗颗粒的缓慢输移，当t=105 s时，纵坡坡脚处形成大范围粗颗粒沉积，纵坡度受沉积影响逐渐趋于一致，导致溃决流量出现短暂波动现象，并逐渐稳定。

4 结 论

基于3维水动力模型和泥沙输运模型，开展一系列不同结构的非均匀结构坝溃坝数值模拟，以研究非均匀结构对溃决过程及溃决流量的影响，结论如下。

1）材料抗侵蚀能力决定了溃决初期的断面侵蚀特征与溃决流量。抗侵蚀能力的增加导致溃决初期的纵剖面形状由坡度不变的层状冲刷发展到单级陡坎，最终可能出现多级陡坎侵蚀现象，溃决流量峰值也随之降低，达峰时间逐渐延后。

2）竖向非均匀坝的溃决特征受不同土层组合的控制。上层影响溃口的形成，中层影响溃口下切速率，下层影响坝体稳定性。上、中层的材料组合主要影响溃决历时与溃决流量，其中，“上粗下细”的反序结构容易延缓溃决过程，引起突然溃坝，导致流量出现较大峰值。中、下层间材料组合决定了侵蚀特征，“上细下粗”的顺序结构会引发陡坎侵蚀特征，而反序结构的非均匀溃坝则以层状冲刷为主。

3）水平非均匀坝可反映粗颗粒分选对溃坝的影响。过流侧上部是溃坝的核心区，主要影响溃决特征与溃决流量。过流侧下方与对岸侧上方分别影响溃口下切和展宽，决定溃口形状。这两个区域共同决定了溃决优先发育方向。过流侧下部区域对溃决影响不大。

4）粗颗粒沉积与溃决流速降低形成恶性循环，对溃口溯源侵蚀的阻碍作用放大。大量粗颗粒沉积导致溃口纵剖面坡度变缓，溃决水流流速降低，进一步阻碍粗粒材料的冲刷启动。同时，缓坡溃口的过流能力较差，溃口的溯源侵蚀过程受阻，溃决流量峰值较低。

后续将对当前3维非均质溃坝模型持续优化，并考虑土颗粒黏性与渗流作用，实现统一平台下考虑侧坡失稳的溢流冲刷‒渗流潜蚀耦合堰塞坝溃决模拟。

参考文献

[1] Casagli N, Ermini L, Rosati G.Determining grain size distribution of the material composing landslide dams in the Northern Apennines:Sampling and processing methods[J].Engineering Geology,2003,69(1/2):83‒97. doi:10.1016/s0013-7952(02)00249-1

[2] Costa J E, Schuster R L.The formation and failure of natural dams[J].Geological Society of America Bulletin,1988,100(7):1054‒1068. doi:10.1130/0016-7606(1988)100<1054:tfafon>2.3.co;2

[3] Fan Xuanmei, Dufresne A, Whiteley J,et al.Recent technological and methodological advances for the investigation of landslide dams[J].Earth‒Science Reviews,2021,218:103646. doi:10.1016/j.earscirev.2021.103646

[4] Peng M, Zhang L M.Analysis of human risks due to dam-break floods:Part 1:A new model based on Bayesian networks[J].Natural Hazards,2012,64(1):903‒933. doi:10.1007/s11069-012-0275-5

[5] Zhou Yuanyuan, Shi Zhenming, Qiu Tong,et al.Experimental study on morphological characteristics of landslide dams in different shaped valleys[J].Geomorphology,2022,400:108081. doi:10.1007/s10346-022-01907-6

[6] Xu Fugang, Yang Xingguo, Zhou Jiawen,et al.Experimental research on the dam‒break mechanisms of the Jiadanwan landslide dam triggered by the Wenchuan earthquake in China[J].The Scientific World Journal,2013,2013:272363. doi:10.1155/2013/272363

[7] Huang Wei, Qi Zijie, Duan Wengang,et al.Large‒scale experiment on impact of downstream slope on breaching mechanisms of landslide dams[J].Advanced Engineering Sciences,2022,54(3):14‒24.

[8] 黄卫,齐子杰,段文刚,等.堰塞坝背水面坡度对溃决过程影响机理大尺度试验[J].工程科学与技术,2022,54(3):14‒24.

[9] Cao Zhixian, Yue Zhiyuan, Pender G.Landslide dam failure and flood hydraulics.Part I:Experimental investigation[J].Natural Hazards,2011,59(2):1003‒1019. doi:10.1007/s11069-011-9814-8

[10] Yang Yang, Cao Shuyou, Yang Kejun,et al.Experimental study of breach process of landslide dams by overtopping and its initiation mechanisms[J].Journal of Hydrodynamics,Ser B,2015,27(6):872‒883. doi:10.1016/s1001-6058(15)60550-9

[11] Chen S C, Lin T W, Chen C Y.Modeling of natural dam failure modes and downstream riverbed morphological changes with different dam materials in a flume test[J].Engineering Geology,2015,188:148‒158. doi:10.1016/j.enggeo.2015.01.016

[12] Jiang Xiangang, Wei Yunwei, Wu Lei,et al.Experimental investigation of failure modes and breaching characteristics of natural dams[J].Geomatics,Natural Hazards and Risk,2018,9(1):33‒48. doi:10.1080/19475705.2017.1407367

[13] Jiang Xiangang, Huang Jiahua, Wei Yunwei,et al.The influence of materials on the breaching process of natural dams[J].Landslides,2018,15(2):243‒255. doi:10.1007/s10346-017-0877-9

[14] Shi Zhenming, Wang Youquan, Peng Ming,et al.Landslide dam deformation analysis under aftershocks using large-scale shaking table tests measured by videogrammetric technique[J].Engineering Geology,2015,186:68‒78. doi:10.1016/j.enggeo.2014.09.008

[15] Shi Zhenming, Wang Youquan, Peng Ming,et al.Characteristics of the landslide dams induced by the 2008 Wenchuan earthquake and dynamic behavior analysis using large‒scale shaking table tests[J].Engineering Geology,2015,194:25‒37. doi:10.1007/978-3-319-09057-3_202

[16] Itoh T, Ikeda A, Nagayama T,et al.Hydraulic model tests for propagation of flow and sediment in floods due to breaking of a natural landslide dam during a mountainous torrent[J].International Journal of Sediment Research,2018,33(2):107‒116. doi:10.1016/j.ijsrc.2017.10.001

[17] Peng Ming, Jiang Qiulu, Zhang Qingzhao,et al.Stability an-alysis of landslide dams under surge action based on large-scale flume experiments[J].Engineering Geology,2019,259:105191. doi:10.1016/j.enggeo.2019.105191

[18] Peng Ming, Ma Chenyi, Chen Hongxin,et al.Experimental study on breaching mechanisms of landslide dams composed of different materials under surge waves[J].Engineering Geology,2021,291:106242. doi:10.1016/j.enggeo.2021.106242

[19] Shi Zhenming, Zhang Gongding, Peng Ming,et al.Experimental investigation on the breaching mechanisms of landslide dams with heterogeneous structures[J].Advanced Engineering Sciences,2023,55(1):129‒140.

[20] 石振明,张公鼎,彭铭,等.非均质结构堰塞坝溃决机理模型试验[J].工程科学与技术,2023,55(1):129‒140.

[21] Wu Weiming.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

[22] Shan Yibo, Chen Shengshui, Zhong Qiming,et al.Development of an empirical model for predicting peak breach flow of landslide dams considering material composition[J].Landslides,2022,19(6):1491‒1518. doi:10.1007/s10346-022-01863-1

[23] Shen Danyi, Shi Zhenming, Peng Ming,et al.Longevity analysis of landslide dams[J].Landslides,2020,17(8):1797‒1821. doi:10.1007/s10346-020-01386-7

[24] Peng M, Zhang L M, Chang D S,et al.Engineering risk mitigation measures for the landslide dams induced by the 2008 Wenchuan earthquake[J].Engineering Geology,2014,180:68‒84. doi:10.1016/j.enggeo.2014.03.016

[25] Chen Chen, Zhang Limin, Xiao Te,et al.Barrier lake bursting and flood routing in the Yarlung Tsangpo Grand Canyon in October 2018[J].Journal of Hydrology,2020,583:124603. doi:10.1016/j.jhydrol.2020.124603

[26] Cai Yaojun, Cheng Haiyun, Wu Shuaifeng,et al.Breaches of the baige barrier lake:Emergency response and dam breach flood[J].Science China Technological Sciences,2020,63(7):1164‒1176. doi:10.1007/s11431-019-1475-y

[27] Mei Shengyao, Chen Shengshui, Zhong Qiming,et al.Effects of grain size distribution on landslide dam breaching—Insights from recent cases in China[J].Frontiers in Earth Science,2021,9:658578. doi:10.3389/feart.2021.658578

[28] 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

[29] Zhou G G D, Lu Xueqiang, Xie Yunxu,et al.Mechanisms of the non‒uniform breach morphology evolution of landslide dams composed of unconsolidated sediments during overtopping failure[J].Journal of Geophysical Research (Earth Surface),2022,127(10):e2022JF006664. doi:10.1029/2022jf006664

[30] Bui H H, Nguyen G D.Smoothed particle hydrodynamics (SPH) and its applications in geomechanics:From solid fracture to granular behaviour and multiphase flows in porous media[J].Computers and Geotechnics,2021,138:104315. doi:10.1016/j.compgeo.2021.104315

[31] Feng Ruofeng, Fourtakas G, Rogers B D,et al.Two‒phase fully-coupled smoothed particle hydrodynamics(SPH) model for unsaturated soils and its application to rainfall-induced slope collapse[J].Computers and Geotechnics,2022,151:104964. doi:10.1016/j.compgeo.2022.104964

[32] Issakhov A, Borsikbayeva A, Issakhov A.Dam-break flow on mobile bed through an idealized city:Numerical study[J].Water Resources Management,2022,36(11):4425‒4446. doi:10.1007/s11269-022-03253-7

[33] Mei Shengyao, Chen Shengshui, Zhong Qiming,et al.Detailed numerical modeling for breach hydrograph and morphology evolution during landslide dam breaching[J].Landslides,2022,19(12):2925‒2949. doi:10.1007/s10346-022-01952-1

[34] Zhang Lucheng, Zhong Qiming, Yang Meng,et al.Centrifugal model tests and numerical modeling on overtopping-induced breach processes of landslide dams[J].Frontiers in Earth Science,2023,10:1062981. doi:10.3389/feart.2022.1062981

[35] Soulsby R.Ch. 9:Bedload transport[M].London:Thomas Telford Publications,1997. doi:10.1016/j.coastaleng.2005.04.003

[36] Mastbergen D R, Van Den Berg J H.Breaching in fine sands and the generation of sustained turbidity currents in submarine canyons[J].Sedimentology,2003,50(4):625‒637. doi:10.1046/j.1365-3091.2003.00554.x

[37] Meyer‒Peter E, Müller R.Formulas for bed‒load transport[C]//Proceedings of the 2nd Meeting of the International Association for Hydraulic Structures Research,1948.

基金资助

国家自然科学基金项目(42061160480)

国家自然科学基金项目(42071010)

国家自然科学基金项目(U23A2044)