黄河上游德恒隆古滑坡-堰塞湖-溃决洪水灾害链全过程模拟与灾害放大效应

Runout Simulation and Disaster Amplification of Dehenglong Paleolandslide-Dammed Lake-Outburst Flood Chain in the Upper Yellow River

摘要

位于青藏高原东北缘的黄河上游峡谷段发育了大量巨型古滑坡-堰塞湖-溃决洪水灾害链，厘清其运动演化过程对防控此类灾害链具有重要意义.在野外调查基础上，通过构建深度平均滑坡-堵河运动学模型与考虑颗粒物运移及物理侵蚀沉积机制的溃决洪水动力学耦合模型，以其中典型的德恒隆巨型滑坡灾害链为例开展模拟研究，系统分析了滑坡堵河-堰塞湖蓄水-溃决洪水的完整运动学过程.结果表明：德恒隆滑坡体积达35亿m³，峰值平均速度约44.5 m/s，总滑移时长约120 s，形成堰塞坝高达234 m，堵塞黄河形成面积约280 km²、库容约232.3亿m³的巨型堰塞湖；堰塞坝溃决过程持续约100 h，洪峰流量达23.8万m³/s，洪水沿黄河河道波及下游900 km流域，最高洪水位达158 m.相比于滑坡本身的运动距离，本次巨型滑坡-堵河-溃坝洪水灾害链波及了上下游约1 000 km河道范围内的广大区域，具有显著的灾害放大效应.研究成果可为黄河流域滑坡灾害链的风险评估与防灾决策提供科学依据.

Abstract

The upper Yellow River canyon on the northeastern margin of the Tibetan Plateau has developed numerous giant paleo-landslide-dammed lake-outburst flood disaster chains. Understanding their kinematic evolution is crucial for mitigating such cascading hazards. Based on field investigations, this study constructs a depth-averaged landslide damming kinematic model and a coupled outburst flood dynamics model incorporating sediment transport and physical erosion/deposition mechanisms, using the typical Dehenglong giant landslide disaster chain as a case study. The simulation systematically analyzes the entire kinematic process, including landslide damming, lake formation, and dam breaching. The results indicate that the Dehenglong landslide had a volume of 3.50 billion m³, reached a peak average velocity of 44.5 m/s, and had a total sliding duration of 120 s. This formed a 234 m-high dam that blocked the Yellow River, creating a giant dammed lake with an area of 280 km² and a storage capacity of 23.23 billion m³. The dam breach lasted 100 h, with a peak discharge of 238 300 m³/s. The resulting flood propagated 900 km downstream along the Yellow River, reaching a maximum water level of 158 m. Compared to the landslide movement distance itself, this giant landslide-damming-outburst flood disaster chain affected a vast area spanning approximately 1 000 km of the river channel, exhibiting a significant disaster amplification effect. The findings of this study provide a scientific basis for risk assessment and disaster prevention decision-making regarding landslide disaster chains along the Yellow River Basin.

Graphical abstract

关键词

深度平均法 / 巨型滑坡 / 灾害链 / 全过程模拟 / 工程地质学.

Key words

depth-averaged model / giant landslide / disaster chain / full-process simulation / engineering geology

0 引言

位于青藏高原东北缘的黄河上游峡谷段构造活动强烈、地貌演化迅速，该区域内重大地质灾害类型多样、分布广泛且突发性强，灾害链效应显著，严重威胁下游城市安全（彭建兵等， 2004，2025； 兰恒星等， 2022）.近年黄河科考研究表明，黄河流域上游峡谷段发育多处巨型古滑坡，部分滑坡处于局部蠕变状态，在极端气候或地震扰动下存在复活堵河风险，对区域防洪安全构成严重威胁（金钊等， 2024）.特别是在黄河上游龙羊峡至刘家峡段，由于特殊的地形地貌、地层岩性、地质作用以及古气候条件，历史上形成了多达30余个巨型滑坡（李小林等， 2011）.其中，德恒隆滑坡位于尖扎盆地出口，体积达数十亿立方米，在滑坡滑动过程中，其前缘刨蚀河床，最终形成高度逾200 m的堰塞坝体（郭小花， 2018），完全堵塞黄河干流，导致上游形成长度超过45 km的堰塞湖（殷志强等， 2021），最终因坝体溃决引发洪水，构成了“滑坡堵河-堰塞湖蓄水-溃决洪水”的完整灾害链.当前对黄河上游历史巨型滑坡灾害链运动演化过程的研究较为缺乏，制约了潜在灾害链的预测精度.

数值模拟是对此类灾害链运动过程进行反演和预测分析的重要手段.近几十年已经发展出了较为成熟的可模拟单独滑坡堵河、堰塞坝溃决、洪水运动演化过程的预测模型，现有研究通过集成这些模型实现了滑坡灾害链模拟，取得了显著进展，但在大流域灾害链模拟方面仍存在一些局限性.在滑坡堵河与溃决过程模拟中，已有研究使用DAN^3D和FLOW^3D模拟滑坡和洪水演进过程，较为精细地反演了滑坡堵河和堰塞坝的溃决过程（戴兴建等 2019）；或使用离散元法（DEM）模拟滑坡过程，之后利用光滑粒子流法（SPH）分析了堰塞坝的溃决过程（Zhu et al.， 2023）；或使用水土耦合的SPH方法模拟滑坡-堵江-成坝灾害链全过程，将侵蚀和侧坡失效与SPH耦合，可以更精确地模拟坝体溃决过程（Li et al.， 2024； 李爽等， 2025）.然而，此类方法需对溃决过程开展精细模拟，计算量大，限制了其在大规模滑坡灾害链模拟中的应用.为提高计算效率，部分研究利用PFC3D进行滑坡动力学分析，并通过HEC-RAS模拟滑坡坝溃决后的洪水过程（贾珂程等， 2023）；或采用Massflow对滑坡的运动过程进行模拟，之后结合DABA模型模拟堰塞坝的溃决过程，最后利用HEC-RAS模型模拟溃坝洪水过程（Fan et al.， 2020）；亦或使用CFD方法模拟滑坡坝溃决过程，可准确刻画坝体受洪水侵蚀而引发的阶段性演化规律（Li et al.， 2025）.此类方法可快速估算滑坡坝溃决洪水的演进过程，适用于中小尺度滑坡灾害链的分析，但对堰塞坝的实际溃决过程考虑尚不充分.近年来，有研究建立了考虑尺寸选择性侵蚀的滑坡坝溃决模型（Ma et al.， 2025）；还有研究利用GPU加速的ED-SHP模拟土石坝漫顶破坏过程中的侵蚀和沉积现象（Zhou et al.， 2024）；此类方法可以精确模拟侵蚀沉积与侧岸侵蚀，但仅适用于小尺度的精细模拟.大尺度模拟上，有研究使用双层深度平均模型，通过Savage-Hunter类型的深度平均模型描述滑坡运动，并结合浅水方程模拟洪水演进（刘威和何思明， 2020）.此类方法的计算效率较高，理论上适用于大尺度滑坡灾害链模拟，但现有模型采用的颗粒侵蚀沉积算法多基于经验，若能在其中考虑侵蚀沉积的物理机制将有助于提高其精度和适用性.

基于以上现状，本研究通过耦合深度平均滑坡-堵河运动学模型和考虑堰塞坝物理溃决过程的洪水动力学模型，建立了适用于大流域巨型滑坡-堰塞湖-溃决洪水灾害链演进过程模拟的运动学模型，并以德恒隆滑坡为例，模拟分析了其堵河成坝、堰塞坝溃决和洪水演进的全过程，量化了此类灾害的链式放大效应，为黄河流域滑坡灾害链的风险预警与防控提供依据.

1 研究区概况

研究区位于青藏高原东北缘黄河上游龙羊峡至石嘴山段（图1a），根据灾种类型可将其划分为堰塞湖段、滑坡段、洪水段三部分，全长约1 000 km.该区域黄河总体呈西北-东南走向，地势西北高、东南低，地跨青海、甘肃、宁夏三省（区）.其中，滑坡段位于青海省黄南藏族自治州尖扎县境内，地处公伯峡水库库尾右岸，地理坐标为北纬35°56′、东经102°04′，属黄河上游峡谷区向宽谷盆地的过渡带.堰塞湖段（龙羊峡-公伯峡）以深切峡谷为主（如龙羊峡、李家峡、公伯峡等），河谷狭窄，岸坡陡峭（纪一奇， 2024）.洪水段（公伯峡-石嘴山）过渡为宽谷盆地（如兰州盆地、银川平原），河道蜿蜒，发育多级阶地，局部形成河漫滩与冲积平原.正是此区域上游河道窄地势高下游河道宽地势低的地形，为德恒隆滑坡诱发堵河-堰塞湖溃决-洪水链生灾害提供了条件.

德恒隆滑坡发生于约8万a前，形成堰塞湖并于滑坡发生3 000 a后溃决（郭小花等， 2015； 白世彪等， 2021）.郭小花等（2015）曾对德恒隆滑坡的特征及形成机制进行详细的研究，其研究表明德恒隆滑坡主滑方向为WSW255°，滑坡后缘高程为2 940 m，高出黄河（2 040 m）900 m.现黄河河道两岸仍清晰可见大量残留的堰塞坝堆积体，残留体积约14亿m³，坝长东北岸约2 150 m，西南岸约4 000 m，厚度左岸约120 m，右岸约25 m.卫星影像解译（图2）和现场调查结果表明，该滑坡位于尖扎盆地出口处，滑距约5 km，根据残余地形与郭小花等复原滑坡剖面估计滑坡滑前总体积约35亿m³，属于巨型滑坡.滑坡启动区位于黄河东北岸，启动后越过古黄河，与河对岸山体发生碰撞，完全堵塞河道并形成堰塞湖，滑坡堰塞坝长约4.4 km，高约234 m.堰塞湖沿滑坡坝中前部发生溃决，后在河水冲刷侵蚀下逐渐形成如今U字形的河道.

滑坡区地质图和滑坡主断面地质剖面图如图1所示.德恒隆滑坡主体为新太古-古元古界片麻岩，河对岸基岩为全新统砂岩.滑坡后缘沿一断层发育，断层面产状为255°∠48°.砂岩与片麻岩沿另一断层面呈不整合接触，断层面产状为72°∠65°（郭小花等， 2015）.

2 模型构建

针对滑坡堵河过程，基于浅水理论建立了其运动学模型：由于巨型滑坡在高速运移中因剧烈剪切作用呈现类流体特性，采用等效流体假设与深度平均法（DAM）将三维问题简化为二维平面运动，通过质量守恒方程和动量守恒方程描述滑体厚度、速度与地形的动态耦合.这一简化在保证物理机制真实性的同时显著提升了计算效率，结合高精度Godunov型有限体积法，可精确捕捉滑体前缘激波和干湿边界，并利用重采样至50 m的三维地形数据平衡了地貌细节保留与计算资源消耗，为预测滑坡坝的几何形态与堵河风险提供了可靠工具.

2.1　滑坡过程模拟

滑坡堵河过程运动学模型的控制方程由滑体质量守恒和动量守恒方程组成（Xia and Liang， 2018），该模型建立在全局笛卡尔坐标系下，采用深度积分法推导得出，模型考虑了复杂弯曲地形对滑坡运动过程的影响.控制方程如下：

\frac{\partial q}{\partial t}+\frac{\partial f\left(q\right)}{\partial x}+\frac{\partial g\left(q\right)}{\partial y}=S_{\mathrm{b}}+S_{\mathrm{f}} ， (1)

式中：q为与流量变量相关的矩阵，f（q）和g（q）是x和y方向上的流量矩阵，S_b是地形项，S_f是摩擦项，分别表示为：

q=\left[\begin{array}{c}h\\ uh\\ vh\end{array}\right]f\left(q\right)=\left[\begin{array}{c}uh\\ u^{\mathrm{2}}h+\frac{\mathrm{1}}{{\varphi }^{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{g}{h}^{\mathrm{2}}\\ uvh\end{array}\right]\begin{array}{c}\begin{array}{c}g\left(q\right)\end{array}\end{array}=\left[\begin{array}{l}vh\\ uvh\\ v^{\mathrm{2}}h+\frac{\mathrm{1}}{{\varphi }^{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{g}{h}^{\mathrm{2}}\end{array}\right] ， (2)

S_{\mathrm{b}}=\left[\begin{array}{c}\mathrm{0}\\ -ah\frac{\partial b}{\partial x}+\frac{\mathrm{1}}{\mathrm{2}}\mathrm{g}{h}^{\mathrm{2}}\frac{\partial \left(\mathrm{1}/{\varphi }^{\mathrm{2}}\right)}{\partial x}\\ -ah\frac{\partial b}{\partial y}+\frac{\mathrm{1}}{\mathrm{2}}\mathrm{g}{h}^{\mathrm{2}}\frac{\partial \left(\mathrm{1}/{\varphi }^{\mathrm{2}}\right)}{\partial y}\end{array}\right] ， (3)

S_{\mathrm{f}}=\left[\begin{array}{c}\mathrm{0}\\ \frac{\mu ahu\varphi +c/\rho }{\sqrt[]{u^{\mathrm{2}}+v^{\mathrm{2}}+{\left(u\frac{\partial b}{\partial x}+v\frac{\partial b}{\partial y}\right)}^{\mathrm{2}}}}\\ \frac{\mu ahv\varphi +c/\rho }{\sqrt[]{u^{\mathrm{2}}+v^{\mathrm{2}}+{\left(u\frac{\partial b}{\partial x}+v\frac{\partial b}{\partial y}\right)}^{\mathrm{2}}}}\end{array}\right] ， (4)

式中：h为滑体厚度，u和v是x和y方向上滑体的深度平均速度，g是重力加速度，b是河床高程，μ是摩擦系数，c是滑体粘聚力，ρ是滑体的密度，a是与弯曲地形离心效应相关的参数，反映复杂地形曲率引起的离心效应，ϕ是与滑面地形坡度相关的参数，反映地形坡度对滑面支持力的影响.a和ϕ的表达式如下：

a=\frac{\mathrm{1}}{{\varphi }^{\mathrm{2}}}\left(\mathrm{g}+V^{T}HV\right) ， (5)

\varphi =\left\{{\left(\frac{\partial b}{\partial x}\right)}^{\mathrm{2}}+{\left(\frac{\partial b}{\partial y}\right)}^{\mathrm{2}}+\mathrm{1}\right\} ， (6)

式中：

V={\left(u,v\right)}^T ， (7)

H=\left[\begin{array}{cc}\frac{{\partial }^{\mathrm{2}}b}{\partial x^{\mathrm{2}}}& \frac{{\partial }^{\mathrm{2}}b}{\partial x\partial y}\\ \frac{{\partial }^{\mathrm{2}}b}{\partial x\partial y}& \frac{{\partial }^{\mathrm{2}}b}{\partial y^{\mathrm{2}}}\end{array}\right] . (8)

2.2　溃坝洪水过程模拟

滑坡坝堵塞河道后，上游来水在坝后淤积形成堰塞湖.随着堰塞湖水位上升并逐渐漫过坝顶，滑坡坝在流水侵蚀作用下可能发生溃决（梅胜尧等， 2023； 杨泽文等， 2025）.绝大多数堰塞湖以漫顶溃坝的形式发生溃决，在坝顶洪水的下切和侧蚀作用下，坝体形成溃口，上游湖水通过溃口迅速排放，形成洪水波及下游.由于溃决洪水的传播距离和宽度远大于其深度，本项目采用深度平均法建立其运动学控制方程，并利用有限体积法构建溃决洪水动力学模型.模型将堰塞坝视为可侵蚀物，引入侵蚀与沉积项以考虑洪水对坝体颗粒物质的搬运过程，可模拟堰塞坝的漫顶溃决过程，预测洪水的流量、流速及其在下游区域的传播和淹没范围.

溃决过程洪水动力学模型的控制方程可根据流体的质量守恒和动量守恒关系以及颗粒的质量守恒关系推导得到（Xia et al.， 2023），表示成矩阵形式如下：

\frac{\partial q}{\partial t}+\frac{\partial f\left(q\right)}{\partial x}+\frac{\partial g\left(q\right)}{\partial y}=S_{\mathrm{b}}+S_{\mathrm{f}}+S_{\mathrm{e}} ， (9)

式中：q为流体变量矩阵，f和g分别为x和y方向的流量向量， S_b、S_f分别为与河床地形和流动阻力相关的源项，S_e为与侵蚀和沉积作用相关的源项，这些矩阵可展开如下：

q=\left[\begin{array}{c}{h}_{\mathrm{w}}\\ h_{\mathrm{w}}u\\ h_{\mathrm{w}}v\\ \begin{array}{c}{h}_{\mathrm{w}}C\\ h_s\end{array}\end{array}\right]f\left(q\right)=\left[\begin{array}{c}{h}_{w}u\\ h_{\mathrm{w}}{u}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}\mathrm{g}{h_{\mathrm{w}}}^{\mathrm{2}}\\ h_{\mathrm{w}}uv\\ \begin{array}{c}{h}_{\mathrm{w}}uC\\ \mathrm{0}\end{array}\end{array}\right]g\left(q\right)=\left[\begin{array}{c}{h}_{w}v\\ h_{\mathrm{w}}uv\\ h_{\mathrm{w}}{v}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}\mathrm{g}{h_{\mathrm{w}}}^{\mathrm{2}}\\ \begin{array}{c}{h}_{\mathrm{w}}vC\\ \mathrm{0}\end{array}\end{array}\right] ， (10)

S_{\mathrm{b}}=\left[\begin{array}{c}\mathrm{0}\\ -\mathrm{g}{h}_{\mathrm{w}}\frac{\partial b}{\partial x}\\ -\mathrm{g}{h}_{\mathrm{w}}\frac{\partial b}{\partial y}\\ \begin{array}{c}\mathrm{0}\\ \mathrm{0}\end{array}\end{array}\right]S_{\mathrm{f}}=\left[\begin{array}{c}\mathrm{0}\\ \frac{{\tau }_x}{\rho }\\ \frac{{\tau }_y}{\rho }\\ \begin{array}{c}\mathrm{0}\\ \mathrm{0}\end{array}\end{array}\right] ， (11)

S_{\mathrm{e}}=\left[\begin{array}{c}E+i\\ -\frac{\left({\rho }_{\mathrm{s}}-{\rho }_{\mathrm{f}}\right)\mathrm{g}{h_{\mathrm{w}}}^{\mathrm{2}}}{\mathrm{2}\rho }\frac{\partial C}{\partial x}+\frac{\left({\rho }_{\mathrm{s}}-{\rho }_{\mathrm{f}}\right)uE}{\rho }\\ -\frac{\left({\rho }_{\mathrm{s}}-{\rho }_{\mathrm{f}}\right)\mathrm{g}{h_{\mathrm{w}}}^{\mathrm{2}}}{\mathrm{2}\rho }\frac{\partial C}{\partial y}+\frac{\left({\rho }_{\mathrm{s}}-{\rho }_{\mathrm{f}}\right)vE}{\rho }\\ \begin{array}{c}(\mathrm{1}-e)E\\ -E\end{array}\end{array}\right] ， (12)

式中：h_w和h_s分别表示洪水深度和可侵蚀体厚度，u和v表示x和y轴的深度平均速度，E表示侵蚀沉积速率，C表示水中固体颗粒浓度，e表示可侵蚀物的孔隙比，τ是洪水基底阻力，g是重力加速度，ρ是洪水的密度，ρ_s是可侵蚀物的密度，ρ_f 是水的密度.

当剪切力τ大于抗侵蚀力τ_c时，可侵蚀物将会被洪水侵蚀，并在洪水中扩散，反之，当洪水不足以搬运其中的颗粒物质时，颗粒物会在河床发生沉积.侵蚀和沉积速率可表示为：

， (13)

\tau =\rho \mathrm{g}{n}^{\mathrm{2}}{h_{\mathrm{w}}}^{-\frac{\mathrm{1}}{\mathrm{3}}}\left(u^{\mathrm{2}}+v^{\mathrm{2}}\right)+\left({\rho }_s-{\rho }_f\right)\mathrm{g}{h}_{\mathrm{w}}C\mathrm{t}\mathrm{a}\mathrm{n}{\phi }_{\mathrm{d}} ， (14)

{\tau }_{\mathrm{0}}=\left({\rho }_s-{\rho }_f\right)\mathrm{g}{h}_{\mathrm{w}}C\mathrm{t}\mathrm{a}\mathrm{n}{\phi }_{\mathrm{s}} . (15)

式中：α和β为无量纲因子，分别用于调节侵蚀和沉积速率，τ为施加在可侵蚀物上的剪应力，τ₀为侵蚀阈值应力，表示可侵蚀物的抗侵蚀力，n为曼宁系数，φ_d为水流与基底之间的动摩擦角，φ_s是可侵蚀物的静摩擦角.

3 滑坡堵河灾害放大效应

为反演德恒隆滑坡堵河并形成堰塞湖的过程，参考郭小花等（2015）复原的滑坡剖面，并结合现存的地质特征，如滑面特征、剪出口位置、周边地形特征、等高线特征以及滑移后的形态、坝高、堰塞坝位置等，反推还原了滑坡前的初始滑体形态，据此建立滑坡运动学模型.

为了得到最佳的滑坡形态，开展了不同工况下的滑坡堵河过程模拟.段颖星和张天芳（2022）指出，岩体密度是诸多因素（地震加速度、内摩擦角、粘聚力、库水位变化、岩体密度）中敏感性最小的，其对滑坡的影响主要是通过质量影响重力势能，但实际运动中，势能在转化为动能的过程中克服摩擦阻力而被消耗，导致密度的影响被削弱.因此，本文仅将摩擦系数（μ）和粘聚力（c）作为工况变量.德恒隆滑坡滑体为片麻岩且多破碎，其密度参考了唐家山滑坡（李敬和苏爱军， 2018），具体参数选取见表1.

图3为不同摩擦系数和粘聚力下得到的模拟结果.结果表明，摩擦系数和粘聚力对滑坡堵河的过程和形态均有一定影响.随着摩擦系数的减小，滑体滑移距离增大，物质向前缘堆积.随着粘聚力增大，滑坡滑移距离减小，且后部堆积变厚.通过与现存地形进行对比，发现μ=0.11时其滑移距离较远，坝体位置较近；而μ=0.13和μ=0.14时，由于摩擦系数较大，二者整体滑移距离较近，堆积靠后，坝体位置靠前；在μ= 0.12时，滑移距离和坝体位置对应较好，且不同粘聚力下其结果相近，但在坝体形态上仍有差异.当c=20 kPa时，其坝体形态较其他两种工况堆积体的厚度及分布情况与实际最为接近，且坝体下凹更加明显，便于后续的溃坝过程模拟，如图2中剖面图所示.在此工况下，溃决处坝体高度234 m，黄河改道约向西移动800 m，与实测接近.

图4为滑体的平均速度曲线.从曲线可以看出，滑体启动后迅速加速，27 s时就达到了44.5 m/s的峰值速度.随后，平均速度以相对加速过程稍缓的速率下降，在75 s时，降速进一步减缓，并在120 s时基本停止运动.

根据模拟结果，滑体启动后前缘首先开始滑动，24 s时，沟道两侧滑体发生碰撞，滑坡整体下移，跨越黄河，覆盖了现今尖扎县部分区域，最大厚度由812.8 m锐减至550.9 m；滑体前部区域堆积加厚，滑体厚度由后缘向前缘逐渐递减，后缘的厚度变化较小，而滑体速度则呈现从后缘到前缘递增的趋势，前缘中部的最大速度达到117.0 m/s.48 s时，滑体继续前移，此时已覆盖尖扎县的大部区域，最大厚度降至410.3 m；滑体前缘接近滑移的最远距离，前部逐渐堆积，滑体后缘较厚的区域下降约130 m；此时最大速度位于滑体中部边缘.72 s时，滑体前缘已基本稳定，最大厚度为305.4 m；前缘中部的少量滑体继续向前扩展；此时，最大速度为70.1 m/s. 96 s时，滑体前部已趋于稳定，最大厚度为280.0 m，前部最大厚度为187 m；滑体后部少量滑体仍在滑动.120 s时，滑体已基本稳定，形成堰塞坝高度约234 m.

滑坡坝形成后，由于河道堵塞，往往会在其上游形成堰塞湖，淹没较大区域并造成一定危害.而不同规模的滑坡堰塞坝形成的堰塞湖大小不同，对上游的危害程度也不同.为建立滑坡堰塞坝坝高与堰塞湖规模的关系，假设在德恒隆滑坡处形成不同坝高的堰塞坝，通过GIS软件分析了不同坝高对应的堰塞湖面积和体积.形成的堰塞湖区域面积和体积趋势曲线如图5a和5b所示.

堰塞湖的面积和体积与坝高关系的拟合公式如下：

S=\mathrm{4}\mathrm{237.02}{h}^{\mathrm{2}}+\mathrm{204}\mathrm{198}h ， (16)

V=\mathrm{1}\mathrm{600.3}{h}^{\mathrm{3}}-\mathrm{3}\mathrm{953.6}{h}^{\mathrm{2}}+\mathrm{1.515}\mathrm{0}\times {\mathrm{10}}^{\mathrm{7}}h.(17)

式中：S为堰塞湖的面积，m²；V为堰塞湖的体积，m³；h为堰塞坝的高度，m.

当坝高为100 m时，堰塞湖体积为30.9亿 m³.当坝高增加到150 m时，堰塞湖体积为76.2 亿 m³.在德恒隆滑坡实际坝高（约234 m）的条件下，计算得到的堰塞湖面积约280 km²，体积约232.3亿 m³.据此趋势估计，当坝高达到300 m时，堰塞湖体积可达475亿 m³，是当前模拟堰塞湖体积的2倍.

4 堰塞坝溃决灾害放大效应

4.1　堰塞坝溃决过程

由于堰塞湖形成至溃坝时间跨度过长，为节省计算量，假设在模拟时堰塞湖已蓄满，即以堰塞湖水深为233.4 m作为初始条件，同时在计算域上游河道边界处，按黄河上游年均径流量900 m³/s施加常流量边界，启动堰塞湖溃决过程模拟.模拟过程中，在坝体的前、中、后各设置1个监测点，用以监测堰塞坝的溃决过程.其中监测点S1同时也是坝后侵蚀速率监测点.

模型模拟所需基本参数取值见表2.动静摩擦系数和孔隙比参考唐家山堰塞坝案例取值，静摩擦系数为0.577，动摩擦系数取静摩擦系数的0.8倍，孔隙比为0.32（李敬和苏爱军， 2018）.α和β取值参考Xia et al.（2023）的推荐值.曼宁系数参考黄河上游河渠的取值（王慧敏等， 2024）.

模拟所得溃口流量曲线如图6所示.Peng和Zhang基于大量案例数据，提出了经典的堰塞坝峰值流量预测Peng和Zhang公式，被国内外学者广为

采用，其预测精度得到了大量案例验证，计算结果和实测数据吻合良好（Peng and Zhang， 2012）.本文计算峰值流量与其吻合，表明模拟结果基本合理.Peng和Zhang模型使用坝高、堰塞湖库容以及侵蚀率系数等估计堰塞坝溃口的峰值流量，计算公式如下：

Q_{\mathrm{p}}={\mathrm{g}}^{\mathrm{1}/\mathrm{2}}{H}_{\mathrm{d}}^{\mathrm{5}/\mathrm{2}}{\left(\frac{h_{\mathrm{d}}}{H_{\mathrm{r}}}\right)}^{-\mathrm{1.371}}{\left(\frac{V_{\mathrm{l}}^{\mathrm{1}/\mathrm{3}}}{H_{\mathrm{d}}}\right)}^{\mathrm{1.536}}{\mathrm{e}}^a,(18)

式中：Q_p为溃决洪峰流量；H_d为坝高，m；V_l为堰塞湖库容，m³；g为重力加速度，m/s²；H_r为单位长度，1 m；e为自然常数；a为侵蚀率系数，与坝体材料性质有关，坝体可蚀性为高、中、低时，a分别取1.236、-0.380、-1.615.各参数取值见表3.该模型计算得到的溃口峰值流量为23.76万m³/s，而模拟值为23.8万m³/s，二者仅相差0.31%.

为定量研究堰塞湖的溃决特征，将堰塞湖开始溃决后任意时刻的最大水深与初始最大水深的比值定义为溃决程度B_d，计算公式如下：

B_{\mathrm{d}}=\left(\mathrm{1}-\frac{h_{\mathrm{t}}}{h_{\mathrm{m}\mathrm{a}\mathrm{x}}}\right)\times \mathrm{100}\mathrm{\%} ， (19)

式中：h_t表示t时刻堰塞湖最大水深，h_max表示堰塞湖溃决前最大水深.

根据模拟结果，绘制堰塞坝溃决程度时程曲线（图7a）及3处监测点（图7b）的侵蚀速率时程曲线.结果表明，堰塞坝溃决过程持续约100 h.溃决开始，坝体后部最先受到侵蚀，并逐渐扩展至中前部；坝后侵蚀速率逐步增大，至31 h达到7.4 m/h的峰值.随着溃口向上游扩展，坝中16 h开始发生剧烈侵蚀，至19 h时最大侵蚀速率达到8.2 m/h. 21 h时，侵蚀延伸至坝体前部，侵蚀通道贯通.25 h，由于前部和中部的侵蚀物堆积，测得坝前监测点的侵蚀速率减缓.60 h，坝后和坝中监测点处的坝体已基本稳定.85 h，坝前有少量侧岸坍塌，至100 h坝体侵蚀过程基本结束.总体而言，坝体前部侵蚀过程最快，中部次之，后部最慢但持续时间最长.

图8展示了堰塞坝不同时刻的溃决形态.T=0 h时，堰塞湖水位正在上升，尚未漫过堰塞坝.T=8 h时，水流刚漫过坝顶，溃决过程开始.T=16~30 h时，堰塞坝发生快速溃决，溃口最初在坝体后部形成，之后在水流侵蚀下迅速下切并向上游扩展，同时两侧出现坍塌，导致溃口加宽；溃口处已侵蚀的滑体被水流带走，并在下游不远处沉积.T=40 h时，溃口基本贯通，形成上游窄下游宽的喇叭形溃决通道，通道中水位开始下降.T=300 h时，溃口水流趋于稳定，水位缓慢下降.溃决后的河道为凸向滑动方向的弧形，与现今黄河河道位置接近，验证了模拟结果的合理性.

4.2　洪水演进过程

为监测洪水变化过程，在模拟全域共设置23个监测点和监测剖面，如图9所示.监测点位于监测剖面上的地势最低处.图中蓝线为盆地处的监测剖面，红线为各盆地出口处的沟道剖面.

由图10可知，盆地和沟道处的峰值流量随剖面到溃口距离增大均呈现明显的衰减趋势.溃决过程从0 h开始，13 h时流量迅速上升，至33 h达到峰值流量，为23.8万m³/s，后迅速下降，至100 h降至0.94万m³/s. 59 h时洪水前缘到达剖面3，该处接近兰州市，峰值流量仍达到5.9 万m³/s. 242 h时到达接近银川市的剖面6，峰值流量降为2.1万m³/s. 403 h到达出口处，峰值流量为1.0万m³/s，该处流量仍为黄河年平均流量0.14万m³/s （采用2021年水利部公开数据计算）的7倍左右.

根据剖面流量统计结果，绘制洪水峰值流量-距离关系曲线，如图11.拟合结果显示，洪水峰值流量随距离增大以幂指数形式衰减，可表达为：

q=\mathrm{5.546}\mathrm{4}\times {\mathrm{10}}^{\mathrm{6}}{D_{\mathrm{f}}}^{-\mathrm{0.853}\mathrm{9}} ， (20)

式中：q表示峰值流量，D_f为到溃决口的距离.

根据模拟结果绘制洪水淹没区域的历史最大水深图，如图9.由图可知，沿滑坡坝向下游，洪水的最大淹没深度有减小的趋势：循化盆地最高水位达127 m，位于盆地出口；受地形影响官亭盆地最高水位135 m，位于盆地出口；喇家遗址处最高水位77.07 m；在兰州市范围内，洪水最高水位仍有30.43 m；到中卫市时，最高水位降7.26 m；而银川市最高水位进一步衰减至3.47 m；抵达银川盆地出口时，最高水位降至2 m左右，洪水影响基本消减.由此可见，相比于滑坡本身的运动距离（约5 km），本次巨型滑坡-堵河-溃坝洪水灾害链形成的堰塞湖淹没至上游90 km尼那水库出水口处，溃决洪水更是波及了坝体下游超过900 km河道范围内的广大区域，具有显著的灾害放大效应.

5 结论

本文通过模拟德恒隆滑坡-堰塞湖-溃决洪水灾害链的运动过程，分析了滑坡的运动特征及其堵塞黄河形成堰塞湖的规模，反演了洪水演进的过程，得到以下结论：

（1）德恒隆滑坡持续了120 s，其在27 s时达到最大平均速度44.5 m/s，为高速滑坡.滑坡最大堆积厚度为289.4 m，堵塞黄河河道，形成堰塞坝高达234 m.在上游形成面积约280 km²，体积约232.3亿 m³的堰塞湖.

（2）参数敏感性分析结果表明，摩擦系数和粘聚力对滑坡堵河的过程和形态均有一定影响.随着摩擦系数的减小，滑体滑移距离增大，物质向前缘堆积.随着粘聚力增大，滑坡滑移距离减小，且后部堆积变厚.

（3）溃坝持续约100 h，洪水演进过程超过600 h，高水位超过150 m，波及下游约900 km流域范围，具有显著的灾害放大效应.

（4）狭窄的地形会对洪水灾害进一步放大.洪水流量在演进过程整体呈现出逐渐下降的趋势，但由于地形变化，洪水会在盆地处流速减缓，水深降低，但在经过盆地出口时洪水流速增加，水深升高，危害性增加.

（5）构建了深度平均滑坡-堵河运动学模型与考虑颗粒物运移及物理侵蚀沉积机制的溃决洪水动力学耦合模型，模型对此类巨型滑坡灾害链的运动过程模拟具有较好的适用性.

参考文献

[1] Bai, S. B., Cui, P., Ge, Y. G., 2021. ¹⁰Be and ²⁶Al Exposure Dating Study of the Dehenglong Giant Landslide in the Upper Reaches of the Yellow River. In Proceedings of the 2021 Chinese Geoscience Union Meeting (Volume 37), 97-100 (in Chinese).

[2] Dai, X. J., Yin, Y. P., Xing, A. G., 2019. Simulation and Dynamic Analysis of Yigong Rockslide-Debris Avalanche-Dam Breaking Disaster Chain. The Chinese Journal of Geological Hazard and Control, 30(5): 1-8 (in Chinese with English abstract).

[3] Duan, Y.X., Zhang, T.F., 2022. The Effects of Interaction on Landslide Stability. Statistics and Application, 11(6): 1313-1317 (in Chinese with English abstract).

[4] Fan, X. M., Yang, F., Siva Subramanian, S., et al., 2020. Prediction of a Multi-Hazard Chain by an Integrated Numerical Simulation Approach: The Baige Landslide, Jinsha River, China. Landslides, 17(1): 147-164. https://doi.org/10.1007/s10346-019-01313-5

[5] Guo, X. H., 2017. Relationship between Landslide-Dammed Lakes with the Evolution of Upper Stream of The Yellow River from Longyang Gorge to Liujia Gorge on Northeast Margin of Tibet Plateau (Dissertation). Chang’an University, Xi’an (in Chinese with English abstract).

[6] Guo, X. H., Lu, Y. D., Li, X. L., et al., 2015. Event of Block up of Upper Yellow River by Dehenglong-Suozi Landslides. Journal of Jilin University (Earth Science Edition), 45(6): 1789-1797 (in Chinese with English abstract).

[7] Ji, Y. Q., 2024. Tectonic Geomorphology and Genetic Mechanism of Giant Landslides from Lagan Gorge to Liujia Gorge in the Upper Reaches of the Yellow River (Dissertation). Chang’an University, Xi’an (in Chinese with English abstract).

[8] Jia, K. C., Zhuang, J. Q., Zhan, J. W., et al., 2023. Reconstruction of the Dynamic Process of the Holocene Gelongbu Landslide-Blocking-Flood Geological Disaster Chain Based on Numerical Simulation. Earth Science, 48(9): 3402-3419 (in Chinese with English abstract).

[9] Jin, Z., Yu, X., Li, Z. H., et al., 2024. Scientific Investigation on the Entire Yellow River Basin: Understanding and Reflections. China Water Resources, (23): 7-11, 6 (in Chinese with English abstract).

[10] Lan, H. X., Peng, J. B., Zhu, Y. B., et al., 2022. Research on Geological and Surfacial Processes and Major Disaster Effects in the Yellow River Basin. Scientia Sinica (Terrae), 52(2): 199-221 (in Chinese).

[11] Li, J., Su, A. J., 2018. Numerical Simulation of Dam-Break Flood for Tangjiashan Landslide Dam. Water Resources and Power, 36(8): 56-59, 147 (in Chinese with English abstract).

[12] Li, S., Peng, M., Gao, L., et al., 2024. A 3D SPH Framework for Simulating Landslide Dam Breaches by Coupling Erosion and Side Slope Failure. Computers and Geotechnics, 175: 106699. https://doi.org/10.1016/j.compgeo.2024.106699

[13] Li, S., Peng, M., Shi, Z. M., et al., 2025. Simulation and Analysis of Cascading Hazard Based on Fluid-Soil Coupled SPH Method. Earth Science, 50(10): 3967-3981 (in Chinese with English abstract).

[14] Li, X. L., Guo, X. H., Li, W. H., 2011. Mechanism of Giant Landslides from Longyangxia Valley to Liujiaxia Valley along Upper Yellow River. Journal of Engineering Geology, 19(4): 516-529 (in Chinese with English abstract).

[15] Li, Y. C., Guo, S. F., Zheng, B. W., et al., 2025. Analysis of Ancient Rongcharong Landslide Dam Failure Events in the Suwalong Reach of the Upper Reaches of the Jinsha River. Journal of Earth Science, 36(5): 2005-2022. https://doi.org/10.1007/s12583-024-0052-x

[16] Liu, W., He, S. M., 2020. Numerical Simulation of the Evolution Process of Disaster Chain Induced by Potential Landslide in Woda of Jinsha River Basin. Advanced Engineering Sciences, 52(2): 38-46 (in Chinese with English abstract).

[17] Ma, C. Y., Peng, M., Yang, J. T., et al., 2025. A Breaching Model for Landslide Dams Involving Size-Selective Erosion and Breach Evolution. Engineering Geology, 355: 108235. https://doi.org/10.1016/j.enggeo.2025.108235

[18] Mei, S. Y., Zhong, Q. M., Chen, S. S., et al., 2023. Numerical Simulation of Breach Hydrograph and Morphology Evolution during Landslide Dam Breaching. Earth Science, 48(4): 1634-1648 (in Chinese with English abstract).

[19] Peng, J. B., Ma, R. Y., Lu, Q. Z., et al., 2004. Geological Hazards Effects of Uplift of Qinghai-Tibet Plateau. Advance in Earth Sciences, 19(3): 457-466 (in Chinese with English abstract).

[20] Peng, J. B., Wang, F. Y., Xu, J. S., 2025. Basin Layered Structure and Urban Geological Safety. Earth Science, 50(8): 3117-3131 (in Chinese with English abstract).

[21] Peng, M., Zhang, L. M., 2012. Breaching Parameters of Landslide Dams. Landslides, 9(1): 13-31. https://doi.org/10.1007/s10346-011-0271-y

[22] Wang, H. M., Luan, J. K., Ma, N., 2024. Hydrodynamic Modelling of Flash Flood Processes in the Typical Regions in the Upper Yellow River Basin. Journal of Beijing Normal University (Natural Science), 60(4): 519-528 (in Chinese with English abstract).

[23] Xia, X. L., Jarsve, K. T., Dijkstra, T., et al., 2023. An Integrated Hydrodynamic Model for Runoff-Generated Debris Flows with Novel Formulation of Bed Erosion and Deposition. Engineering Geology, 326: 107310. https://doi.org/10.1016/j.enggeo.2023.107310

[24] Xia, X. L., Liang, Q. H., 2018. A New Depth-Averaged Model for Flow-Like Landslides over Complex Terrains with Curvatures and Steep Slopes. Engineering Geology, 234: 174-191. https://doi.org/10.1016/j.enggeo.2018.01.011

[25] Yang, Z. W., Wu, B. B., Liu, W. M., et al., 2025. Progress in Erosion Mechanism and Geomorphological Effects of High-Energy Outburst Floods. Earth Science, 50(2): 718-736 (in Chinese with English abstract).

[26] Yin, Z. Q., Wei, G., Qin, X. G., et al., 2021. Research Progress on Landslides and Dammed Lakes in the Upper Reaches of the Yellow River, Northeastern Tibetan Plateau. Earth Science Frontiers, 28(2): 46-57 (in Chinese with English abstract).

[27] Zhou, M. J., Shi, Z. M., Peng, C., et al., 2024. Two-Phase Modelling of Erosion and Deposition Process during Overtopping Failure of Landslide Dams Using GPU-Accelerated ED-SPH. Computers and Geotechnics, 166: 105944. https://doi.org/10.1016/j.compgeo.2023.105944

[28] Zhu, Q. Y., Jiang, N., Chen, Q., et al., 2023. DEM-SPH Simulation for the Formation and Breaching of a Landslide-Dammed Lake Triggered by the 2022 Lushan Earthquake. Landslides, 20(9): 1925-1941. https://doi.org/10.1007/s10346-023-02096-6

基金资助

长安大学中央高校基本科研业务费专项资金项目(300102264203)

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

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

中国博士后科学基金项目(2023M730353)