基于水土耦合SPH方法的滑坡-堵江-成坝灾害链全过程动力演化模拟

Simulation and Analysis of Cascading Hazard Based on Fluid-Soil Coupled SPH Method

摘要

采用一种双向耦合的SPH数值模型，精确模拟滑坡-堵江-成坝灾害链的全过程.模型以Drucker–Prager准则描述滑体大变形行为，结合混合物理论与非线性渗流拖曳力实现水土耦合.通过室内试验验证后，成功重现白格滑坡灾害链演化，模拟结果与实测高度吻合.结果表明，滑坡入水引发的涌浪及成坝过程可依据滑体速度与能量变化清晰划分阶段.量化分析显示，内摩擦角φ增大（5°~20°）导致堰塞坝长度线性减小，高度呈幂函数增长，涌浪峰值高度显著降低.涌浪峰值与滑体入水弗劳德数呈线性正相关.上述发现揭示了滑体参数对灾害链演化路径的系统性影响，为高风险山地河流域灾害预测与风险评估提供理论支撑.

Abstract

This study adopts a bidirectionally coupled SPH numerical model to accurately simulate the full evolution of a landslide-dammed lake disaster chain. The model captures large deformation of the landslide body using the Drucker-Prager criterion and achieves water–soil coupling through mixture theory and nonlinear seepage drag forces. Validated against laboratory experiments, the model successfully reproduces the Baige landslide disaster chain, with simulation results closely matching field observations. Results show that the processes of landslide motion, impulse wave generation, and dam formation can be clearly delineated by the evolution of landslide velocity and energy. Quantitative analysis reveals that increasing the internal friction angle φ from 5° to 20° leads to a linear decrease in dam length, a power-law increase in dam height, and a significant reduction in wave height. The peak wave height exhibits a linear correlation with the landslide Froude number at impact. These findings highlight the systematic influence of landslide material properties on disaster chain dynamics and offer theoretical support for hazard prediction and risk assessment in mountainous river basins.

Graphical abstract

关键词

堰塞坝 / 滑坡灾害链 / 光滑粒子流体动力学 / 水土耦合 / 工程地质学.

Key words

landslide dam / landslide hazard chain / smoothed particle hydrodynamics / fluid-soil coupling / engineering geology

0 引言

滑坡-堵江-成坝是由滑坡体失稳运移，滑坡体入江涌浪-滑坡体堆积成堰塞坝构成的一种天然地质灾害链，广泛分布于我国的西部山区（Fan et al.， 2019）.除滑坡体本身运动过程中对沿途建设的损害外（郭长宝等，2024），滑坡体入水产生的涌浪也会影响周边区域.滑坡成坝形成的堰塞坝会堵塞河道，随着上游水位的不断上升，发生堰塞坝溃决（贾珂程等，2023），更会对下游居民生命财产安全带来严重威胁（石振明等， 2023）.如金沙江流域的白格地区，在2018年10月和11月，连续发生了两次滑坡堵江，分别都造成堵塞和堰塞坝溃决，分别造成了超20 000人和近75 000人疏散（Cai et al.， 2020）.最大滑坡堵江预期库容近8×10⁸ m³，人工开挖泄流槽后，仍有不可避免的洪水导致下游金沙江大桥7 跨桥面被完全冲毁，部分道路、学校等被淹，为金沙江下游沿岸带来巨大经济损失.随着地震活动和极端降雨频发，滑坡堵江事件呈增长趋势（吴昊等，2023），深入理解其形成演进机制对保障山区人民生命财产安全有重要意义.

目前一系列的现场调查（许强等，2018； Ouyang et al.， 2019）和室内试验（Luo et al.， 2025）加深了对该灾害链的理解，但是受到成本和观测手段的限制，难以大范围开展相关研究.随着计算机硬件及相关算法的发展，数值模拟方法作为一种低成本的方法，逐渐被应用于滑坡堵江灾害链的研究中（周礼等，2019；Zhou et al.， 2020； 杜文杰等，2022；李东阳等，2023）.然而，目前的连续性研究多基于网格法的研究（Long et al.， 2024），如有限差分法（FDM），有限体积法（FVM）等，在模拟滑坡体大变样方面会面临网络畸变等问题（杜文杰等，2022），其中广泛采用的深度平均算法忽略了材料在高度上的差异及复杂的三维形态，影响模拟结果的精度.在非连续法方面，离散元法（DEM）是模拟滑坡涌浪过程最常用的方法之一.将有限元法（FEM）与DEM结合形成的FDEM方法，能够更精细地模拟岩体破裂行为.DEM和FDEM可与流体模拟方法（如计算流体力学CFD和光滑粒子流体力学SPH）耦合，实现水土耦合模拟，从而再现滑坡涌浪过程，例如CFD-DEM（Nian et al.， 2021；Li et al.， 2023a）、DEM-SPH（Li et al.， 2024a）以及FDEM-SPH（Bao et al.， 2023）等方法均在室内试验尺度和现场尺度滑坡涌浪进行了应用.但是DEM及FDEM计算过程需要计算每个粒子间的接触作用，相比于连续性方法计算成本高昂（Li et al.， 2024b），同时需要标定的参数较多（Huang et al.， 2023），限制了其在现场尺度模拟的广泛应用.尽管效率更高的非连续变形分析法（DDA）也被用来模拟滑坡运动，进而开发DDA-SPH等算法模拟滑坡涌浪（Peng et al.， 2020），但是不同算法之间的耦合会带来不可避免的信息传递误差，影响计算精度.因此，统一平台的、高效率的、可直接考虑滑坡体大变形的水土耦合方法有待被开发应用于滑坡堵江三维模拟中.

光滑粒子流体动力学方法（SPH）是一种无网格连续性方法，由一系列SPH粒子携带物理信息，并且可以自由运动，兼顾了能够模拟滑坡体大变形和计算效率的优势，逐渐被应用于滑坡灾害的模拟中（Bui et al.， 2008； Bui and Nguyen， 2021； Li et al.， 2023b）.此外，不同的对象（如水、土）可以直接采用不同的SPH粒子进行模拟，并且不同粒子的交界面可以根据邻近关系自动识别（Peng et al.， 2024），近些年逐渐有相关的水土耦合算法被开发和应用（Bui and Nguyen， 2021； Li et al.， 2024b）.基于水土耦合方式，相关的模型可以分为单向耦合模型和双向耦合模型两类（Zhu et al.， 2022）.单相模型通常隐式计算土体，认为土体是不可运动的孔介质，根据渗流理论模拟土体对水流的阻碍作用.相关的模型对孔隙率动态更新、渗流力计算及本构方程等方面进行优化，被应用于坝体渗流（Feng et al.， 2022）、潜堤消浪（Ren et al.， 2016）等案例中，但是这些方法忽略了土体的运动，无法应用于滑坡堵江的模拟.双向耦合模型则是在此基础上，认为土体可以自由大变形，基于水土两相相互作用，实现完全耦合.基于本构方程改写和渗流力作用模拟相互作用，被应用于射流冲击（Bui and Nguyen， 2021）、海底滑坡（Zhang et al.， 2024a）等案例的模拟中，在滑坡-堵江-成坝灾害链有很好的应用前景.

在双向耦合方法中，如何定义水土混合模型的控制方程是水土耦合的核心.其中混合物理论是最常用的理论模型（Bui and Nguyen， 2021； Zhu et al.， 2022），该方法认为水土和土体共同充满整个计算空间，在运动过程中，两者体积分数发生变化，从而影响控制方程的计算.Zhu et al.（2022）对混合物理论模型的相关方法进行了整理，对比了本构方程的区别，并改进了现有方程，通过一系列室内试验案例证明所改进的方程在理论上和数值结果上有更好的体积守恒性和更准确的结果.

本文基于水土两向耦合SPH算法（Zhu et al.， 2022），将其应用于滑坡-堵江-成坝模拟中，基于室内滑坡入水试验进行了验证.借助现场调查和数据收集，该算法被进一步应用到2018年“10·11”金沙江特大白格滑坡堵江事件（许强等，2018），基于数值模拟结果对该灾害链的动力学过程进行分析，并讨论不同内摩擦角下运动过程的区别及其对成坝形态的影响.本文旨在开发的滑坡-堵江-成坝灾害链的全过程数值模拟方法，并加深对其动力学过程的理解，为山地地区相关灾害的防灾减灾提供支撑.

1 SPH水土耦合方法

1.1　控制方程

在本文，主要有SPH水粒子和SPH土粒子来分别模拟水体和滑坡体，如图1所示.基于混合物理论，水土和土体充满计算的物理空间，两者的体积分数满足以下条件：

{\varphi }_{\mathrm{f}}+{\varphi }_{\mathrm{s}}=\mathrm{1},(1)

式中：{\varphi }_{\mathrm{f}}和{\varphi }_{\mathrm{s}}分别为流体和土体的体积分数，在拉格朗日形式下基于材料密度求解两相质量守恒方程（Zhu et al.， 2022），从而得到：

\frac{{\mathrm{d}}^s{\rho }_{\mathrm{s}}}{\mathrm{d}t}=-{\rho }_{\mathrm{s}}\nabla \cdot {\mathit{u}}_{\mathrm{s}}，(2)

\frac{{\mathrm{d}}^f{\rho }_{\mathrm{f}}}{\mathrm{d}t}=-\frac{{\rho }_{\mathrm{f}}}{{\varphi }_{\mathrm{f}}}\nabla \cdot \left({\varphi }_{\mathrm{f}}{\mathit{u}}_{\mathrm{f}}+{\varphi }_{\mathrm{s}}{\mathit{u}}_{\mathrm{s}}\right)，(3)

式中：{\rho }_{\mathrm{s}}和{\rho }_{\mathrm{f}}分别为土体和流体的材料密度.{\mathit{u}}_{\mathrm{s}}和{\mathit{u}}_{\mathrm{f}}分别为土体和流体的速度向量.基于线性动量守恒，水体和土体的动量守恒方程分别为：

\frac{{\mathrm{d}}^s{\mathit{u}}_{\mathrm{s}}}{\mathrm{d}t}=\frac{\nabla \cdot {\mathbf{\sigma }}_{\mathrm{s}}}{{\rho }_{\mathrm{s}}}+\mathbf{g}-\frac{{\varphi }_{\mathrm{s}}\nabla p_{\mathrm{f}}}{{\rho }_{\mathrm{s}}}-\frac{{\mathit{f}}_{\mathrm{d}}}{{\rho }_{\mathrm{s}}}，(4)

\frac{{\mathrm{d}}^f{\mathit{u}}_{\mathrm{f}}}{\mathrm{d}t}=-\frac{\nabla p_{\mathrm{f}}}{{\rho }_{\mathrm{f}}}+\mathbf{g}+\frac{{\mathit{f}}_{\mathrm{d}}+{\mathit{f}}_{\mathrm{v}}}{{\varphi }_{\mathrm{f}}{\rho }_{\mathrm{f}}}，(5)

式中：g=9.8 m/s²为重力加速度， σ_s为土体的柯西应力张量， f_d和 f_v分别为渗流力和粘性力，其公式为：

\begin{array}{l}{\mathit{f}}_{\mathrm{d}}=-\mathrm{150}\frac{\mu (\mathrm{1}-{\varphi }_{\mathrm{f}}{)}^{\mathrm{2}}}{{\varphi }_{\mathrm{f}}{d}_{\mathrm{50}}^{\mathrm{2}}}({\mathit{u}}_{\mathrm{f}}-{\mathit{u}}_{\mathrm{s}})\\ -\mathrm{1.75}\frac{{\rho }_{\mathrm{f}}(\mathrm{1}-{\varphi }_{\mathrm{f}})}{d_{\mathrm{50}}}‖{\mathit{u}}_{\mathrm{f}}-{\mathit{u}}_{\mathrm{s}}‖({\mathit{u}}_{\mathrm{f}}-{\mathit{u}}_{\mathrm{s}})\end{array}，(6)

{\mathit{f}}_{\mathrm{v}}=\eta {\nabla }^{\mathrm{2}}{\mathit{u}}_{\mathrm{f}}，(7)

式中：\eta为动力粘度.在SPH计算中，水土都会被离散为不同的SPH粒子，如图1所示.采用SPH核函数W_{ij}=W(r_i-r_j,h_{\mathrm{s}})可以对物理方程进行离散化，从而对每个SPH粒子进行计算，其中i，j为粒子的序号，h_s为光滑半径.对流体的质量守恒和动量守恒方程进行离散化后，可以得到流体的SPH形式控制方程为：

\frac{{\mathrm{d}}^f{\rho }_i}{\mathrm{d}t}=-\frac{{\rho }_i}{{\varphi }_i}\sum _j({\varphi }_j{\mathit{u}}_j-{\varphi }_i{\mathit{u}}_i)\cdot {\nabla }_i{W}_{ij}{V}_j-\frac{{\rho }_i}{{\varphi }_i}\sum _j(\mathrm{1}-{\varphi }_j){\mathit{u}}_{\mathrm{j}}^s\cdot {\nabla }_i{W}_{ij}{V}_j+\widehat{\rho }{\nabla }_i{W}_{ij}V，（8）

\frac{{\mathrm{d}}^f{\mathit{u}}_i}{\mathrm{d}t}=-\frac{\mathrm{1}}{{\rho }_i}\sum _j(p_i+p_j){\nabla }_i{W}_{ij}{V}_j+\frac{\mathrm{1}}{{\varphi }_{\mathrm{i}}{\rho }_i}\sum _j\mathrm{2}\eta \frac{{\mathbf{v}}_{ij}}{r_{ij}}\frac{\partial W}{\partial r_{ij}}{V}_j-\frac{\mathrm{1}}{{\varphi }_{\mathrm{i}}{\rho }_i}{\displaystyle \sum _a^{}}{f}_{ia}{W}_{ia}{V}_a+\widehat{P}{\nabla }_i{W}_{ij}{V}_j+{\mathbf{g}}_i，（9）

式中：i，j表示水粒子的序号，a，b表示土粒子的序号，V为粒子体积，{\mathit{u}}_{ij}={\mathit{u}}_i-{\mathit{u}}_j是流体粒子的相对速度.

基于黎曼相互作用求解（Zhang et al.， 2017），\widehat{\rho }和\widehat{P}分别为：

\widehat{\rho }={\rho }_{\mathrm{L}}\frac{P_{\mathrm{L}}-P_{\mathrm{R}}}{{\rho }_{\mathrm{L}}{c}_{\mathrm{L}}+{\rho }_{\mathrm{R}}{c}_{\mathrm{R}}}，(10)

\widehat{P}=\frac{\mathrm{2}{\rho }_{\mathrm{L}}{c}_{\mathrm{L}}{\rho }_{\mathrm{R}}{c}_{\mathrm{R}}\left(U_{\mathrm{L}}-U_{\mathrm{R}}\right)\beta }{{\rho }_{\mathrm{L}}({\rho }_{\mathrm{L}}{c}_{\mathrm{L}}+{\rho }_{\mathrm{R}}{c}_{\mathrm{R}})}，(11)

式中：\beta =\mathrm{m}\mathrm{i}\mathrm{n}\left(\xi \mathrm{m}\mathrm{a}\mathrm{x}\left(\left(U_{\mathrm{L}}-U_R\right)/\tilde{c},\mathrm{0}\right),\mathrm{1}\right)，其中耗散系数ζ=3.0，\tilde{c}=\left({\rho }_{\mathrm{L}}{c}_{\mathrm{L}}+{\rho }_{\mathrm{R}}{c}_{\mathrm{R}}\right)/\left({\rho }_{\mathrm{L}}+{\rho }_{\mathrm{R}}\right).下标L，R分别代表黎曼相互作用中的左右项，计算方式为：

\left\{\begin{array}{c}({\rho }_{\mathrm{L}}^{},U_{\mathrm{L}}^{},P_{\mathrm{L}}^{},c_{\mathrm{L}}^{})=({\rho }_i^{},-v_i^{}\cdot {\mathrm{e}}_{ij}^{},p_i^{},c_i^{})\\ ({\rho }_{\mathrm{R}}^{},U_{\mathrm{R}}^{},P_{\mathrm{R}}^{},c_{\mathrm{R}}^{})=({\rho }_j^{},-v_j^{}\cdot {\mathrm{e}}_{ij}^{},p_j^{},c_j^{})\end{array}\right.，(12)

式中：c为声速.采用同样的方式，可以对土体的质量守恒方程和能量守恒方程进行黎曼SPH离散化，得到：

\frac{{\mathrm{d}}^s{\rho }_{\mathrm{a}}}{\mathrm{d}t}={\rho }_{\mathrm{a}}\sum _{\mathrm{b}}(u_{\mathrm{a}}-u_{\mathrm{b}}){\nabla }_{\mathrm{a}}{W}_{\mathrm{a}\mathrm{b}}{V}_{\mathrm{b}}+\widehat{\rho }{\nabla }_i{W}_{ij}V，(13)

\begin{array}{l}\frac{{\mathrm{d}}^s{u}_{\mathrm{a}}}{\mathrm{d}t}=\frac{\mathrm{1}}{{\rho }_{\mathrm{a}}}\sum _{\mathrm{b}}({\sigma }_{\mathrm{a}}+{\sigma }_{\mathrm{b}}){\nabla }_{\mathrm{a}}{W}_{\mathrm{a}\mathrm{b}}{V}_{\mathrm{b}}+{\mathrm{g}}_a-\\ \frac{{\varphi }_{\mathrm{a}}}{{\rho }_{\mathrm{a}}}\sum _i(p_i-p_a){\nabla }_{\mathrm{a}}{W}_{\mathrm{a}i}{V}_i+\\ \frac{\mathrm{1}}{{\rho }_{\mathrm{a}}}{\displaystyle \sum _i^{}}{f}_{\mathrm{a}i}{W}_{\mathrm{a}i}{V}_i+\widehat{P}{\nabla }_i{W}_{ij}V\end{array}，(14)

式中：\widehat{\rho }和\widehat{P}计算方式与水体相同，但根据Zhang et al.（2024b）的研究，在土体中耗散系数ζ取20D， 其中D代表维度，在二维算例下取2，三维算例下取3.

1.2　本构方程

流体基于弱可压缩牛顿流体模型进行求解，其本构方程为采用状态方程：

p=c_{\mathrm{0}}^{\mathrm{2}}(\rho -{\rho }_{\mathrm{0}}^{}),(15)

式中：c₀为数值声速，ρ₀为初始密度，为保证流体弱可压缩条件，取c₀=10U_max，其中U_max为期望最大速度.

土体采用Drucker-Prager（DP）本构的弹塑性体模型，其屈服条件为：

f(I_{\mathrm{1}},J_{\mathrm{2}})={\alpha }_{\varphi }{I}_{\mathrm{1}}+\sqrt[]{J_{\mathrm{2}}}-k_c，(16)

式中：I_{\mathrm{1}}=\mathrm{t}\mathrm{r}(\mathit{\sigma })是第一应力不变量，J_{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}\mathit{s}:\mathit{s}是第二偏应力张量，{\alpha }_{\varphi }和k_c是DP模型参数，其定义为：

{\alpha }_{\varphi }=\frac{\mathrm{t}\mathrm{a}\mathrm{n}\varphi }{\sqrt[]{\mathrm{9}+\mathrm{12}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\varphi }},k_c=\frac{\mathrm{3}c}{\sqrt[]{\mathrm{9}+\mathrm{12}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\varphi }}，(17)

式中：\varphi和c分别代表土体的内摩擦角和粘聚力.土体的应力采用以下方式进行积分：

\mathit{\sigma }(t+\Delta t)=\mathit{\sigma }(t)+\dot{\mathit{\sigma }}\Delta t，(18)

式中：∆t是时间步，\dot{\mathit{\sigma }}是土体的应力\mathit{\sigma }的变化率，相同的标注方式在后文也同样表示物理量的变化率.采用非关联流动法则计算应力率：

\begin{array}{c}\dot{\mathit{\sigma }}\end{array}=\mathrm{2}G{\dot{\mathit{\epsilon }}}^s+Ktr\left(\dot{\mathit{\epsilon }}\right)\mathit{I}-\left[\dot{\lambda }\mathrm{3}{\alpha }_{\psi }K\mathit{I}+\frac{G}{\sqrt[]{J_{\mathrm{2}}}}{\mathit{\sigma }}^{\mathrm{s}}\right]+\mathit{\sigma }•{\dot{\mathit{\omega }}}^T+\dot{\mathit{\omega }}•\mathit{\sigma }，(19)

式中：K和G分别为体积模量和剪切模量，\dot{\mathbf{\epsilon }}是应变张量，{\dot{\mathit{\epsilon }}}^{\mathrm{s}}=\dot{\mathit{\epsilon }}-\mathrm{t}\mathrm{r}\left(\dot{\mathit{\epsilon }}\right)\mathit{I}/\mathrm{3}是偏应变率张量，\dot{\mathit{\omega }}是旋转张量，λ是塑性应变乘子，它们的计算公式为：

\dot{\lambda }=\frac{\mathrm{3}{\alpha }_{\varphi }Ktr(\dot{\mathit{\epsilon }})+(G/\sqrt[]{J_{\mathrm{2}}}){\mathbf{\sigma }}^s:\dot{\mathit{\epsilon }}}{\mathrm{9}{\alpha }_{\varphi }K{\alpha }_{\psi }+G}，(20)

\dot{\mathit{\omega }}=\frac{\mathrm{1}}{\mathrm{2}}\left[\nabla {\mathit{v}}^{\mathrm{s}}-(\nabla {\mathit{v}}^{\mathrm{s}}{)}^{\mathrm{T}}\right]，(21)

\dot{\mathit{\epsilon }}=\frac{\mathrm{1}}{\mathrm{2}}\left[\nabla {\mathit{v}}^{\mathrm{s}}+(\nabla {\mathit{v}}^{\mathrm{s}}{)}^{\mathrm{T}}\right].(22)

此外，应力拉回算法（Bui et al.， 2008）被采用，以保证应力状态不超过屈服面，具体为： ，（23）

式中：\tilde{\mathbf{\sigma }}是拉回后的应力张量.

本文基于开源C++代码库SPHinXsys （Zhang et al.， 2021）自主开发了上述代码，并基于英特尔至强Platinum 8377C处理器进行并行加速计算.在所有的SPH模拟中均采用Wendland C2核函数（Wendland， 1995），并且设置光滑半径h_{\mathrm{s}}=\mathrm{1.3}dp，其中dp为初始粒子间距.

2 数值验证

经典的滑坡入水涌浪室内试验（Viroulet et al.， 2013，2014）被用来验证本文所采用的水土耦合SPH算法准确性.图2展示了该模型的初始设置，一个体积为2 kg的散粒体被放置在坡角为45°的斜坡上，在一个挡板的限制下保持静止，以模拟天然滑坡体.下方设置高度为0.15 m的静止水体，以模拟天然河流.突然撤下挡板后，滑坡体会在自重作用下发生失稳大变形，从而落入水中，堆积停止并产生涌浪.在距离滑坡体不同距离处设置了共4个浪高仪W1~W4用以监测水位变化， 距离滑坡体的位置分别为0.45 m，0.75 m，1.05 m和1.35 m.

SPH模型按照试验进行设置，如图2所示，其中墙体用墙粒子表示，不会发生运动，流体和土体均采用无滑移边界（Bui et al.， 2008； Adami et al.， 2012）.水体和土体的参数基于室内试验进行设置，SPH参数参考Huang et al.（2023）的研究设置.具体如表1所示.

图3展示了室内试验（Viroulet et al.， 2013，2014）所拍摄的不同时刻照片.在t=0 s时，滑坡体挡板被快速撤掉，此时滑坡体开始失稳变形.在t=0.3 s时，滑坡体前缘进入水中，与水体发生剧烈混合和冲击，造成局部的涌浪.在t=0.9 s时，滑坡体完全落入水中，并且在冲击处产生了二次涌浪.滑坡体逐渐趋于稳定，所产生的波浪逐渐向远处传播，如t=1.5 s时所示.在图3的右侧展示了不同时刻下SPH模拟得到的结果，包括滑坡运动、滑坡入水、部其滑坡体运动形式及涌浪的形态与室内试验基本一致，证明了所采用的数值模型在滑坡体与水相互作用时的准确性.

图4进一步展示了不同监测点的浪高曲线，滑坡所导致的涌浪最早于0.40 s传播到W1，并且在约0.26 s后到达峰值，最大峰值浪高为2.1 cm.随后首波的波峰分别于0.95 s，1.18 s和1.44 s传递到浪高仪W2，W3及W4上.从图3中的结果可以得到，这些首波主要是由于滑坡体前缘入水导致的.而滑坡体后缘入水后，同样产生了涌浪，但是其幅值远小于首波，如图4a~4b所示.随着首波传播到尾端并发生反射，与次波发生了叠加，其幅值在W4甚至超过了首波.SPH模拟的结果也同样被展示在图中，可以看出所模拟得到波形与室内试验结果一致，说明数值模拟可以复现首波，次波以及首波的反射等现象，证明了该模型的准确性.

3 白格滑坡模拟与分析

3.1　白格概况及数值模型设置

在我国西藏地区的白格地区， 分别于2018年的10月和11月发生了两次巨型天然滑坡.如前文介绍，此次滑坡及次生灾害严重毁坏了下游流域，并带来了巨大的财产损失，本文针对于占主要体积的”10·11”白格滑坡体进行了SPH全过程模拟与分析.图5a展示了白格滑坡的位置，它位于金沙江流域的右岸（31° 4′56.41″ N， 98° 42′ 17.98″ E），下游有巴塘、丽江等多个县市.图5b进一步展示了白格滑坡堵江后的谷歌地球影像， 滑坡体方量超2.3×10⁷ m³，其滑动区宽度范围为500~700 m，巨大的滑坡堵塞了金沙江，形成了白格堰塞坝.

本文基于高精度现场数字高程模型（Bao et al.， 2023）及相关现场调查和影像（许强等，2018），建立了白格地区的地形、滑坡体和水体模型.其中滑坡体的体积为2.316×10⁷ m³，初始水位高度设置为25 m（Bao et al.， 2023），河流的上下游人工设置刚性边界进行限制.其他的SPH模型参数基于已有现场调查及数值模拟标定结果（Zhang et al.， 2019； Fan et al.， 2020； Li et al.， 2023b），具体如表2所示.

图6展示了SPH初始模型的设置，除上文介绍的形态和参数设置外，一系列监测点被设置在不同位置，以分析白格滑坡体滑动堵江过程中的动力学过程，可以得到3类监测信息：（1）滑坡体速度监测.设置了4个不同高度的速度监测点，分别位于z=195 m， 290 m， 555 m和750 m，以监测在滑坡体前缘，中部及后缘不同位置的滑坡速度演化.（2）浪高监测点，在预期滑坡堆积区的上下游布设了不同的浪高仪，其具体坐标列在图6，用以监测不同位置下水位变化，以定量分析滑坡涌浪的形成.（3）能量监测.通过计算每个粒子动能并求和，对滑坡体整体的动能演化进行了监测和记录.

3.2　白格灾害链全过程模拟

图7展示了SPH模拟得到不同时刻下的白格滑坡体运动-堵江-成坝的全过程，同时绘制了滑坡体速度和流体速度的变化.在t=10 s时，滑坡体在重力作用下开始运动，表面速度最高约40 m/s；在t=20 s时，滑坡体加速运动，其表面最高速度近80 m/s，同时滑坡体前缘进入水中，在局部与水体发生相互作用，产生局部的涌浪.在t=40 s时，滑坡体前缘已完全冲击进入河道，产生剧烈的涌浪并开始向远处传播；在t=60 s时，大部分滑坡体堆积在河道上，形成河道的堵塞，同时涌浪继续向远处传播.最后，滑坡体逐渐停止运动，形成堰塞坝堆积体，同时波浪也随着传播距离的变远逐渐减弱，滑坡堵江成坝过程完成，如图7e~7f所示.

图8进一步展示了SPH模拟得到的白格滑坡堆积形成的堰塞坝，及其与现场调查结果（Ouyang et al.， 2019； Li et al.， 2023a）的对比.其中堆积区域指形成堰塞坝体的区域，影响区域包括滑坡运动区域及涌浪冲击区域.同时，CFD-DEM模型（Li et al.， 2023a）和深度平均模型（Ouyang et al.， 2019）的模拟结果也被绘制在图上.对于不同算法结果，深度平均算法由于无法充分获得滑坡体三维形态，模拟结果缺少堆积体四周堆积高度较低的区域，其堆积区域略小于CFD-DEM与SPH的模拟结果.CFD-DEM的模拟结果和本文SPH模拟结果基本一致，其与现场调查结果相比，下游堆积区形态一致，上游堆积区略大于现场调查所得堆积区.这一差异可能是由于模拟没有考虑部分上游松散滑坡体被流体带走的现象，但模拟结果都处于现场调查所得灾害影响区域范围内.此外，模拟结果显示，堰塞坝的高度变化在河道方向上呈现中间高、上下游低的特征，最高高度约80 m，这也与目前的现场调查结果一致（Fan et al.， 2020）.

3.3　灾害链全过程分析

图9展示了不同监测点下，监测得到的滑坡体速度随时间的演化.对于S1，从滑坡启动开始速度逐渐增加，于t=11 s达到峰值后迅速减小至0，说明此时滑坡体后缘已经低于S1，如图7b所示.对于S2，其形式与S1类似，于t=17 s时达到峰值速度32 m/s后迅速下降.但是由于在S2附近有少量滑坡体残留和持续运动，使其监测到的速度出现多次波动.对于滑坡体前缘的监测点S3和S4，均在20~30 s间达到峰值后衰减，此时峰值区域对应滑坡体入水发生冲击阶段，如图7b~7c所示.图10展示了涌浪过程中最高水面的变化及不同测点下W1~W4的浪高数据，可以看出，滑坡体推移水体冲击对岸，如图7c所示，其最大浪高可达125 m，与现场调查所得到大于100 m一致（许强等，2018），并且与FDEM-SPH模拟得到的120 m相近（Bao et al.， 2023）.对于不同测点的涌浪高度，距离冲击点最近的W3在t=23 s后观测到了涌浪，证明了该阶段的剧烈冲击涌浪现象.而在发生冲击后，滑坡体的速度增长速率变缓出现速度转折，缓慢达到峰值后慢慢降低.滑坡运动过程中的动能演化同样被绘制在图中，在滑坡启动滑移后，滑坡体动能快速增加，于23 s时达到峰值，最大动能超过了3.5×10¹³ J.由于河道和水体的阻碍，滑坡体堆积成坝，动能迅速减小，在约60 s后基本降为0.河道相比于该巨型滑坡体，其尺寸较小，滑坡体的高速冲击带来了巨大的涌浪，如图10所示.在距离滑坡体较近的监测点W3，最大涌浪高度达40 m.即使是距离堆积体中心约1 000 m的W1，其涌浪高度也可达23 m，与其初始水位25 m相当，体现了该滑坡对相关流域的重大影响.

结合滑坡体速度与能量的演化，以及滑坡涌浪的观测，可以将滑坡堵江过程分为3个主要的阶段：（1）滑坡运移阶段——从初始时刻到前缘速度发生转折前，滑坡体在山体上运移，滑坡体速度快速上升，重力势能快速转变为动能.（2）涌浪阶段——从前缘速度转折到动能归零之间，此时滑坡体冲击河道水体，与水体发生相互作用，在河道中产生涌浪.（3）堵江成坝阶段——动能归零后的阶段.此时滑坡体已经基本停止，仅有少量的坡面残余体继续运动，形成堵塞河道的堰塞坝.

3.4　滑坡体内摩擦对堰塞坝堆积形态的影响

滑坡体的强度参数（内摩擦角和粘聚力）是影响滑坡体运动和堆积特征的重要因素，并且现有研究证明内摩擦角φ在其中占主控因素（Li et al.， 2023b），但是目前关于φ对堰塞坝三维堆积形态的讨论相对较少.为此，本文模拟了一系列不同内摩擦角φ分别为5°，10°，15°，20°下白格滑坡堆积成坝的过程.图11a~11d展示了不同内摩擦角下的最终滑坡堆积体形态.随着φ从5°增长到20°，滑坡体的强度增加，流动性降低，运动放缓.

图11e进一步展示了不同φ下滑坡体动能随时间的演化.可以看出，不同φ下滑坡体的动能呈类似的变化趋势，说明其保持着与前文介绍的类似运动形式.随着φ的增加，整体动能的大小显著降低，这使得滑坡体的长度持续降低，从1 818 m最终降低至824 m，使堰塞坝的堵塞区域显著降低.相对应地，随着φ从5°增长到15°，滑坡体高度从48 m持续增长至96 m.但是当φ增长至20°时，由于滑坡体的运动性过差，过多的滑坡残余堆积在山体上，使得滑坡体没有完全堵塞河道，其形态与低摩擦角下出现明显差异.相对应地，由于更少体积的滑坡体成功堆积在河道上，其高度也基本上停止增长，在φ=20°时高度仅为98 m.堰塞坝形态与内摩擦角φ之间的定量关系被绘制在图11f中，并且拟合了l和h随φ的变化曲线.可以看出，随着φ的增加，l呈现近似直线的下降趋势，拟合直线为：l=-\mathrm{3}\mathrm{674}\mathrm{t}\mathrm{a}\mathrm{n}\left(\phi \right)+\mathrm{2}\mathrm{103}，而h则呈现逐渐放缓的增长趋势，其关系拟合曲线为h=-\mathrm{7.5}\mathrm{t}\mathrm{a}\mathrm{n}{\left(\phi \right)}^{-\mathrm{0.92}}+\mathrm{118}，呈幂函数形式的非线性特征.

3.5　滑坡体内摩擦对涌浪特征的影响

除堆积形态外，滑坡体内摩擦角φ通过影响滑坡体运动特征，对滑坡涌浪特征也有着重要的影响.图12a~12d展示了不同内摩擦角的滑坡冲击入水时的形态，选自不同内摩擦角下的测点W3的浪高峰值时刻，如图12e所示.如前文所述，由于内摩擦角的增加导致滑坡体动能耗散增加，其入水速度显著降低，使得水土之间的相互作用力显著减小，进而减弱涌浪.图12e展示了不同φ时W3的浪高变化曲线，随着φ的增加，滑坡导致的涌浪过程明显延后并且涌浪高度降低，其峰值时刻由φ=5°时的28 s延后到φ=20°时的37 s，其峰值浪高也从最高的43.2 m降低至36.9 m，体现了滑坡流动性降低对涌浪的显著减弱效果.

为进一步定量分析涌浪高度与滑坡运动特征关系，图12f展示了入水时的峰值浪高H_max与滑坡入水时弗劳德数\mathrm{F}\mathrm{r}=\sqrt[]{v_{\mathrm{i}\mathrm{n}}/\mathrm{g}H}的关系，其中v_in是入水前的滑坡前缘速度，即入水速度，g是重力加速度，H为初始水面高度.弗劳德数Fr体现了滑坡体惯性力与重力的比值，是衡量滑坡体入水动能的重要参数（Heller and Hager， 2011）.是随着Fr的增加，滑坡体入水动能增加，涌浪峰值高度呈近似线性增加趋势.线性拟合得到公式为H_max=21.18Fr-2.5，其决定系数R²大于0.97.该公式拟合得到的系数21.18体现了浪高对于滑坡体Fr的高敏感性，即滑坡体入水动能的增加显著影响涌浪强度，同时其定量拟合关系有助于提高对滑移涌浪强度的预测.

4 结论

本文采用SPH水土耦合算法对滑坡-堵江-成坝全过程灾害链进行模拟，并进一步模拟和分析了金沙江白格“10·11”滑坡堵江事件.具体结论如下：

（1）所采用的双层水土双向耦合SPH算法可以准确复现滑坡体入水过程，对滑坡灾害链模拟有重要应用前景.

（2）对金沙江白格“10·11”特大滑坡堵江案例进行了模拟，基于速度突变和能量归零时刻及涌浪特征，可以准确划分该过程的滑坡运移，入水涌浪和滑坡堵江3个阶段.

（3）随着内摩擦角的增加，滑坡运移的动能明显减小，使得滑坡堆积堰塞坝长度显著降低，同时高度变高，在摩擦角φ到达20°时无法完全堵江.定量研究表示，滑坡堆积堰塞坝的长度随tan（φ）的增加呈近似线性减小，高度随tan（φ）呈幂函数形式增加，最终趋于平缓.

（4）随着内摩擦角的增加，由于其动能降低，滑坡入水时产生的涌浪强度显著降低，并且峰值时刻延后.峰值浪高H_max与入水时的滑坡前缘弗劳德数Fr之间的定量关系显示H_max对Fr的变化敏感，并且随着Fr的增加呈近似线性增加趋势.

参考文献

[1] Adami, S., Hu, X. Y., Adams, N. A., 2012. A Generalized Wall Boundary Condition for Smoothed Particle Hydrodynamics. Journal of Computational Physics, 231(21): 7057-7075. https://doi.org/10.1016/j.jcp.2012.05.005

[2] Bao, Y. D., Su, L. J., Chen, J. P., et al., 2023. Dynamic Process of a High-Level Landslide Blocking River Event in a Deep Valley Area Based on FDEM-SPH Coupling Approach. Engineering Geology, 319: 107108. https://doi.org/10.1016/j.enggeo.2023.107108

[3] Bui, H. H., Fukagawa, R., Sako, K., et al., 2008. Lagrangian Meshfree Particles Method (SPH) for Large Deformation and Failure Flows of Geomaterial Using Elastic-Plastic Soil Constitutive Model. International Journal for Numerical and Analytical Methods in Geomechanics, 32(12): 1537-1570. https://doi.org/10.1002/nag.688

[4] Bui, H. H., Nguyen, G. D., 2021. Smoothed Particle Hydrodynamics (SPH) and Its Applications in Geomechanics: From Solid Fracture to Granular Behaviour and Multiphase Flows in Porous Media. Computers and Geotechnics, 138: 104315. https://doi.org/10.1016/j.compgeo.2021.104315

[5] Cai, Y. J., Cheng, H. Y., Wu, S. F., et al., 2020. Breaches of the Baige Barrier Lake: Emergency Response and Dam Breach Flood. Science China Technological Sciences, 63(7): 1164-1176. https://doi.org/10.1007/s11431-019-1475-y

[6] Du, W. J., Sheng, Q., Yang, X. H., et al., 2022. Chain Generation Process of Landslide Blocking River Based on Two-Phase Double-Point Material Point Method. Advanced Engineering Sciences, 54(3): 36-45 (in Chinese with English abstract).

[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. 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

[9] Feng, D. L., Neuweiler, I., Huang, Y., 2022. Numerical Modeling of Wave-Porous Structure Interaction Process with an SPH Model. Scientia Sinica Physica, Mechanica & Astronomica, 52(10): 104715. https://doi.org/10.1360/sspma-2022-0216

[10] Guo, C. B., Wu, R. A., Zhong, N., et al., 2024. Large Landslides along Active Tectonic Zones of Eastern Tibetan Plateau: Background and Mechanism of Landslide Formation. Earth Science, 49(12): 4635-4658 (in Chinese with English abstract).

[11] Heller, V., Hager, W. H., 2011. Wave Types of Landslide Generated Impulse Waves. Ocean Engineering, 38(4): 630-640. https://doi.org/10.1016/j.oceaneng.2010.12.010

[12] Huang, C., Hu, C., An, Y., et al., 2023. Numerical Simulation of the Large-Scale Huangtian (China) Landslide-Generated Impulse Waves by a GPU-Accelerated Three-Dimensional Soil-Water Coupled SPH Model. Water Resources Research, 59(6): e2022WR034157. https://doi.org/10.1029/2022wr034157

[13] 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).

[14] Li, D. Y., Nian, T. K., Tiong, R. L. K., et al., 2023a. River Blockage and Impulse Wave Evolution of the Baige Landslide in October 2018: Insights from Coupled DEM-CFD Analyses. Engineering Geology, 321: 107169. https://doi.org/10.1016/j.enggeo.2023.107169

[15] Li, D. Y., Nian, T. K., Wu, H., et al., 2023. Coupled DEM–CFD Method for Landslide-River Blockage-Impulse Wave Disaster Chain Simulation and Its Application. Advanced Engineering Sciences, 55(1): 141-149 (in Chinese with English abstract).

[16] Li, S., Peng, M., Gao, L., et al., 2024a. 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

[17] Li, S., Tang, H., Peng, C., et al.,2023b. Sensitivity and Calibration of Three-Dimensional SPH Formulations in Large-Scale Landslide Modeling. Journal of Geophysical Research: Solid Earth, 128(2): e2022JB024583. https://doi.org/10.1029/2022jb024583

[18] Li, Y., Liu, H. Q., Yang, L., et al.,2024b. An Optimized DEM-SPH Model for Surge Waves Induced by Riverside Landslides. International Journal for Numerical and Analytical Methods in Geomechanics, 48(1): 270-286. https://doi.org/10.1002/nag.3638

[19] Long, X. Y., Hu, Y. X., Gan, B. R., et al., 2024. Numerical Simulation of the Mass Movement Process of the 2018 Sedongpu Glacial Debris Flow by Using the Fluid-Solid Coupling Method. Journal of Earth Science, 35(2): 583-596. https://doi.org/10.1007/s12583-022-1625-1

[20] Luo, H. W., Zhou, G. G. D., Lu, X. Q., et al., 2025. Experimental Investigation on the Formation and Failure of Landslide Dams Considering the Landslide Mobility and River Flow. Engineering Geology, 346: 107873. https://doi.org/10.1016/j.enggeo.2024.107873

[21] Nian, T. K., Wu, H., Takara, K., et al., 2021. Numerical Investigation on the Evolution of Landslide-Induced River Blocking Using Coupled DEM-CFD. Computers and Geotechnics, 134: 104101. https://doi.org/10.1016/j.compgeo.2021.104101

[22] Ouyang, C. J., An, H. C., Zhou, S., et al., 2019. Insights from the Failure and Dynamic Characteristics of Two Sequential Landslides at Baige Village along the Jinsha River, China. Landslides, 16(7): 1397-1414. https://doi.org/10.1007/s10346-019-01177-9

[23] Peng, M., Li, S., Gao, L., et al., 2024. A Novel Local-Drag-Force-Based Approach for Simulating Wave Attenuation by Mangrove Forests Using a 3D-SPH Method. Ocean Engineering, 306: 118001. https://doi.org/10.1016/j.oceaneng.2024.118001

[24] Peng, X. Y., Yu, P. C., Chen, G. Q., et al., 2020. Development of a Coupled DDA–SPH Method and Its Application to Dynamic Simulation of Landslides Involving Solid–Fluid Interaction. Rock Mechanics and Rock Engineering, 53(1): 113-131. https://doi.org/10.1007/s00603-019-01900-x

[25] Ren, B., Wen, H. J., Dong, P., et al., 2016. Improved SPH Simulation of Wave Motions and Turbulent Flows through Porous Media. Coastal Engineering, 107: 14-27. https://doi.org/10.1016/j.coastaleng.2015.10.004

[26] Shi, Z. M., Zhang, G. D., Peng, M., et al., 2023. Experimental Investigation on the Breaching Mechanisms of Landslide Dams with Heterogeneous Structures. Advanced Engineering Sciences, 55(1): 129-140 (in Chinese with English abstract).

[27] Viroulet, S., Sauret, A., Kimmoun, O., 2014. Tsunami Generated by a Granular Collapse down a Rough Inclined Plane. EPL (Europhysics Letters), 105(3): 34004. https://doi.org/10.1209/0295-5075/105/34004

[28] Viroulet, S., Sauret, A., Kimmoun, O., et al., 2013. Granular Collapse into Water: Toward Tsunami Landslides. Journal of Visualization, 16(3): 189-191. https://doi.org/10.1007/s12650-013-0171-4

[29] Wendland, H., 1995. Piecewise Polynomial, Positive Definite and Compactly Supported Radial Functions of Minimal Degree. Advances in Computational Mathematics, 4(1): 389-396. https://doi.org/10.1007/BF02123482

[30] Wu, H., Nian, T. K., Shan, Z. G., 2023. Research Progress on the Formation Mechanism and Risk Assessment Method of River Blocking Induced by Landslide. Chinese Journal of Rock Mechanics and Engineering, 42(S1): 3192-3205 (in Chinese with English abstract).

[31] Xu, Q., Zheng, G., Li, W. L., et al., 2018. Study on Successive Landslide Damming Events of Jinsha River in Baige Village on Octorber 11 and November 3, 2018. Journal of Engineering Geology, 26(6): 1534-1551 (in Chinese with English abstract).

[32] Zhang, C., Hu, X. Y., Adams, N. A., 2017. A Weakly Compressible SPH Method Based on a Low-Dissipation Riemann Solver. Journal of Computational Physics, 335: 605-620. https://doi.org/10.1016/j.jcp.2017.01.027

[33] Zhang, C., Rezavand, M., Zhu, Y. J., et al., 2021. SPHinXsys: An Open-Source Multi-Physics and Multi-Resolution Library Based on Smoothed Particle Hydrodynamics. Computer Physics Communications, 267: 108066. https://doi.org/10.1016/j.cpc.2021.108066

[34] Zhang, G. B., Tang, D. L., Wen, H. J., et al.,2024a. An Improved Two Phases-Two Points SPH Model for Submerged Landslide. Computers and Geotechnics, 176: 106802. https://doi.org/10.1016/j.compgeo.2024.106802

[35] Zhang, S. H., Zhang, C., Hu, X. Y.,et al.,2024b. A Riemann-Based SPH Method for Modelling Large Deformation of Granular Materials. Computers and Geotechnics, 167: 106052. https://doi.org/10.1016/j.compgeo.2023.106052

[36] Zhang, L. M., Xiao, T., He, J., et al., 2019. Erosion-Based Analysis of Breaching of Baige Landslide Dams on the Jinsha River, China, in 2018. Landslides, 16(10): 1965-1979. https://doi.org/10.1007/s10346-019-01247-y

[37] Zhou, G. G. D., Roque, P. J. C., Xie, Y. X., et al., 2020. Numerical Study on the Evolution Process of a Geohazards Chain Resulting from the Yigong Landslide. Landslides, 17(11): 2563-2576. https://doi.org/10.1007/s10346-020-01448-w

[38] Zhou, L., Fan, X. M., Xu, Q., et al., 2019. Numerical Simulation and Hazard Prediction on Movement Process Characteristics of Baige Landslide in Jinsha River. Journal of Engineering Geology, 27(6): 1395-1404 (in Chinese with English abstract).

[39] Zhu, C. W., Peng, C., Wu, W., et al., 2022. A Multi-Layer SPH Method for Generic Water-Soil Dynamic Coupling Problems. Part I: Revisit, Theory, and Validation. Computer Methods in Applied Mechanics and Engineering, 396: 115106. https://doi.org/10.1016/j.cma.2022.115106

基金资助

国家自然科学基金-联合基金重点项目(U23A2044)

广西重点研发计划项目(No.桂科AB25069121)

国家自然科学基金-青年基金(42207238)