基于带符号距离函数的摩擦接触机制，利用双相双点物质点方法模拟滑坡引发的海啸

《Computers and Geotechnics》：Frictional contact based on signed distance functions for simulating landslide-induced tsunamis by two-phase double-point material point method

　　摘要
山体滑坡引发的海啸会产生局部性强但破坏性极大的冲击波，对周边基础设施构成威胁。为此，采用了两相双点材料点方法（MPM）来分析滑坡体侵入水体过程中的固液相互作用。然而，以往通常简化了波浪流动与斜坡表面之间的摩擦接触，导致对涌浪的产生和传播预测不准确。本研究引入了符号距离函数（SDF）来定义基岩、滑坡体和水体之间的界面。通过为流体-斜坡和固体-斜坡界面分配不同的摩擦系数，该方法能够同时模拟滑坡失稳机制和波浪爬升过程。通过对两个实验案例的验证，该方法被应用于1958年利图亚湾巨型海啸的模拟。模拟结果再现了由于大量岩石以超过100米/秒的速度撞击利图亚湾而产生的极端波浪爬升高度（模拟值526米，实际观测值524米）。模拟不仅展现了从初始落石到颗粒状滑坡的转变过程，还捕捉到了滑坡侵入水体和波浪回流过程中的双向固液相互作用。此外，对流体-斜坡摩擦系数的敏感性分析揭示了斜坡表面粗糙度在评估山体滑坡海啸灾害中的重要性。总之，本研究证明，结合SDF的双相双点MPM框架为评估峡湾、湖泊和水库中的山体滑坡海啸灾害提供了可靠有效的工具。

1. 引言
山体滑坡引发的海啸是由滑坡体快速进入水体所形成的（Heller和Ruffini，2023）。与通常源自线源的地震海啸不同，山体滑坡引发的海啸会产生局部性强且极端的冲击波和波浪爬升高度，对周边社区和基础设施构成严重威胁。例如，一次体积为30×10^6立方米的陆地滑坡，其滑动中心位于610米高度处，引发了迄今为止记录中利图亚湾最高的524米波浪爬升高度（Fritz等人，2009；Ward和Day，2010）。山体滑坡引发的海啸包含多个阶段：滑坡体的初始运动和向下移动、侵入水体并与之相互作用，以及随后海啸波的产生和传播（González-Vida等人，2019；Sabeti和Heidarzadeh，2022；Jenkins等人，2025）。在此过程中，滑坡体的势能大部分转化为波浪的动能，同时伴有滑坡体/涌浪与斜坡表面之间的摩擦能量耗散（Yavari-Ramshe和Ataie-Ashtiani，2016；Romano等人，2025）。为了研究这些大变形现象，人们采用了多种数值方法，包括平滑粒子流体动力学（SPH）（Mao和Guan，2023；Fang等人，2024；Wu等人，2024a；Liu等人，2025；Yuan等人，2025）、计算流体动力学（CFD）（Li等人，2023；Romano等人，2023；Latifah等人，2024；Tofany和Latifah，2025）以及两相双点材料点方法（MPM）（Zhao等人，2023a；He等人，2024；Feng等人，2025）等。

两相双点MPM是一种混合拉格朗日-欧拉方法，能够捕捉固液相互作用（Zheng等人，2024；Xie等人，2025a），如图1所示。在该框架中，滑坡体和水分别被离散为固体和流体材料点，并为每个阶段建立两个独立的背景网格（Gs和Gf）。材料点的物理信息（如质量、速度）通过插值函数映射到相应的网格上。同时包含两阶段信息的节点被称为Gsf（图1中的红色虚线）。在解决网格节点上的动量平衡问题后（包括Gsf处的固液相互作用），更新后的节点信息被投影回材料点，从而更新整个系统的速度和动量场（Zhan等人，2023；Rébillout等人，2024）。与单相模拟不同，这种方法激活了多层背景网格。例如，Bandara和Soga（2015）以及Feng等人（2021）使用水材料点在渗流边界分析水分浸入对土壤斜坡剪切带发展的影响。尽管该框架能够结合降雨强度和3D地形模型进行区域滑坡分析，但在模拟整个滑坡过程时计算成本较高（Yamaguchi等人，2023；Zhan等人，2023），因为降雨和渗透过程可能持续数小时甚至数天。

图1. 两层网格MPM的计算过程及其在固体-斜坡和流体-斜坡界面中的影响区域。
在利用两相双点MPM模拟山体滑坡引发的海啸的背景下，Du等人（2021；Du等人，2022）研究了灾害链，并从系统能量转换的角度揭示了机制演变。Feng等人（2025）引入了用于描述海底颗粒流中非线性大变形流体-土壤相互作用的弹塑性相变模型，指出纯弹塑性土壤模型往往高估了滑坡的蠕变距离。为了解决强固液相互作用期间流体相的体积锁定问题，提出了投影误差补偿方法（Li等人，2026；Liu等人，2026）和高阶插值函数（He等人，2024；Alsardi等人，2026；Cheng等人，2026）。然而，与关注滑坡起始和蠕变沉积的研究不同，山体滑坡海啸的灾害评估重点在于波浪传播和爬升高度，因为波浪能量主要在复杂地形上的爬升过程中耗散。为了捕捉这种受摩擦主导的效果，在时间步n时，将固体-斜坡和流体-斜坡的摩擦接触力纳入山体滑坡海啸系统的动量平衡中，具体表达如下偏微分方程：

(1)
m_n vn+1 - vn Δt = m_ng + σ_n - fsfn + fn_con，
其中 m = diag msf 是网格节点上的质量矩阵，下标s和f分别代表滑坡体和水体；
v = v_sv f_t 是两相的速度矩阵；
σ = σ_s σ_f_t 是由于局部变形产生的内力；
fsf 是滑坡体侵入水体期间的固液相互作用力；
fn_con = fsn_con f_n_con t 是滑坡体/水体与基底边界之间的摩擦接触力。

方程（1）表明，固体和流体相的动量变化（左侧）由重力、惯性应变能变化、固液相互作用以及边界处的摩擦能量耗散（右侧）驱动。然而，在两相双点MPM中纳入摩擦边界带来了两个挑战：（a）在两个背景网格上准确定义复杂的流体-斜坡和固体-斜坡界面；（b）为这些界面分配相应的摩擦系数。第一个挑战在于定义固体-斜坡和流体-斜坡边界，以解决滑坡体和涌浪在复杂地形上施加的接触力（Lei等人，2022；Zhao等人，2023a；Xie等人，2025b）。在MPM中，边界条件通过两种主要方法施加：粒子边界方法和网格边界方法（Zhou等人，2025），如图2所示。在粒子边界方法（图2a）中，沿边界定义了专门的材料点（红色点），摩擦接触力通过基于位置的惩罚公式来控制。该方法可以适应不规则表面和移动的位移边界（Chandra等人，2021；Singer等人，2024；Liang等人，2024；Zhang等人，2025）。然而，对于涉及固体破碎和波浪腔体坍塌的山体滑坡海啸情景，这种方法效率较低，因为需要在两个阶段的新创建表面上生成新的粒子。网格边界方法（图2b）在背景网格上明确定义了边界附近的特定网格节点（红色节点）作为接触边界。这种方法避免了引入新的边界粒子，但将复杂地形整合到双层背景网格中成为实际应用的关键问题。Zhao等人（2023b）和Zheng等人（2026）将MPM与水平集方法结合，研究了滑坡体与复杂地形之间的接触力，但他们的方法仅适用于单层MPM。Feng等人（2025）在两相双点MPM中引入了简化的摩擦边界，但未能有效处理复杂地形。因此，在两相双点MPM框架内准确确定复杂摩擦边界仍然是有效分析山体滑坡海啸的关键问题。

图2. 粒子边界和网格边界：(a) 在接触边界中随机分布红色点；(b) 引入符号距离场（SDF）以确定每条线上的最小SDF值的网格节点。

第二个关键挑战是为山体滑坡海啸中的流体-斜坡界面分配合适的摩擦系数，特别是在狭窄的峡湾或峡谷水道中。在这种情况下，当滑坡体进入水体时，其动能传递给水体，驱动水流向对侧斜坡快速移动。这种动能随后转化为波浪的爬升能量，并通过水流在斜坡表面的流动耗散。因此，地表流动阻力在决定涌浪爬升高度中起关键作用。Nicosia等人（2024）指出，斜坡表面粗糙度和植被显著影响摩擦阻力。然而，与研究较为充分的固体-斜坡摩擦系数（0.1至0.25）不同，流体-斜坡摩擦系数通常设置为零或与固体-斜坡界面相同（Du等人，2021；Feng等人，2021；Feng等人，2025）。如何考虑斜坡表面粗糙度和植被对山体滑坡海啸中涌浪的影响仍是一个未解决的问题。

为了解决上述挑战，本研究引入了符号距离场（SDF）来定义固体-斜坡和流体-斜坡界面的摩擦接触，为每个界面分配不同的摩擦系数。SDF(i, j, k)表示点(xi, yj, zk)到形状或物体表面的最短距离（Tang等人，2025），这种技术之前用于不规则粒子之间的接触检测（Lai等人，2022；Lai等人，2023）。在最近的离散元方法（DEM）研究中（Lai等人，2024），两个接触点之间的位移是从SDF中得出的，接触力通过将该位移乘以接触弹簧刚度获得（van der Haven等人，2023；Feldfogel等人，2023；Cai等人，2025）。然而，此类SDF-DEM方法依赖于拉格朗日类型的接触检测（节点到节点、节点到面或面到面）。在本研究中，SDF用于表示固体-斜坡和流体-斜坡界面，分别计算两个相α中每个网格节点的SDF(Xgα)。如果SDF(Xgα) < 0，则该节点位于基岩内部（图2(b)中的白色点）；如果SDF(Xgα) = 0，则该节点恰好位于边界上；如果SDF(Xgα) > 0，则具有最小正SDF的节点被定义为特定相的接触边界（图2(b)中的红色点）。因此，在后两种情况下，网格节点被定义为接触界面，并根据斜坡表面粗糙度为固体-斜坡和流体-斜坡界面分配不同的摩擦系数。本研究的研究创新在于将SDF有效集成到两相双点MPM框架中，同时允许为固体-斜坡和流体-斜坡界面分配不同的摩擦系数。值得注意的是，每次迭代时更新全局SDF在计算上代价较高。相比之下，传统的基于拉格朗日网格的接触检测方法避免了这种开销，因为接触是通过网格连通性和邻元素查找自然嵌入拉格朗日框架中解决的（Neto，2023；Tang等人，2025）。因此，基于SDF的方法适用于涉及静态、刚性和相对平滑边界的情景，其中SDF可以预先计算并在整个模拟过程中重用。这种特性使得基于SDF的方法在研究山体滑坡海啸时特别有效，因为滑坡体与静止的水库边界相互作用。

本文旨在将基于符号距离场（SDF）的摩擦接触边界整合到两相双点材料点方法中，用于分析山体滑坡引发的海啸。文章的结构如下：第2节介绍了结合SDF的两相双点MPM的控制公式，用于确定固体-斜坡和流体-斜坡摩擦界面；第3节提供了两个详细的验证案例：岩石滑坡和颗粒滑坡引发的海啸。第4节进一步将这种方法应用于利图亚湾（Lituya Bay）滑坡引发的海啸的三维分析中，其中分析了巨浪产生和传播过程中的双向流体-固体相互作用，以及回流水作用下滑坡堆积物的重新分布。第5节讨论了流体-坡面摩擦系数对滑坡引发的海啸的影响。最后，第6节总结了研究的主要结论。

2. 结合SDF的MPM控制方程
为了考虑滑坡引发的海啸过程中的固体-流体相互作用，将固体、液体和气体相视为受混合理论控制的多组分连续介质。因此，单位体积定义为这三个组成相的总和：
(2) 1 - n??solid + n??liquid + n??air = 1
其中n和Sr分别代表多孔介质的孔隙度和饱和度。每个相都有其自身的速度场，分别为vα（α代表固体、液体和气体相）。这些相在空间上作为重叠的连续介质共存，形成混合物。因此，混合物的总密度和总动量由各自的相贡献之和给出（Atkin和Craine, 1976; Ding等人, 2019）：ρ = ∑αρα 和 ρv = ∑αραvα。每个相的密度定义为混合物体积中的相应质量 ρα = mα/V。在一个单位体积内，固体和液体的密度可以表示为 ρs = 1 - nρg 和 ρf = n??ρw，其中ρg和ρw分别是固体颗粒和水的密度，n和Sr表示孔隙度和饱和度（Borja, 2006）。为了降低计算成本，通过忽略气体压力和密度来简化气体相（Ceccato等人, 2021; Feng等人, 2021）。因此，在当前的公式中省略了气体相的平衡方程，未饱和条件通过空隙中的液体体积分数来考虑（Tampubolon等人, 2017）。因此，控制方程包括固体和液体相的动量平衡，以及混合物的质量平衡（Ceccato等人, 2021; Zhan等人, 2023）：
(3) ραDvαDt - ?·σα - ραg - fα = 0
(4) ?ρα?t + ?·ραvα = 0
其中σ和g分别是柯西应力张量和重力；D/Dt = ?/?t + v·?是欧拉物质导数；fα是固体和液体相之间的相互作用力，其中∑αfα = 0。通过对方程(3)和方程(4)中的固体和液体相求和，可以得到混合物的质量和动量守恒：
(5) ρDvαDt - ?·σ - ρg = 0
(6) ?ρ?t + ?·ρv = 0
其中σ = ∑ασα代表混合物的柯西应力。
在方程(3)中，考虑了固体和液体相之间的两种相互作用力：压力梯度力ff→sp和阻力ff→sd：
(7) fα = ff→sp + ff→sd
当固体相浸没在流体中时，作用压力梯度力ff→sp：
(8) ff→sp = -1 - n?p
其中p是由混合理论中的饱和度Sr加权的流体压力pf和空气压力pa组成的孔隙压力（Borja, 2006）：
(9) p = Srpf + (1 - Sr)pa
阻力ff→sd被视为当两个相之间存在速度差时固体和液体相之间的相互作用力（Bandara和Soga, 2015）：
(10) ff→sd = ceVw - vs
其中ceE = n?ρfg/k^，k^是颗粒材料的渗透率。

2.1. 两相双点物质点方法
为了模拟滑坡产生的冲击波期间固体和液体相之间的相互作用，采用两相双点物质点方法来求解控制方程(3)和方程(4)。在这个框架中，地质材料（例如岩石和土壤）由固体物质点表示，而水由液体物质点表示。在每个时间步长内，两组物质点的质量和动量分别使用插值函数映射到固体和液体背景网格上。通过考虑内部应力、边界条件和固体-流体相互作用力来更新网格节点的动量和速度。然后将更新后的节点信息映射回物质点，之后推进固体和液体点的速度和位置。

2.2. 物质点到网格的插值
在时间步长n期间，固体和液体相的粒子质量mp,αn通过权重wip,αn = Nipnxp,αn - xi转移到它们各自的网格mi,αn，其中Nipn是二次B样条插值核，xp,αn和xi分别是物质点和网格节点的位置：
(11) mi,αn = ∑pmp,αnwip,αn
(12) Nipnx = 3/4 - x2 for 0 ≤ x < 1/2
3/2 for 1/2 ≤ x < 3/2
2/3 for 2 ≤ x
其中x是无量纲距离，x = xp,αn - xi/Δx，Δx是网格间距。
采用Affine-augmented Separable FLIP（ASFLIC）方法（Fei等人, 2021）来在粒子到网格的动量传输过程中保持物质点的能量动态：
(13) mi,αnvi,αn = ∑pwip,αnmp,αnvp,αn + Cp,αnxi - xp,αn
(14) Cp,αn = ∑iwip,αn - 1vi,αnxi - xp,αn - 1
TDp,αn - 1
其中Cpα,n的初始值为Cpα,0 = 0，Dp,αn = 1/4Δx2。因此，网格节点的速度为：
(15) vi,αn = 1/mi,αn∑pwip,αnmp,αnvp,αn + Cp,αnxi - xp,αn
网格节点由于Δt的变化而产生的变形vi,αn会导致由于势能ψα的变化而产生的内部弹力σi,αn：
(16) σi,sn = -?ψs?xi,s = -∑pVp,s0?ψs?FsFp,sE,nT?wip,sn
(17) σi,fn = -?ψf?xi,f = -∑pVp,f0?ψf?JfJp,fn?wip,fn
其中Fp,αE,n是变形梯度的弹性分量，Jp,αn是Fp,αE,n的行列式。
水粒子被视为几乎不可压缩的（Becker和Teschner, 2007），势能ψf和孔隙压力确定为：
(18) ψfJf = -KfJf1 - γ
(19) σf = -?ψf?JfJf = Kf1Jfγ?1
其中Kf是水的体积模量，γ是惩罚水相体积变形的项（本研究中的γ = 3）。
对于固体相，岩石块和颗粒滑体质量被建模为使用Klár等人（2016）提出的弹性势能密度和Drucker-Prager屈服条件的弹塑性地质材料。弹性势能密度在方程(20)中用对数应变表示，Lamé参数由杨氏模量和泊松比确定。该能量密度在方程(21)中的导数为计算方程(16)中所需的弹性力σi,sn的项：
(20) ??ψs = μtr?2 + λ2tr?2
(21) ?ψs?Fs = U2μΣ?1lnΣ + λtrlnΣΣ?1VT
其中??是对数应变，??=lnΣ，U、Σ和V是变形梯度Fp,sE,n的分量。
在本研究中，使用Drucker-Prager屈服条件来描述颗粒粒子之间的库仑摩擦相互作用，其中剪切应力应小于黏聚强度和正应力乘以摩擦系数（Klár等人, 2016）：
(22) CFtrσs + ‖σs - trσsdmI‖ ≤Cc
其中CF和Cc分别是地质材料的内摩擦系数和黏聚强度，CF = ξ/2 + 2ξ2，ξ = tan?f，?f是颗粒粒子的内摩擦角；dm是模型的维度，I是单位矩阵。

2.3. 基于有符号距离场（SDF）的摩擦接触
在粒子到网格和网格到粒子的传输过程中，物理信息通过插值核映射，其半径通常超过物质点的特征尺寸。例如，方程(12)中的B样条插值核的半径为1.5Δh，Δh是网格间距。这样的核可能在滑体/水体与下层基岩之间的接触界面附近产生不利的共享影响节点（Mast等人, 2012; Nakamura等人, 2021）。这些共享节点显著降低了受影响区域内的节点速度，导致滑体和涌浪的运动锁定。一个实际的方法是根据复杂地形定义激活的边界网格节点，并在滑体/水体与基岩之间存在相对运动趋势时在这些节点上应用摩擦力。
在这项工作中，使用SDF来确定滑体/水体与基岩之间的边界网格节点，称为固体坡面和液体坡面界面，其中SDF值是为每个相α的网格节点Xgα计算的。如果SDF(Xgα) < 0，则认为节点在基岩内部；如果SDF(Xgα) = 0，则认为节点在其表面；如果SDF(Xgα) > 0，则认为节点是与SDF值最小的正SDF值的接触边界节点i。确定接触网格边界i后，通过计算网格节点的轨迹动量p～in+1/2和速度v～in+1/2来获得两个物体之间的相对运动趋势：
(23) p～in+1/2 = pin - 1/2 + Δtfin
(24) v～in+1/2 = vin - 1/2 + Δtfin
如果固体和液体相的轨迹速度v～i,αn+1/2和边界v～i,bn+1/2满足：
(25) v～i,αn+1/2 - v～i,bn+1/2niαn ≥ 0
其中niαn是接触网格节点的法向量，通过Huang等人（2011）提出的质量梯度获得：
(26) n^iαn+1 = ∑pwip,αnmp,α
n^ibn+1 = ∑pwip,bnmp,b
∑pwip,bnmp,b
为了避免基岩和滑体/水相之间的法向量不共线，然后修改法向量：
(27) niαn+1 = -niβn+1 = n^iαn+1 - n^ibn+1
两个物体之间的接触产生表面fin处的接触力con，网格节点的速度和动量可以更新为：
(28) pin+1/2 = p～in+1/2 + Δtfin,con
(29) vin+1/2 = vin-1/2 + Δtfin,con
对于无滑动接触，方程(23)和方程(24)可以给出无滑动接触力：
(30) fiαn,con = -fibn,con = mi,αnmi,bnmi,αn+mi,bn
因此，可以得到接触力的法向和切向分量：
(32) fiαn,nor = fiαn,con?niαn
(33) fiαn,tan = fiαn,con - fiαn,con?niαn
当切向摩擦超过接触面的最大静摩擦力时，‖fiαn,tan‖ > μα‖fiαn,nor‖，会发生滑动摩擦，可以通过库仑摩擦模型描述：
(34) fiαn,con = fiαn,nor + μα‖fiαn,nor‖fiαn,tan‖
其中μα是固体坡面和液体坡面界面的摩擦系数。

2.4. 更新网格动量和物质点
在粒子将动量传输到网格节点后，网格节点xi,α的位置将发生Δtvi,α的变形，更新后的位置为x^i,α = xi,α + Δtvi,α。因此，每个相中网格节点的动量守恒可以离散化为：
(35) mi,snvi,sn+1 - vi,sn = Δtσi,sn + fi,sn
(36) mj,fnvj,fn+1 - vj,fn = Δtσj,fn + fjfn,con + mj,fng - fsfn
忽略压力梯度力后，可以得到固体-流体相互作用项的离散形式：
(37) fsfn = -cEmi,snmj,fnvi,sn - vj,fn
假设质量矩阵m = ms00，速度v = vsv，内力σ = σsσf，接触力fn,con = fsn,con，以及固体-流体相互作用fsf，可以得到固体-流体耦合系统：
(38) mn+Δtfsfnvn+1 = mnvn+Δtmng + σn+fn,con
通过求解这个系统，可以获得每个相中网格节点的速度，从而更新变形梯度、位置和粒子的速度。

2.5. 更新物质点
然后将网格节点的速度场根据ASFLIP方案（Fei等人, 2021）传输回两个相中的物质点，速度空间导数矩阵Cp,αn为：
(39) vp,αn+1 = ∑iwip,αnvi,αn+1 + apvp,αn - ∑iwip,αnvi,αn
(40) xp,αn+1 = xp,αn + Δt∑iwip,αnvi,αn+1 + apbpvp,αn - ∑iwip,αnvi,αn
(41) Cp,αn+1 = ∑iwip,αnvi,αn+1xi - xp,αnTDp,αn-1
其中ap和bp是PIC-FLIP混合比率（Bridson 2015）和位置的调整比率（Fei等人, 2021），在本研究中ap = 0.999，bp = 0.3。值ap = 0.999对应于FLIP主导的更新，有助于在高动态滑坡-海啸模拟中减少过度的数值耗散，同时保留少量的PIC贡献以增强数值稳健性。系数bp = 0.3在ASFLIP方案中用作位置调整参数，遵循Fei等人（2021）。流体粒子的行列式通过以下方程更新：
(42) Jp,fn+1 = I + Δttr?vp,fn+1
固体粒子的变形梯度通过（Stomakhin等人, 2013）更新：
(43) Fp,sn+1 = I + Δt?vp,sn+1
对于沙粒上的水饱和度，采用Tampubolon等人（2017）的近似：
(44) Srn+1 = ∑iwip,sn?i,sn+1
其中?i,sn+1是一个项，当这个节点对水和沙粒都有非零质量时?i,sn+1 = 1，否则?i,sn+1 = 0。
总体而言，在方程(34)中定义的固体坡面和液体坡面界面的摩擦接触力被纳入两相双点物质点方法框架中。这种配置对于模拟滑坡引发的海啸至关重要，因为滑体与坡面的相互作用涉及动态摩擦、碰撞和固体材料的破碎，而水波在相反坡面上的上升特性是由过渡性和湍流性地表流（Nicosia等人, 2024）表征的。固体坡面和液体坡面的摩擦系数可能有几个数量级的差异，例如，液体坡面的摩擦系数为10?3（Nicosia等人, 2024），而固体坡面的摩擦系数为10?1（Zhao等人, 2023a）。与依赖单一摩擦系数的现有MPM模型相比，采用两个不同的系数能够更全面地研究与滑坡引发的海啸相关的全部危险链。

3. 用模型测试验证当前方法
3.1. S?levik等人（2009）的滑坡引发的海啸
在1934年由滑坡引发的Tafjord海啸之后，调查人员在Akneset fjord支流中识别出多个不稳定的岩石质量（S?levik等人, 2009）。为了研究这种撞击产生的冲击波，S?levik等人（2009）在尺寸为25米（长度）× 0.51米（宽度）× 1米（深度）的波浪池中进行了滑坡入水的物理模型测试。该水池具有透明的玻璃底部和侧面，可以使用高速相机和粒子图像测速（PIV）捕捉波浪运动。实验设置在1:500的模型比例下，如图5所示。水深为0.6米，滑面以35°的角度与平衡水面相交，固体坡面的摩擦系数为0.14。由于流体-坡面界面是玻璃的，其摩擦系数设为0.01，以考虑玻璃平面上水的粘度（Nicosia等人, 2024）。滑块尺寸为1.0米（长度）× 0.45米（高度）× 0.16米（厚度），前端角度为45°。沿斜面的冲击速度达到了3.38米每秒，由此产生的波高（η）通过三个声波测量仪记录下来，这些测量仪分别位于距离点O 3.25米、4.24米和6.55米的位置。为了分别重建MPM中的固体和流体域，采用了一种数字图像处理技术（Wu等人，2024b）。首先将相位域转换为带有RGB值的像素阵列，RGB值定义颜色，灰度级别表示亮度，从而区分固体相和流体相以及斜面边界。然后通过粒子采样将固体和流体相的像素场转换为物质点。根据相应的RGB值为这些物质点分配机械属性，如密度、内聚力和摩擦角。斜面由图3(a)中的粗斜线表示，被设定为边界。使用了一个开源代码仓库（Yu和Shang，2022）来确定斜面边界附近的网格节点，从而为两相双点MPM提供了固体-斜面和流体-斜面界面，如图3(b)所示。该模型中的网格间距是由基于图像的离散化确定的，总长度15米被分为256个间隔，得到0.0586米的精确间距和20,865个物质点。岩石块参数列在表1中，而水的属性包括密度1000千克/立方米和体积模量2×10^8帕斯卡。图3(c)展示了水箱中的孔隙压力，这与通过ρfgz计算出的理论静水压力非常吻合。

下载：下载高分辨率图像（278KB）
下载：下载全尺寸图像

图3. 岩石滑坡-海啸案例中滑移面和波波仪的示意图，其中d表示水深，η表示波高，H表示绝对波高。

表1. 用于岩石滑坡-海啸案例的MPM建模的滑体参数。
高度（米）
长度（米）
密度（千克/立方米）
速度（米/秒）
内聚力（千帕）
弹性模量（兆帕）
摩擦角（°）
泊松比
0.16
1
12
4
8
3.38
2
30.77
2×10^4
4
0.3

图4通过选定时间步长的水平速度等值线展示了滑动岩石块引发的海啸生成过程。当与水面碰撞时，滑体引发了流速的迅速增加，并形成了沿正x方向传播的第一波，其中包含一个高度超过0.5米的水舌（图4b）。在大约t=0.75秒时，水舌重新塌陷回水箱中，在正负两个方向上产生了次级流动。向前流动（正向）与第一波合并并放大它，而向后流动（负向）从边界反射并返回岩石块，随后转向正向（图4c）。这种反射促成了第二波的形成，在t=1.75秒后变得可见（图4d）。第一波和第二波的峰值之间的距离大约为3.8米。在波浪传播过程中，波高主要由重力（作为主要恢复力）和惯性控制。当波浪能量将水体推向波峰时，重力将其拉回下方；然后惯性使其超过平衡位置，形成一个低谷。这两种波浪的相互作用和传播在水面上引发了进一步的振荡。随着时间的推移，这些波动逐渐减弱，直到水面恢复到静止水位。

下载：下载高分辨率图像（357KB）
下载：下载全尺寸图像
下载：下载高分辨率图像（120KB）
下载：下载全尺寸图像

图4. 本发明方法模拟的在不同时间步长的岩石滑坡引发的海啸的水平速度等值线。

图5比较了三个观测点监测到的波高，这些数据包括实验数据（S?levik等人，2009）、SPH结果（Qi等人，2022）以及本发明方法，使用无量纲波高η/d和无量纲时间tg/d进行比较。在图5(a)中，岩石块进入水中后，第一个观测点的波高在第一波过后上升至η/d=0.21，然后由于重力作用降至-0.12。第二波在同一位置引发了第二次上升，达到η/d=0.08。与第一点相比，第二点和第三点的波长较长但振幅较小（图5b和5c）。需要注意的是，尽管采用了ASFLIP方案来改善粒子-网格传递过程中的动态特征守恒，但数值耗散无法完全消除（Fei等人，2021）。随着模拟的进行，这种累积的数值耗散导致残余的水面波动比实际实验中观察到的衰减得更快。因此，模拟的波高在tg/d>20后趋于更加平滑，并接近静止水位。同时，记录结果和模拟结果之间的峰值波高的均方根误差（RMSE）为RMSEk=∑i=1nηi-η^i,k^2/d^2n，其中k表示数值方法（SPH或MPM）。对于MPM模拟，14个峰值产生的RMSEMPM=0.03597；对于SPH，5个峰值产生的RMSESPH=0.06434。总体而言，MPM预测的波高与实验测量结果吻合得很好，从而验证了本发明的数值方法。

下载：下载高分辨率图像（522KB）
下载：下载全尺寸图像

图5. 实验观测和数值结果中三个观测点的无量纲波高比较。

图6比较了岩石块切入水箱过程中的水舌形态，左侧列显示了实验中的PIV图像（S?levik等人，2009），右侧列显示了MPM模拟结果。当滑体进入水面时，其动能传递给水体，使其向上形成明显的水舌。模拟的水舌与实验观察结果非常吻合，特别是在滑体进入水面时形成的接触角方面。

下载：下载高分辨率图像（679KB）
下载：下载全尺寸图像

图6. 左侧列显示了实验结果（S?levik等人，2009），右侧列显示了MPM模拟结果，用于描述岩石块切入水体过程中的水舌。

3.2. Fritz等人（2001）研究的颗粒滑坡引发的海啸
1958年，阿拉斯加海岸的Fairweather断层发生了一次M8.3级地震，引发了Lituya湾的岩石滑坡。由此产生的岩石雪崩进入水中并引发海啸，冲击了对岸，记录到的浪高达到了524米，对沿海森林造成了严重破坏（图7a中的森林修剪线）。岩石滑坡起源于海拔915米、平均坡度为45°的活跃塌陷区。滑坡物质的固有密度为2700千克/立方米，长度约为970米，最大厚度约为92米，重心高度H约为610米。水深d为122米，静止水面与滑坡斜面的角度为45°。为了研究这一事件，Fritz等人（2001）进行了一个实验室实验，模拟了峡湾的NE-SW方向的剖面。实验在一个长11米、宽0.5米、高1米的矩形水箱中进行。模型比例（ms）为675，总模型长度为4.2米，包括1.98米的水面长度。模型中的水深为0.18米，静止水面与滑坡斜面的角度保持在45°。滑坡质量使用平均粒径为4毫米的颗粒材料建模，孔隙率为0.4，密度为2640千克/立方米，与原位岩石密度非常接近。模型滑坡的长度为0.61米，平均厚度为0.18米。监测点位于滑体入水点1.31米处。实验装置如图7(b)所示，其中ηR表示岬角坡上的波浪上升高度。网格间距为0.0164米，该数值模型中有43,036个物质点。颗粒材料的模拟参数如下：杨氏模量200兆帕，泊松比0.3，内聚力0.5千帕，内摩擦角43°，渗透率为1×10^-3米/秒，固体-斜面界面的摩擦系数为0.1（Fritz等人，2009），流体-斜面界面的摩擦系数为0.01。

下载：下载高分辨率图像（360KB）
下载：下载全尺寸图像

图7. 基于NE-SW方向剖面和Lituya滑坡参数的模型测试示意图（Fritz等人，2001）。

图8展示了模拟的滑坡引发的海啸及其对应的水平速度等值线。释放后，滑体变形，形成一个平缓倾斜的前沿，然后在t=0.45秒时撞击水面，并在t=0.65秒时抵住沟底发生偏移，导致显著的滑体膨胀。在此阶段，最大的水平流速（1.65米/秒）集中在滑体-水体界面。快速进入水中导致在滑体肩部形成水流分离，在t=0.75秒时在颗粒质量后面形成了一个大的水舌和空气腔。相关的动量交换使最大水速在固体-流体相互作用界面升高至1.78米/秒。随着滑体沿着沟底继续移动，崩溃的水舌促进了水体和碎片的混合。领先波，即一个能量脉冲，在t=1秒和t=1.25秒时传播并上升到对岸坡。这种上升代表了水流从动能向势能的转换，降低了流速，并在t=2秒时导致水面向岬角坡上升形成一层薄水层。图9展示了颗粒滑坡在渗透过程中饱和度的相应演变。最初干燥的滑体在水冲击后其接触边界处的饱和度增加。t=0.75秒时的高速冲击引发了滑体前端的沉积物运输，表现为水舌中的饱和颗粒，这是剥离效应的表现（P?htz等人，2020）。腔体崩溃后，持续的水-碎片混合进一步将饱和度向滑体内部推进。

下载：下载高分辨率图像（280KB）
下载：下载全尺寸图像

图8. 本发明方法获得的滑坡引发的海啸案例的水面模式和速度等值线。

图9. 颗粒物质在滑体渗入水中过程中的饱和度变化。

图10进一步比较了Fritz等人（2001）的实验结果和本发明数值模拟之间的颗粒滑坡进入过程，展示了从滑体渗透到随后波浪生成和传播的运动传递。实验中的PIV图像揭示了过程中的两个不同阶段：（1）冲击和分离阶段：冲击发生后，滑体渗透速度超过波浪传播速度，导致颗粒质量与水分离，形成明显的水舌和相关的空气腔，如图10中t=0.39秒到t=0.58秒所观察到的。（2）腔体崩溃和波浪超越阶段：随着滑体继续沿沟道移动，空气腔体崩溃，传播的波浪超越了滑体肩部，如t=0.65秒到t=0.715秒所示。应当注意的是，在该研究中，模型设置过程中忽略了气体相，这可能会影响腔体形成、气泡破裂、相混合以及水进入时的冲击力的分析（Fritz等人，2009）。这种简化可以显著降低计算成本，因为滑体和质量只占计算域的一小部分。总体而言，MPM模拟得到的结果与实验观察结果吻合得很好。

下载：下载高分辨率图像（451KB）
下载：下载全尺寸图像

图10. 颗粒滑坡引发的海啸模型测试中滑体进入和舌形态的比较。

图11比较了（a）x=1.31米处的波高（原型x=885米）和（b）对岸坡上的波浪上升高度的数值和实验结果。在这个二维设置中，滑坡引发的波浪向前传播并在滑体和岬角坡之间反复反射。为了与Fritz等人（2001）的结果直接比较，模拟的波高和上升高度都按模型比例（ms=675）进行了缩放，而时间按tg/d进行了归一化。滑体渗透引起的动量传递驱动波浪在正x方向上传播，产生了第一个波高峰，其缩放值为η·ms=152米（图11a）。然后波浪上升到对岸坡，达到最大缩放上升高度ηR·ms=526米（图11b）。这个值与实际事件后记录的森林破坏线（524米）非常吻合。在上升流动之后，水从坡道上回流到河道中，产生一个向负x方向传播的波浪，并在图11(a)中形成一个第二个峰值（η·ms = 85米）。值得注意的是，在实验中的涌浪上升过程中，一些水也会绕过滑动质量的后部流动。这种环绕水流阻碍了涌浪的回流，从而减缓了上升高度的减小（Fritz等人，2009年）。与三维实验不同，这种环绕流动在二维数值模拟中不会发生，这导致了在本研究中使用的三维物理实验和二维数值模拟之间上升高度的差异（tg/d > 15）。总体而言，数值结果与实验数据非常吻合，验证了目前方法在模拟颗粒滑动引发的海啸方面的准确性。

下载：下载高分辨率图像（322KB）
下载：下载全尺寸图像

图11. 当前数值结果与Fritz等人（2001年）的实验结果相比，观察点的波高和岬角坡道上的波浪上升高度。

4. 利图亚湾巨型海啸的案例研究
在本节中，我们将第3.2节中利的图亚湾滑坡引发的海啸的二维模型扩展为三维研究，通过引入美国地质调查局的数字高程模型作为摩擦边界。滑动质量在横截面中被表示为一个近似三角棱柱体，沿斜坡延伸730至915米，最大厚度约为92米，质心高度为610米。滑动质量的总体积为33×10^6立方米，与Miller（1960年）报告的30.6×10^6立方米估计值一致。使用了美国海岸和大地测量局的水深数据，揭示了一个U形沟槽，其墙壁陡峭，底部平坦，最大深度为120米（González-Vida等人，2019年），如图12所示。

下载：下载高分辨率图像（105KB）
下载：下载全尺寸图像

图12. 基于美国地质调查局数字高程模型的利图亚湾海啸的三维模型。
该模拟包括7.4×10^6个物质点，网格间距为20米，时间步长为10^-3秒。它是在Intel i9-12900 CPU（16核@5.1 GHz，64 GB RAM）和GeForce RTX 3090 GPU（24 GB RAM）上运行的。在GPU上每个步骤的平均运行时间为26.5秒，而边界网格节点的初始SDF评估的计算开销为27.5秒。一旦在第一步中确定了边界网格节点，在接下来的计算步骤中只会激活边界内的计算域。滑动质量的属性包括密度2650 kg/m3、杨氏模量5.4×10^9 Pa、泊松比0.2、内摩擦角30°和渗透系数0.1 m/s。固体-斜面界面的摩擦系数为0.15。水的模型密度为1000 kg/m3，体积模量2×10^8 Pa，流体-斜面界面的摩擦系数为0.12，这考虑了表面粗糙度和植被的影响（Nicosia等人，2024年）。为了监测速度演变，在滑动质量上选择了三个点Ps1（915米）、Ps2（610米）和Ps3（230米）。此外，还跟踪了三个水粒子Pw1、Pw2和Pw3，它们位于Gilbert Inlet西南侧的小山上，海拔分别为524米、450米和200米（见图12）。

图13展示了Gilbert Inlet中巨浪的生成和演变过程，固体和水粒子的速度等值线范围从0 m/s到100 m/s。图14展示了三个固体粒子（Ps1、Ps2和Ps3）和三个水粒子（Pw1、Pw2和Pw3）的速度历史。从静止状态释放后，由于势能转换为动能，滑动质量加速，在水进入之前达到70 m/s（图13，t = 13秒）。当碎片进入水中时，由于水的相互作用力，浸没的部分减速，如图14(a)中t = 14.5秒时粒子Ps3的速度突然降至51 m/s所示。相比之下，水面以上的碎片（例如Ps1）继续加速，到t = 25.5秒时达到115 m/s，这是由于前方粒子的重新加速。值得注意的是，质心粒子Ps2在水进入之前达到了108.2 m/s的峰值速度，这一速度与Miller（1960年）报告的总体滑坡速度110 m/s、Fritz等人（2009年）的实验结果以及González-Vida等人（2019年）的CFD模拟结果一致。滑动质量进入水中后，水粒子的速度急剧增加，在t = 28秒时达到最大波高225米，这与Franco等人（2020年）的数值估计208米和Pastor等人（2008年）的226米一致，流速范围为60至70 m/s（图13）。如图14(b)所示，三个被跟踪的水粒子在上升之前迅速加速到大约65 m/s。然后波浪冲击对面的岬角，产生最大上升高度526米，这与现场观察到的 trimming line高度524米非常吻合（图13，t = 70秒）。在41秒到52秒之间，部分涌浪越过小山顶部，例如图14(b)中t = 40秒时Pw3的速度达到76.6 m/s，t = 46秒时Pw1的速度达到56.3 m/s，其余部分在49秒后回流到Gilbert Inlet。溢出的水粒子随后沿斜坡向下流动，势能重新转换为动能，Pw3、Pw2和Pw1的速度分别在t = 46秒、49秒和64秒时达到103.4 m/s、110.8 m/s和114.6 m/s。

下载：下载高分辨率图像（366KB）
下载：下载全尺寸图像

图13. Gilbert Inlet中由滑坡引发的巨浪的生成和传播。

下载：下载高分辨率图像（204KB）
下载：下载全尺寸图像

图14. 三个固体粒子Ps1、Ps2和Ps3位于滑动质量的海拔915米、615米和230米处，以及三个水粒子Pw1、Pw2和Pw3分别越过斜坡顶部海拔524米、450米和200米处的速度发展。

图15进一步展示了图12中所示的I-Ⅰ截面在滑坡进入水中和随后的涌浪回流到Gilbert Inlet过程中的流体-固体相互作用。滑坡质量进入后，不仅产生了225米的波高，还使整个水墙的过渡速度加快到60 m/s。到t = 36秒时，水表面高度达到对面斜坡脚下的海平面以上150米，然后海啸波涌向岬角斜坡，部分波浪在t = 46秒和53秒时越过小山顶部。值得注意的是，水前端以及用红色虚线标出的部分滑动质量也在t = 28秒时沿着冰川渠道向上流动。当波浪达到最大上升高度后，水以超过50 m/s的峰值速度冲回Gilbert Inlet，冲刷并将滑坡沉积物推向滑坡斜坡。因此，这两种方向的流体-固体相互作用不仅产生了巨型海啸和巨浪的传播，还重新分配了水中的滑坡沉积物。

下载：下载高分辨率图像（188KB）
下载：下载全尺寸图像

图15. 滑坡质量进入水中和涌入Gilbert Inlet过程中的双向流体-固体相互作用。

图16展示了整个事件中滑动质量内累积塑性应变的发展。由Fairweather断层上的Mw 8.3级地震触发，Miller（1960年）和Fritz等人（2009年）提出，整个岩石质量作为连贯单元迅速进入海洋是产生最大记录海啸的主要机制。然而，这种质量运动的精确分类仍存在争议：Pararas-Carayannis（1999年）将其归类为落石，Fritz等人（2001年、2009年）归类为滑坡，而Miller（1960年）认为两者兼有。在本研究中，图16展示了质量运动类型的时间变化，从最初的落石阶段（0-15秒）发展到完全发展的滑坡（30秒后），这一过程由岩石碎裂的逐步进行所驱动，表现为塑性应变的发展。启动后，由于岩石块的粘聚强度，整个质量作为一个连贯单元沿斜坡向下移动。在基底固体-斜面界面的摩擦接触引发了脱离区的岩石碎裂，这一点在10秒和15秒时沿接触边界积累的塑性应变中得到了证明。随着碎裂的进行，到30秒时明显出现了颗粒滑动的转变。在25秒到40秒之间，颗粒质量滑入利图亚湾并向对面斜坡扩散。随后，在50秒到80秒之间，上升波的回流将颗粒沉积物推回原始斜坡，这表明流体-固体相互作用不仅在初始滑坡穿透期间发生，在回流期间也发生。最终沉积物在80秒时的情况与US Coast and Geodetic Survey（González-Vida等人，2019年报告）提供的此次事件之前（基于1926年数据）和之后（基于1959年数据）的利图亚湾水位变化进行了比较。可接受的一致性表明，该模拟合理地再现了水进入和回流影响下的最终沉积物形态。

下载：下载高分辨率图像（356KB）
下载：下载全尺寸图像

图16. 滑坡开始和进入水中期间滑动质量中的累积塑性应变。

应当指出，尽管典型的岩石接触摩擦系数通常在大约0.3–0.8的范围内（Blanpied等人，1995年），但快速的利图亚湾滑坡涉及到剧烈的碎裂、碰撞、滚动和岩石质量的滑动。这些过程可以降低宏观的基底阻力，并由于滑坡剪切带中的颗粒间摩擦减弱而提高滑坡的移动性（Chen等人，2025年）。在本研究中，进行了一系列具有不同固体-斜面摩擦系数的试验模拟，μs = 0.15提供了与观察到的滑坡高移动性、估计的冲击速度和记录的上升高度的一致性。因此，μs = 0.15被视为一个有效的基底摩擦系数，它结合了颗粒断裂、颗粒间摩擦减弱、滚动和碰撞效应。

5. 讨论
本研究在两相双点MPM框架内发展了一种摩擦接触公式，用于分析由滑坡引发的海啸。所提出的方法使用SDF来识别基岩和滑动质量/水之间的接触边界，从而可以根据斜坡的表面粗糙度和植被情况分配固体-斜面和流体-斜面界面的摩擦系数。该方法通过两个实验基准案例进行了验证，确认了其模拟由连贯岩石块（S?levik等人，2009年）和颗粒滑动（Fritz等人，2001年）引发的海啸的适用性。这些验证证明了该方法捕捉控制滑坡引发海啸中复杂相互作用能量传递机制的能力。利图亚湾巨型海啸的案例研究进一步表明，极端的上升高度是由滑坡的特定运动学和随之产生的波浪动力学造成的。大量岩石体积（约33×10^6立方米）以超过100 m/s的高速冲击进入利图亚湾，产生了最大振幅为225米的巨浪，随后在入口对面的小山脊上产生了大约526米的最大上升高度。这一数值与观察到的森林破坏的trimline高度524米非常一致。

此外，本研究解决了早期几项数值研究中未能准确再现岩石质量破坏过程的限制。虽然通过“更密集的流体”类比获得了有价值的见解来模拟滑坡的高移动性（González-Vida等人，2019年；Franco等人，2020年），但基于CFD的这种方法本质上难以捕捉脱离区内的机械破坏以及滑坡传播过程中的 subsequent 碎裂。相比之下，所提出的方法结合了Drucker-Prager屈服准则来描述岩石质量的破坏机制，允许岩石质量作为一个连贯单元进入水中并产生足够的波浪能量。此外，通过这种两相双点MPM观察到的流体-固体相互作用不仅发生在颗粒质量进入水道期间，也发生在涌浪回流期间，这与利图亚湾中的滑坡颗粒沉积物相互作用。

尽管总体上与现场观察和实验结果一致，但必须承认本研究的某些限制。首先，该模型没有考虑到引发巨大山体滑坡的机制——即Fairweather断层上的8.3级地震，因此没有明确模拟由地震引起的滑动面的形成和启动过程。其次，流体-斜坡界面的摩擦系数主要是通过敏感性分析来确定的，以克服参数不确定性（González-Vida等人，2019年），而不是基于现场规模的实验来获取流体流动的摩擦阻力（Evers等人，2019年；Heller和Hager，2010年；Nicosia等人，2024年）。尽管如此，考虑到这次有记录以来最大的海啸的规模，通过敏感性分析确定流体-斜坡摩擦系数，并根据Miller（1960年）、Fritz等人（2009年）和González-Vida等人（2019年）的研究结果对其进行校准是合理的。图17展示了流体-斜坡摩擦系数ξf（范围从0到1）对滑动质量中质心粒子的速度、最大波高、涌浪的上升高度以及模拟的破坏线的影响。图17(a)显示，在滑动质量穿透前的0到16秒内，不同情况下的速度曲线重叠，表明ξf对滑动质量的势能到动能转换影响很小。18秒后，较高的ξf值导致由于水流摩擦阻力的增加，滑动质量的速度下降得更快。由于水流摩擦阻力导致的滑动质量动能的耗散，最大波高从246米（ξf = 0）降低到197米（ξf = 1），如图17(b)所示。因此，涌浪的上升高度从718米（ξf = 0）降低到142米（ξf = 1），这不仅是因为从滑坡质量传递到水中的动能减少，还因为在上升过程中涌浪与斜坡表面之间的能量耗散。图17(d)中观察到ξf = 1时的森林破坏线明显较低，而ξf = 0.12的结果与观测到的破坏线相符。因此，ξf对上升高度的影响（从ξf = 0到ξf = 1减少了80.2%）大于对最大波高的影响（减少了19.9%）。适当的ξf值可以避免过高估计（ξf = 0）或过低估计（ξf = 1）滑坡引发的海啸灾害。

通过上述敏感性分析，Lituya湾案例研究中流体-斜坡界面的摩擦系数确定为0.12。应当注意的是，与实验设置中的ξf = 0.01（S?levik等人，2009年；Fritz等人，2001年）相比，其中流体-固体界面被建模为光滑玻璃，Lituya湾的摩擦系数高一个数量级。这是合理的，因为在流体-斜坡摩擦系数中考虑了对面岬角的表面粗糙度以及森林对斜坡的摩擦阻力。因此，为了评估滑坡引发的海啸在狭窄峡湾或峡谷水道中的影响，在确定流体-斜坡界面的摩擦系数时，应考虑斜坡的表面粗糙度和植被。

需要指出的是，本研究没有明确考虑关键输入参数的不确定性，包括滑动几何形状（例如体积和前沿角度）、机械性能（例如内摩擦角、凝聚力）以及水深分辨率。这些不确定性会通过模拟传播，并影响预测的上升高度和最大波幅。例如，滑动体积或释放位置的小变化已被证明可以在近场区域使海啸波幅增加高达30%（Heller和Spinneken，2015年）。同样，水的流变参数的不确定性可以使波形变化25%（Romano等人，2023年）。特别是在沿海地区，水深误差可能导致由于共振或浅滩效应而严重高估或低估上升高度（L?vholt等人，2015年）。因此，未来的工作应进行敏感性分析或概率评估（例如蒙特卡洛模拟），以量化输入变化对最终结果的影响。

6. 结论

本研究结合了基于SDF的摩擦接触公式和两相双点MPM方法，提供了一种在复杂地形中定义流体-斜坡和固体-斜坡界面并分配摩擦系数的稳健而准确的方法。所提出的方法通过与岩石块和颗粒质量穿透水道的实验结果进行了验证，在观测点的波高和渗透引起的水波形态方面显示出很好的一致性。将其应用于三维Lituya湾案例，并结合数字高程模型，证明该方法可以预测由不稳定岩坡在峡湾、湖泊或水库中产生的潜在海啸，其中冲击波威胁人类生命和基础设施。初步的敏感性分析进一步表明，流体-斜坡摩擦系数显著影响固体-流体渗透过程中的波浪动力学以及随后的上升涌浪过程。未来的工作将专注于在不同地形和流动条件下系统地量化这一系数，并将该方法扩展到实时灾害评估中。

**作者贡献声明**

朱松：撰写原始草稿、获取资金、概念化。
徐豪：软件、方法学、数据管理。
邓旭新：验证、项目管理、调查。
于冰冰：可视化、验证。
宋先照：形式分析、数据管理。

朱松 | 徐豪 | 邓旭新 | 于冰冰 | 宋先照
南京理工大学安全科学与工程学院（应急管理学院），中国南京210094