2024, Vol. 56
近岸或海底滑坡失稳快速滑移冲击水体可在水库或海洋环境诱发重大涌浪次生灾害[1-2]。滑坡产生涌浪的规模仅次于地震海啸，能够在沿河岸或海岸引发几米至数十米的浪高，具有传播速度快、影响范围广、灾害后果严重的特点[3-4]。目前，国内外已经报道了大量滑坡涌浪引发的灾害性事件，例如1958年美国Alaska州Lituya Bay冰川滑坡涌浪灾害[5-6]，2007年加拿大Chehalis Lake岩崩滑坡涌浪事件[7]，以及中国三峡水库2008年和2015年相继发生的龚家方和红岩子滑坡涌浪灾害[8-9]等。此外，海洋地质调查显示全球在历史上曾发生过大量巨型海底滑坡诱发海啸事件，例如1888年巴布几内亚Ritter Island火山滑坡海啸灾难[10]，2002年意大利Stromboli海底滑坡诱发涌浪灾害事件[11-12]，以及2018年印度尼西亚Palu Bay海底滑坡产生海啸灾害[13]等。随着全球地质环境日益复杂化，滑坡涌浪灾害事件频发且往往伴随灾难性后果，因此，研究滑坡涌浪灾害动力学机制，发展精准有效的预测方法具有重要的意义。
滑坡涌浪的产生过程(wave generation)是滑坡涌浪灾害链研究的关键和难点，涉及到滑坡、水体及空气之间的动态、多相、耦合作用，物理过程比较复杂。在目前的所有研究手段中，数值模拟可以综合考虑三者的相互作用，提供动态的涌浪产生及演化结果。目前主流的滑坡涌浪数值模拟方法主要涉及到3个主要问题：1)控制方程的选择；2)控制方程的数值解析方法；3)运动方程中滑体与水体的力学机理及相应的力学模型。控制滑体和水体运动的方程主要有Navier-Stokes方程，Boussinesq波动方程、浅水波方程等等。例如，文献[14]使用Navier-Stokes方程与本构模型相结合来描述滑坡、水和空气3种类型的流体运动；文献[15]基于Navier-Stokes方程，采用k-ε湍流模型构建流固耦合模型来模拟波浪的产生；文献[16-17]基于能够考虑波浪色散变形及非线性变形等特性的高阶Boussinesq型方程，研究近岸海域波浪产生及传播特性；文献[18]通过在浅水方程中添加滑块运动项建立了一种有限差分滑坡涌浪模型，进行水下滑坡涌浪产生过程研究。常用的数值求解方法有流体体积法(VOF)、有限元法(FEM)、有限差分(FDM)、有限体积法(FVM)、边界元法(BEM)、光滑颗粒流(SPH)、离散元(DEM)等。运动过程中的力学模型往往嵌套在运动方程中，不同的控制方程和解析手段对滑体与水体相互作用的动力学机理也会有不同的考虑。例如，文献[19]采用FLOW-3D中的通用刚体运动模型(GMO)来模拟滑坡体运动，滑坡体所受到的力包括重力、流体的阻力与壁面摩擦阻力。滑体运动至水下阶段，忽略了水体对滑坡体的阻力作用，采用给定滑坡体运动轨迹的方式进行计算；文献[20]在满足连续方程和动量方程的基础上，考虑滑体的运动受到摩擦阻力和紊流阻力，引入了Realizable k-ε紊流模型，通过VOF追踪各相之间的界面，以确定相的体积分数，最后通过每个单元中的密度和黏性系数判定流体的属性，从而更准确地模拟涌浪的全过程；文献[21]提出一种SPH-DEM耦合算法，考虑流体压力、重力、颗粒间碰撞力、表面张力用于处理宏观尺度下离散体与流体之间的相互作用问题；文献[22]采用耦合有限-离散元法(FDEM)与欧拉-拉格朗日法(CEL)的数值方法模拟了可破裂岩体的滑坡涌浪过程，首先采用FDEM考虑岩石的弹塑性、拉压荷载与剪切荷载模拟岩体的破裂过程，其次基于欧拉方法在流体运动方程中考虑了流体压力和流体黏滞性。以上方法多聚焦于如何构建有效的数值方法求解控制方程、或者创新性提出多相耦合方法，但是对滑体-水体相互作用力学机制及模型的精细化研究不足，可能导致对涌浪高度或能量产生更大或更小的估计。
文献[8,23]提出了一种全新的流体及类流体数值模拟方法Tsunami Square(TS)，该方法对运动的控制直接基于质量和动量守恒关系，免解控制方程，因此在TS框架下的数值模拟可以更多的关注模拟对象的动力学机制与模型。文献[23]基于TS提出了非动量传递(NMT)的体积抬升效应(LU)和界面摩擦效应(DA)共同作用产生涌浪的力学模型，二者分别对应于文献[24]提出的解析理论中的体积涌浪和冲击涌浪；文献[25-26]在LU与DA的基础上，补充提出了冲击挤压(push ahead，PA)力学效应及其解析模型，完善了模型的应用场景，提高了滑坡涌浪预测精度；文献[27]基于LU、PA和DA 3种力学效应及模型，建立3种力学耦合作用下的滑坡涌浪产生过程动力学模型及参数，探讨3种力学效应与16种首浪高度影响因素之间的力学关系，进一步验证所提出的动力学模型的可靠性。上述研究表明，对于TS框架下3种滑体-水体相互作用力学效应研究已取得大量成果。然而，目前尚缺乏对LU的精细化研究以及3种力学效应的对比关系的探讨；同时，3种力学效应对波浪规模的影响也尚未进行充分的探究。
基于以上背景，本文进一步解释了TS框架下3种力学效应及三者之间的关联与区别，并补充阐述体积抬升效应修正模型及其应用范围。在此基础上，选取滑坡涌浪物理模型试验中的若干工况，结合嵌入LU+DA+PA动力学模型的TS滑坡涌浪数值模拟方法，设计考虑3种力学效应参数比例变化的数值试验，研究PA、DA和LU&BT参数敏感性，分析PA、DA和LU 3种力学效应分别对涌浪波幅和波能的贡献，揭示滑坡涌浪产生过程中3种力学效应对波浪规模的影响。本研究旨在为滑坡涌浪灾害的预测与评价提供动力学模型及理论依据，对提高滑坡涌浪灾害预测准确性具有重要意义。
1 滑坡涌浪产生过程力学机制及模型 1.1 滑体与水体受力状态分析滑坡涌浪产生全过程主要包括滑体入水前运动、滑体与水体相互作用(冲击水体)、涌浪的传播与爬坡等几个关键阶段[28]。涌浪产生过程同时伴随波浪的反射、折射和色散等强非线性效应，为避免产生过多不确定性参数，且同时能使模型满足灾害评估的精度要求，可在确定主要力学参数及模型的基础上对滑坡涌浪的数值建模做适当假设。
针对滑坡入水之前的运动，考虑滑体受到重力的驱动产生加速度源，经过推导，流体及类流体在侧压力和重力的共同作用下，其加速度受到上表面梯度的控制[8,24]。同时引入基底摩擦(basal friction，BF)和动摩擦(dynamic friction on dry land，DF1)，即分别为滑体与地面的接触摩擦和滑体自身的黏滞阻力。二者产生的综合摩擦阻力如图 1所示，控制流态化或散粒体滑坡的运动和堆积。若滑体在运动过程中不产生解体，仍然保持刚性滑体的性质，则其加速度受下表面梯度的控制，且不考虑黏滞阻力。
随着滑体逐渐进入水体，滑体还受到水体对其产生的动摩擦阻力(dynamic friction under water，DF2)，相当于水体受到滑体冲击挤压和界面摩擦的反作用力，如图 1(c)所示。滑体入水后受到的阻力增加，速度逐渐减小、堆积、直至停止。
水体被滑体侵占而抬升，同时因受到摩擦阻力的反作用力和冲击挤压力而获得相应的动量，因此水体受到冲击而产生涌浪的力学模型可以表达为
| $ \boldsymbol{a}(\boldsymbol{r}, t)=\boldsymbol{a}_{\mathrm{gLU}}(\boldsymbol{r}, t)+\boldsymbol{a}_{\mathrm{DA}}(\boldsymbol{r}, t)+\boldsymbol{a}_{\mathrm{PA}}(\boldsymbol{r}, t) $ | (1) |
式中：agLU(r, t)为考虑流体侧压力与体积抬升效应的重力加速度项；对应于体积抬升效应(lift up，LU)，aDA(r, t)和aPA(r, t)分别表示由界面摩擦(drag along，DA)和冲击挤压(push ahead，PA)效应产生的将水平方向动量由滑体传递给水体的加速度项；vs和vw分别为滑体速度和水体速度，如图 1(a)所示。
|
图 1 滑体入水过程受力分析 Fig. 1 Sketch of the general forces during wave generation |
当滑体运动停止后，水体受到的合加速度为
| $ \boldsymbol{a}(\boldsymbol{r}, t)=\boldsymbol{a}_{\mathrm{g}}(\boldsymbol{r}, t)+\boldsymbol{a}_{\mathrm{d}}(\boldsymbol{r}, t) $ | (2) |
式中：ag(r, t)为考虑流体侧压力的重力加速度项，ad(r, t)是水体受到的除滑体以外的综合摩擦阻力(dynamic friction of water，DFW)产生的加速度项；包括水体自身的黏滞阻力以及与周围河道地形之间的摩擦阻力图 1(b)。由于水体的黏滞阻力很低，许多情况下被视为无黏性流，且涌浪在跨越江海的传播中能量耗散速率极低，一般可以设置ad(r, t)=0。但是，涌浪在爬坡过程中，由于水体处于一个高度非线性状态，且运动形态受地形影响较大，一般应该设置一定量的动摩擦阻力。
1.2 体积抬升效应及其修正体积抬升效应(lift up，LU)描述滑体在冲击水体过程中占据水体体积，同时使河底地形发生变化的力学过程。该过程中水体受到竖直方向的抬升而形成一定水体雍高，一般认为水体被抬升的总量等于滑体入水的总体积。LU不直接向水体传递水平方向动量，但是被抬升的水体在重力作用下会产生加速度从而驱动涌浪传播。传统的海啸模型一般仅采用LU产生的加速度agLU(r, t)来描述体积抬升效应中的动量传递过程。
| $ \boldsymbol{a}_{g \mathrm{LU}}(\boldsymbol{r}, t)=-g \Delta \eta(\boldsymbol{r}, t) $ | (3) |
式中η(r, t)为t时刻位置r处水体的表面高度。
对于库岸滑坡涌浪，通常情况下当滑体入侵水体后，水体表面高度按下式计算
| $ \eta(\boldsymbol{r}, t)=H_{\mathrm{w}}(\boldsymbol{r}, t)+T_0(\boldsymbol{r})+\min \left[T_{\mathrm{s}}(\boldsymbol{r}, t), H_{\mathrm{W}}(\boldsymbol{r}, t)\right] $ | (4) |
式中：Hw(r, t)为t时刻r处的水深，T0(r)为初始河床r处的地面高程，Ts(r, t)为t时刻r处滑体的厚度。
但是针对深海的海底滑坡，被抬升至海面水体的厚度会小于海底滑体的变化厚度，且被抬升至表面水体的覆盖面积要远大于初始滑体面积[29](图 2)。根据线性波理论，在水深为H的海洋中由一个微小的海底抬升uzbot(r0)引起的水面初始竖向位移uzsurf(r)可表达为
| $ \boldsymbol{u}_{\mathrm{z}}^{\text {surf }}(\boldsymbol{r})=\frac{1}{4 \mathtt{π}^2} \int\limits_{\boldsymbol{r}_0} \boldsymbol{u}_{\mathrm{z}}^{\text {bot }}\left(\boldsymbol{r}_0\right) \mathrm{d} \boldsymbol{r}_0 \int\limits_k \frac{\mathrm{e}^{\mathrm{i} \boldsymbol{k} \cdot\left(\boldsymbol{r}-\boldsymbol{r}_0\right)} \mathrm{d} \boldsymbol{k}}{\cosh (\boldsymbol{k} H)} $ | (5) |
|
图 2 海底滑坡涌浪产生过程机制 Fig. 2 Ketch map of submarine landslide tsunami generation mechanisms |
式中：
若底部抬升uzbot(r0)=uzbot(x0, y0)是一个在xy坐标系下厚度为U0且长宽为Wx和Wy的矩形，则式(5)变为
| $ \begin{aligned} \boldsymbol{u}_z^{\text {surf }}(x, y)= & \frac{U_0}{4 \mathtt{π}^2} \int\limits_{-W_x / 2}^{W_x / 2} \mathrm{~d} x_0 \int\limits_{-W_y / 2}^{W_y / 2} \mathrm{~d} y_0 \int\limits_{-\infty}^{\infty} \mathrm{d} \boldsymbol{k}_x \\ & \int\limits_{-\infty}^{\infty} \mathrm{d} \boldsymbol{k}_y \frac{\mathrm{e}^{i \boldsymbol{k}_x\left(x-x_0\right)} \mathrm{e}^{i \boldsymbol{k}_y\left(y-y_0\right)}}{\cosh (\boldsymbol{k }H )} \end{aligned} $ | (6) |
由于在一定程度上
| $ \frac{1}{\cosh (\boldsymbol{k} H)} \approx \frac{1}{\cosh \left(\boldsymbol{k}_x H\right)} \frac{1}{\cosh \left(\boldsymbol{k}_y H\right)} $ | (7) |
则式(6)变更为
| $ \begin{aligned} \boldsymbol{u}_{\mathrm{z}}^{\text {surf }}(x, y)= & \frac{U_0}{4 \mathtt{π}^2}\left[\int\limits_{-W_x / 2}^{W_x / 2} \mathrm{~d} x_0 \int\limits_{-\infty}^{\infty} \mathrm{d} \boldsymbol{k}_x \frac{\mathrm{e}^{\mathrm{i} \boldsymbol{k}_x\left(x-x_0\right)}}{\cosh \left(\boldsymbol{k}_x H\right)}\right] \\ & {\left[\int\limits_{-W_y / 2}^{W_y / 2} \mathrm{~d} y_0 \int\limits_{-\infty}^{\infty} \mathrm{d} \boldsymbol{k}_y \frac{\mathrm{e}^{\mathrm{i} \boldsymbol{k}_y\left(y-y_0\right)}}{\cosh \left(\boldsymbol{k}_y H\right)}\right] } \end{aligned} $ | (8) |
因此有
| $ \begin{gathered} \boldsymbol{u}_{\mathrm{z}}^{\text {surf }}(x, y)=\frac{4 U_0}{\mathtt{π}^2} \tan ^{-1} \xi\left(x, W_x, H\right) \tan ^{-1} \xi\left(y, W_y, H\right) \\ \xi(\boldsymbol{a}, W, H)=\frac{\sinh [\mathtt{π} W / 4 H]}{\cosh [\mathtt{π} \boldsymbol{a} / 2 H]} \end{gathered} $ | (9) |
式(9)代表经过平滑处理后的矩形底部抬升形成的涌浪位置和形态。Wx/H或Wy/H越小，则平滑效果越强，例如水深相对于矩形规模很大，则被平滑化的越多。在任何水深条件下，平滑处理均不改变被抬升水体的总体积，即
| $ U_0 W_x W_y=\int\limits_{\boldsymbol{r}} \boldsymbol{u}_{\mathrm{z}}^{\text {surf }}(\boldsymbol{r}) \mathrm{d} \boldsymbol{r} $ | (10) |
假定Wx=Wy=50 m，U0=1 m，水深H在5 m~50 m之间变动，跟据式(9)可计算出表面抬升高度和形态，如图 3所示(红色方框表示抬升范围)。H/W小于1/5时，抬升瞬间几乎没有产生扩散，表面的抬升非常接近于一个边长为W，高度为U0的矩形体，如图 3(a)所示。当水深逐渐增加，表面抬升范围扩散而抬升高度下降。水深H=W时，中间的被抬升高度仅为底面抬升量的0.2倍，如图 3(j)所示。
|
图 3 方形块体抬升1 m水位导致的水体表面形态变化 Fig. 3 Surface morphology changes caused by lifting a square block 1 m in water |
水底底面的变形会促使非动量传递(NMT)的水体抬升从而产生涌浪。在以往的TS涌浪产生的计算中考虑体积抬升效应时，无论是地震或海底滑坡产生的地底形变，均将水体表面的抬升量视为与水体底面的高度变化相等。然而这种假设在许多例子中不一定严密，因此在式(9)所适用的范围内，有必要解释底部和顶部运动的差异。
为了修正NMT方法，考虑图 3中红色的块体为TS中的一个微块体，边长为W。在时间t到t+dt内，块体i的底面高程的改变量用dU(xi, yi)表示。块体i表面抬升量将向周围块体扩散，当块体j不等于i时，则周围水块体上增加的厚度dH(xi, yi)为
| $ \begin{aligned} \mathrm{d} H\left(x_j, y_j\right)= & \frac{4 \mathrm{~d} U\left(x_i, y_i\right)}{\mathtt{π}^2} \tan ^{-1} \xi\left(x_j-x_i, W, H_i\right) \\ & \tan ^{-1} \xi\left(y_j-y_i, W, H_i\right), j \neq i \end{aligned} $ | (11) |
如果j块体等于i，则从该水块体中减去一个厚度增量dH(xi, yi)，即
| $ \begin{gathered} \mathrm{d} H\left(x_i, y_i\right)=-\mathrm{d} U\left(x_i, y_i\right) \\ {\left[1-\frac{4}{\mathtt{π}^2} \tan ^{-1} \xi\left(0, W, H_i\right) \tan ^{-1} \xi\left(0, W, H_i\right)\right], j=i} \end{gathered} $ | (12) |
式(12)是式(11)在网格单元尺度上的应用，既能满足水体厚度Hi在时间和空间上变化以描述被抬升水体实际表面形态，且保证被抬升水体总体积不变。
当河床发生变形时，式(11)、(12)需要应用于所有在t+dt时间内底部高程发生变化的方形块体上。如果H/W小于1/5，式(11)近似为0，式(12)近似为1，则该修正可以忽略，即表面抬升量直接用水底抬升量表示。当j块体与i块体的距离为数倍的水深时，则式(11)可以忽略(停止计算)，因为此时的修正量非常小。如果H/W很大，则计算量较大，需要增加一个计算距离的约束条件。
该方法是针对滑坡涌浪过程中被抬升水体底部至顶部运动(bottom to top, BT)形态的修正，是LU效应必要且可靠的补充，本文对该体积抬升效应的修正简称LU&BT。
1.3 界面摩擦效应界面摩擦效应(drag along，DA)描述滑体冲击水体过程中由滑体与水体接触界面拖拽摩擦并将动量传递给水体的力学过程。DA以剪切力的形式存在于滑体与水体的接触界面上。滑体与水体受到的界面摩擦力大小相同方向相反，即DA阻碍滑体运动的同时推动水体向前运动。DA在网格单元中的表达式为
| $ \begin{gathered} \boldsymbol{a}_{\mathrm{DA}}(\boldsymbol{r}, t)=c_{\mathrm{DA}}\left[\left(\boldsymbol{v}_{\mathrm{s}}(\boldsymbol{r}, t)-\boldsymbol{v}_{\mathrm{w}}(\boldsymbol{r}, t)\right) \cdot\right. \\ \left.\hat{\boldsymbol{v}}_{\mathrm{s}}(\boldsymbol{r}, t)\right] \frac{\boldsymbol{v}_{\mathrm{s}}(\boldsymbol{r}, t)}{H_{\mathrm{w}}(\boldsymbol{r}, t)} \\ \end{gathered} $ | (13) |
| $ \begin{array}{c} \hat{\boldsymbol{v}}_{\mathrm{s}}(\boldsymbol{r}, t)=\frac{\boldsymbol{v}_{\mathrm{s}}(\boldsymbol{r}, t)}{\left|\boldsymbol{v}_{\mathrm{s}}(\boldsymbol{r}, t)\right|}= \\ \frac{\left[\boldsymbol{v}_{\mathrm{s}}^x(\boldsymbol{r}, t) \hat{\boldsymbol{x}}+\boldsymbol{v}_{\mathrm{s}}^y(\boldsymbol{r}, t) \hat{\boldsymbol{y}}\right]}{\left[\boldsymbol{\boldsymbol { v }}_{\mathrm{s}}^x(\boldsymbol{r}, t) \boldsymbol{v}_{\mathrm{s}}^x(\boldsymbol{r}, t)+\boldsymbol{v}_{\mathrm{s}}^y(\boldsymbol{r}, t) \boldsymbol{v}_{\mathrm{s}}^y(\boldsymbol{r}, t)\right]^{1 / 2}}= \\ \hat{\boldsymbol{v}}_{\mathrm{s}}^x(\boldsymbol{r}, t) \hat{\boldsymbol{x}}+\hat{\boldsymbol{v}}_{\mathrm{s}}^y(\boldsymbol{r}, t) \hat{\boldsymbol{y}} \end{array} $ | (14) |
式中：vs(r, t)和vw(r, t)分别为滑体和水体的速度；Hw(r, t)为网格单元中水体的厚度；
因此，表面光滑和粗糙滑体的界面摩擦系数显然不同，导致水体受到的DA加速度有差异，即稍大的cDA值会产生大的DA加速度进而形成较大的首浪。DA加速度与典型的相互作用力类似，可以降低滑体的速度且加速水体的向前运动。DA与其他只作用在流体反方向的力不同，它能够在滑体运动的任意方向上加速水体运动，将滑体动量传递给水体以形成涌浪。
1.4 冲击挤压效应文献[8]提出了LU+DA产生滑坡涌浪的力学机制，并在三峡库区龚家方滑坡涌浪案例中应用TS进行了较好的反演模拟。在后期进一步的研究中发现，对于厚度/水深较小的滑体而言，仅考虑LU+DA两种力学机制，也能在数值模拟中得到一个接近的实测数据的浪高值；然而对于厚度相对水深较大的厚层滑体诱发的涌浪而言，仅考虑以上两种力学效应不足以准确描述涌浪的产生，数值模拟结果也会产生较大偏差，例如三峡库区的千将坪滑坡涌浪案例，意大利的瓦伊昂水库滑坡案例。
鉴于此，文献[26]基于TS数值模拟框架，补充提出了滑体-水体相互作用中的冲击挤压效应(push ahead，PA)，描述了滑体冲击水体过程中由滑体竖向表面冲击挤压并将动量传递给水体的力学过程。滑体对水体的冲击挤压效应与流场中压差阻力(pressure drag friction)形成的机制相同，是滑体前缘迎水面推动水体时产生的动量传递，该过程由滑体对水体产生了一个水平方向上的加速度，其在单位网格上的表达式为
| $ \begin{gathered} \boldsymbol{a}_{\mathrm{PA}}(\boldsymbol{r}, t)=c_{\mathrm{PA}}\left|\hat{\boldsymbol{v}}_{\mathrm{s}}(\boldsymbol{r}, t) g \nabla T_{\mathrm{s}}(\boldsymbol{r}, t)\right| \\ {\left[\left(\boldsymbol{v}_{\mathrm{s}}(\boldsymbol{r}, t)-\boldsymbol{v}_{\mathrm{w}}(\boldsymbol{r}, t)\right) \cdot \hat{\boldsymbol{v}}_{\mathrm{s}}(\boldsymbol{r}, t)\right] \frac{\boldsymbol{v}_{\mathrm{s}}(\boldsymbol{r}, t)}{H_{\mathrm{w}}(\boldsymbol{r}, t)}} \end{gathered} $ | (15) |
| $ \left|\hat{\boldsymbol{v}}_{\mathrm{s}}(\boldsymbol{r}, t) g \nabla T_{\mathrm{s}}(\boldsymbol{r}, t)\right|\left|\boldsymbol{v}_{\mathrm{s}}(\boldsymbol{r}, t)\right|=\left|\mathrm{d} T_{\mathrm{s}}(\boldsymbol{r}, t)\right| / \Delta t $ | (16) |
| $ \hat{\boldsymbol{v}}_{\mathrm{s}}(\boldsymbol{r}, t)=\frac{\boldsymbol{v}_{\mathrm{s}}(\boldsymbol{r}, t)}{\left|\boldsymbol{v}_{\mathrm{s}}(\boldsymbol{r}, t)\right|} $ | (17) |
式中：dTs(r, t)为滑体在位置r和时间间隔Δt内的厚度变化；
基于流体动力学数值方法Tsunami Squares(TS)，将上述滑坡涌浪产生的动力学效应及参数进行集成，建立考虑滑体-水体相互作用力学效应的滑坡涌浪动力学模型。该方法数值过程采用一种免解流体控制方程的思路，通过网格移动与分配严格遵守质量和动量守恒的第一性原理。计算流体运动过程只需重点关注微元网格平均加速度的来源，扩展不同的外部加速度(来源于不同驱动力包括重力、风、气压、黏滞阻力或边界摩擦力等)，即可描述各种单相或多相流动力学过程。本文前述详细描述了滑坡涌浪产生过程的PA、DA和LU&BT力学效应，提出了各力学效应在网格单元尺度的计算公式。因此，基于TS数值方法能够无缝集成考虑滑体-水体相互作用力学效应的滑坡涌浪动力学模型。
2 物理模型试验及TS数值模拟反演 2.1 物理模型试验设置本文设计矩形水槽试验装置，结合TS数值模拟，开展滑坡涌浪产生过程中的PA、DA和LU&BT效应研究。选取典型河道滑坡涌浪地质原型以1∶ 200相似比建立试验模型，包括矩形河道模型、控制滑槽和散粒体滑坡模型。河道模型长8.0 m、宽4.5 m、高0.75 m、两侧坡面倾角α=45°，河道斜坡面设置10 cm水位标志线，底部设置10 cm×10 cm网格线以方便试验观测。滑槽模型长2.0 m，宽0.5 m，高0.5 m，贴河道斜坡面固定，能够调节滑体下滑高度和距离，保证滑体达到设计的入水速度要求。
如图 4(a)所示，试验采用天然河流鹅卵石作为散粒体滑坡相似材料，粒径为1.0~5.0 cm，且滑体初始形态设置为临空面竖直的三角体以方便TS数值建模。试验中设计三角体滑坡模型，直角边长(EL)0.25 m、0.35 m和0.5 m，分别对应滑坡体积为0.016 m3(K1)、0.031 m3(K2)和0.063 m3(K3)。通过设置滑体距离静水面高度控制散粒体滑坡入水速度，分别为22.6、53.0、81.6 cm，对应编号为V1、V2、V3。试验设置25 cm(W1)和50 cm(W2)两种静水位，每一种工况的编号均由对应信息的符号组合而成。
|
图 4 滑坡涌浪模型试验水槽及其附属设备 Fig. 4 Photos of the rectangle flume apparatus and its auxiliary facilities |
试验过程中采用两部高清摄像机、8个搭载SAD1000传感器的波高仪及水槽网格，监测滑坡涌浪试验过程信息，包括滑体运动特征、滑体堆积特征(包括滑体堆积范围、堆积形状、堆积厚度等)，以及涌浪形态特征(包括波形、波高、波速等)，如图 4(b)所示。
2.2 物理模型试验及数值反演散粒体以速度V2从斜面一定高度下滑至水深25 cm的水槽体中，观测试验中首浪的产生与滑体堆积形态。以试验中的K2V2W1工况为例，将观测到的滑体堆积厚度与范围、下滑过程中每一时刻滑体的位置与滑动距离作为数值模型的强约束条件[27]，利用集成PA、DA和LU&BT的TS方法建立物理试验的数值模型，设置参数CPA=0.4，CDA=0.03，BF=0.2，DF1=0.3，DF2=0.35，反演分析物理试验中滑坡和涌浪的动力学过程与特征。
滑体冲击产生首浪的瞬态动力学过程如图 5所示，图 5(a)~5(c)为滑坡动力学过程；图 5(d)~5(f)为涌浪动力学特征。数值模拟结果矩形框内标识时间和涌浪能量，其中，KE为动能，D-GPE为重力势能，TOT为总能量。滑体冲击挤压水体，将水体抬升至一定高度形成首浪，静水面被冲击后在滑体正前方瞬时形成一个弧形波峰，同时伴随产生大量的气泡和浪花，滑体与水体充分接触且通过复杂相互作用交换大部分能量，首浪能量开始通过波动形式向整个水槽扩散传播。随后，受水槽地形限制，次生波峰和波谷相互叠加、反射、折射形成更加复杂的波浪直至能量完全耗散。此外，滑体在自重作用下沿滑槽加速下滑并与底部发生碰撞后，由于底部摩擦而缓慢堆积，滑体厚度呈现中间厚两边薄的特征。
|
图 5 K2V2W1工况下TS模拟与物理试验结果对比 Fig. 5 Physical test and TS simulation results for the K2V2W1 experiment |
对比该工况下涌浪波形的TS模拟与观测结果(图 6)可知：在波高仪WG1处的最大波峰达到42 mm，该波峰与高清摄像机捕捉到的形态特征基本一致，且波幅随着观测距离的拉远迅速降低至WG8处的6 mm。TS模拟的涌浪波形与试验观测结果在前20 s匹配度较高，但由于涌浪受边界条件影响产生复杂的叠加、反射/折射、破碎等现象的影响，两者波形后续呈现出轻微的离散、错位和波动现象。综上，首浪波幅的迅速降低表明当滑体速度降为零或者没有滑体-水体相互作用时，PA和DA效应则会消失。波形曲线在前20 s匹配度较高表明TS数值模拟与试验观测结果吻合良好，验证了滑坡涌浪产生过程动力学模型的可靠性。
|
图 6 K2V2W1工况下数值模拟与观测涌浪波形对比结果 Fig. 6 Comparison results between numerical simulation and observed surge waveforms in K2V2W1 condition |
经前文对于体积抬升效应及其修正的解析表达可知，LU&BT效应受滑体厚度/水深的影响，为了定量化研究滑体冲击水体产生涌浪过程中的LU&BT效应，本文分别选取厚度/水深最小和最大的两组工况K2V2W2和K3V3W1，应用考虑和不考虑BT修正的TS方法模拟滑坡涌浪。图 7(a)显示考虑BT修正的模拟结果在WG1和WG2处的首浪波幅较不考虑BT修正的分别减小45%(WG1约为35 mm)和35%(WG2约为15 mm)。相比之下，图 7(b)试验工况K3V3W1由有无BT修正产生的模拟结果差异为33%(约为35 mm)且只发生在WG1位置。研究表明，BT修正在深水环境比浅水环境对体积抬升效应的影响更强，即使在滑坡体积和冲击速度较小的情况也同样适用。此外，模拟和观测的涌浪波幅在远离滑源区基本不产生差异，表明BT修正在该区域无影响。实际上，BT修正作用类似一个过滤器，它降低首浪波峰的高度但是不改变被抬升水体的总体积。该结论也从侧面反映出在涌浪灾害强度的大小不能单凭最大初始涌浪高度去评估，而要综合考虑滑体的规模和水深情况。
|
图 7 有无BT修正的TS方法模拟结果对比 Fig. 7 The comparison of TS method simulation results with and without BT correction |
选取速度水深相同，而体积不同的两组工况K2V2W1和K3V2W1，模拟对比了两种滑体-水体相互作用工况：LU&BT+PA+DA；LU&BT。由图 8(a)可知，LU&BT+PA+DA和LU&BT两种工况在WG1处产生涌浪波幅的最大差距可达35.7%(15 mm)，随着远离滑源区，涌浪波幅差距逐渐降低33.3%(WG2，6 mm)至10%(WG8，0.6 mm)。在图 8(b)中，由于滑坡体积的增大导致观测到更大的涌浪波幅差距，但是，最大波幅降低百分数从38.7%(WG1，24 mm)至8.4%(WG8，1 mm)，与前者比较接近。因此，结果表明滑体冲击水体过程中PA和DA的作用至关重要，忽略它们的贡献将导致预测涌浪波幅降低约38.7%。
|
图 8 LU&BT+PA+DA与LU&BT的TS方法模拟结果对比 Fig. 8 Comparison of simulation results between LU&BT+PA+DA and LU&BT based TS methods |
图 9和图 10呈现了试验工况K2V2W1和K3V2W1对应的考虑不同力学效应组合下TS模拟的涌浪波形。图 9(a)和图 10(a)对比分析了两种力学效应工况组合下的滑坡涌浪波形：LU&BT+PA+DA；LU&BT+PA。由图可知，试验工况K2V2W1在两种组合下存在一定的涌浪波幅差异，力学效应组合“LU&BT+PA+DA”产生的涌浪波幅较组合“LU&BT+PA”稍大，且从在WG1处高出4 mm(9.5%)降低至WG8处的0.1 mm(1.7%)。同时，试验工况K3V2W1在两种组合下存在的差异与工况K2V2W1基本一致，从WG1处的6 mm(9.7%)降低至WG8处的0.2 mm(1.7%)。以上两组力学效应对比工况模拟结果表明：界面摩擦(DA)效应对滑坡涌浪的产生具有一定的作用，DA效应对首浪波高的贡献占比达到9.7%，随着传播距离的增大，DA影响率将逐渐减小。
|
图 9 K2V2W1工况下验证DA与PA对涌浪波形的影响 Fig. 9 Verification of DA and PA effects on surge wave form under K2V2W1 condition |
|
图 10 K3V2W1工况下验证DA与PA对涌浪波形的影响 Fig. 10 Verification of DA and PA effects on surge waveform under K3V2W1 condition |
图 9(b)和图 10(b)对比分析了另外两种力学效应工况组合下的滑坡涌浪波形：3)LU&BT+PA；LU&BT。由对比结果可知，LU&BT+PA和LU&BT的数值分析得到与LU&BT+PA+DA和LU&BT+PA相似的涌浪波幅差异化的趋势，但是波幅相差的绝对值要比后者大。由图可知，试验工况K2V2W1在两种组合下的波幅差距由WG1的11 mm(26.2%)降低至WG8处的0.5 mm(8.3%)。试验工况K3V2W1在两种组合下的波幅差距由WG1处的18 mm(29%)降低至WG8处的0.8 mm(6.7%)。以上两组力学效应对比工况模拟结果表明：冲击挤压(PA)效应对滑坡涌浪的产生具有重要作用，PA效应对首浪波高的贡献达到29%，随着传播距离的增大，PA影响率将逐渐减小。
综上，滑坡冲击水体产生涌浪的过程中，冲击挤压(PA)和界面摩擦(DA)都产生重要的作用，即在滑体-水体相互作用影响范围内拖拽、冲击和挤压水体以形成涌浪，忽略二者的贡献将导致预测涌浪波幅偏低，同时，PA较DA对涌浪规模的贡献作用更大。
3.3 PA、DA和LU&BT对涌浪能量的贡献滑坡涌浪产生过程中滑体和水体之间的能量转化一致被认为是个异常复杂的动力学过程。在TS数值方法中，滑坡的动能通过PA、DA和LU&BT转化成水体的动能和势能，总能量为动能与势能之和。为了研究PA和DA在该过程中的作用，本文基于TS方法计算了两者对产生涌浪波能的贡献比例。图 11显示了试验工况K2V2W1和K3V2W1下TS模拟的涌浪势能、动能和总能量。由图 11可知，TS计算的涌浪波能时序曲线基本遵循涌浪传播过程中势能和动能相互转换的基本规律。涌浪在受到滑体冲击后波能快速增加形成一个波能峰值后随着传播距离逐渐降低至耗散。
|
图 11 试验工况K2V2W1和K3V2W1对应TS模拟的水槽涌浪势能、动能和总能量变化时序曲线 Fig. 11 Potential, kinetic, and total energy forms summed over the whole tank area for simulated K2V2W1 and K3V2W1 experiments |
由波能曲线可知，冲击挤压(PA)作用大致贡献给涌浪23%~30%的势能、28%~29%的动能和25%~30%的总能量。相比之下，界面摩擦(DA)作用大致贡献给涌浪17%~25%的势能、14%~16%的动能和16%~23%的总能量。PA和DA共同贡献给涌浪约48%的势能、43%的动能和47%的总能量。同时，通过PA作用的滑体-水体能量转化较DA和LU快速，特别是在滑移时间短的高速运动滑体冲击浅水水体环境时体现更显著。研究表明，对比试验工况K3V2W1和K2V2W1值，前者较大的滑体体积导致PA和LU&BT效应更显著，同时降低了DA效应的相对能量贡献率。因此，PA在试验工况K3V2W1的能量贡献率较K2V2W1大，而DA在试验工况K3V2W1的贡献率较K2V2W1小。总之，PA、DA和LU&BT对滑坡涌浪产生过程中能量转化起到决定性的重要作用，任何建模过程都不可忽视。
4 讨论 4.1 3种力学效应的相互关系由3种力学效应的模型表达式可以看出，LU是基于滑体厚度的空间导数对水体加速，而PA是基于滑体厚度的时间导数对水体加速。虽然“体积抬升(LU)”和“冲击挤压(PA)”最终可能会起到同样的作用，但PA可以更快地达到效果。在浅水环境中，物体的快速移动往往伴随着短暂的滑动。对于这种情境，PA显得尤为重要。
对比PA和DA，二者的表达式在形式上非常相似，即它们都由滑体运动方向上的滑体和水体速度矢量差
PA和DA的相同点在于：1)PA和DA都在滑体的运动方向上作用；2)PA和DA在滑体速度为零的时候都会消失；3)当水体被加速到与滑体速度相等时，PA和DA效应也都消失；4)当水体速度为零vw(r, t)=0时，PA和DA加速度都会降低至|vs(r, t)|vs(r, t)形式；5)当水体厚度较小且滑体速度相比水体速度较大时，即高速滑体冲击浅水水体时，PA和DA效应都非常显著。因此，不考虑PA和DA的模型可能低估滑坡产生涌浪的规模。
PA和DA的区别在于：PA表达式多一个
LU&BT是一种非动量传递法，即将每一个水块体中的水垂直抬升或降低，重力作用扰动水面使水体产生流动或波浪，但水自身的抬升并没有动量的转换。在滑坡涌浪领域通常将非动量传递与水阻力的影响综合考虑。根据LU&BT修正公式可知，如果滑体厚度与水深之比小于0.2，则LU&BT修正可以忽略。实际上，BT修正作用类似一个过滤器，它降低首浪波峰的高度，对抬升水体的形态重分布，但是不改变被抬升水体的总体积。该结论也从侧面反映出在涌浪灾害强度的大小不能单凭最大初始涌浪高度去评估，而要综合考虑滑体的规模和水深情况。
需要指出的是，虽然引入LU&BT修正模型可以使计算结果更趋于准确，但同时也会消耗更多的计算时长，而且修正模型对后续传播浪和爬坡浪高度和能量影响很小，因此在数值模拟计算时可根据不同的模拟对象、应用场景、计算目标和实际需求考虑是否引入。
5 结论本文基于前期滑坡涌浪产生过程研究提出的三种力学效应，进一步解释其在TS框架下的物理意义、解析表达和三者之间的关联与区别，并补充阐述体积抬升效应修正模型及其应用范围，通过物理模型试验和TS数值反演分析，揭示滑坡涌浪产生过程动力学特征以及对首浪生成规模的影响规律。论文的主要结论包括：
1) 体积抬升效应的BT修正能够较好模拟首浪形态，深水环境比浅水环境影响更显著；当水深与块体边长的比值(H/W)小于1/5时，可以忽略体积抬升效应的修正；相反，当H/W比值较大时，计算变得复杂，需要增加一个计算距离的约束条件。
2) 从力学模型表达可知，LU是基于滑体厚度的空间导数对水体加速，而PA是基于滑体厚度的时间导数对水体加速, 虽然二者会起到同样的作用，但PA会更快的达到效果；PA与DA的表达式比较相似，对于高速滑体冲击浅水水体的情境，PA和DA效应都很显著；此外，厚度变化率较大的滑体会产生较大的PA加速度，形态复杂且表面粗糙的滑体可能产生较大的DA加速度。
3) PA与DA效应对滑坡涌浪的产生具有重要的贡献作用，PA和DA对首浪波高的贡献率分别达到29%和9%，且均随着传播距离的增大而降低。
4) TS计算的涌浪波能时序曲线基本遵循涌浪传播过程中势能和动能相互转换的基本规律。PA贡献给涌浪的势能约23%~30%，动能约28%~29%，总能量约25%~30%；DA贡献给涌浪的势能约17%~25%，动能约14%~16%，总能量约16%~ 23%。同时，通过PA作用的滑体-水体能量转化比DA和LU更快速，特别是在滑移时间短的高速运动滑体冲击浅水水体环境时体现更显著。
本文进一步解释了滑坡涌浪产生过程的动力学效应及其解析模型，揭示了滑坡涌浪动力学机理，探讨了PA、DA和LU&BT对首浪规模的影响规律，研究成果对滑坡涌浪防灾减灾研究具有重要科学意义和应用价值。
| [1] |
WARD S N. Slip-sliding away[J]. Nature, 2002, 415(6875): 973. DOI:10.1038/415973a |
| [2] |
HUANG B, YIN Y, LI R, et al. Three-dimensional experimental investigation on hazard reduction of landslide-generated impulse waves in the Baihetan Reservoir, China[J]. Landslides, 2023, 20(9): 2017. DOI:10.1007/s10346-023-02068-w |
| [3] |
胡杰, 王道熊, 胡斌. 库岸滑坡灾害及其涌浪分析[J]. 华东交通大学学报, 2003, 20(5): 26. HU Jie, WANG Daoxiong, HU Bin. Analysis on reservoir landslides and its surge[J]. Journal of East China Jiaotong University, 2003, 20(5): 26. |
| [4] |
汪洋, 殷坤龙. 水库库岸滑坡涌浪的传播与爬高研究[J]. 岩土力学, 2008, 29(4): 1031. WANG Yang, YIN Kunlong. Research on propagation and climb height of surge triggered by landslide in reservoir[J]. Rock and Soil Mechanics, 2008, 29(4): 1031. DOI:10.16285/j.rsm.2008.04.003 |
| [5] |
FRITZ H, MOHAMMED F. Lituya Bay landslide impact generated mega-tsunami 50th anniversary[J]. Pure Appl geophys, 2009, 166(1): 153. DOI:10.1007/s00024-008-0435-4 |
| [6] |
WARD S N, SIMON N D. The 1958 Lituya Bay landslide and tsunami-a tsunami ball approach[J]. Journal of Earthquake and Tsunami, 2010, 4(4): 285. DOI:10.1142/S1793431110000893 |
| [7] |
WANG J, WARD S N, XIAO L. Numerical simulation of the december 4, 2007, landslide-generated tsunami in Chehalis Lake, Canada[J]. Geophysical Journal International, 2015, 201(1): 372. DOI:10.1093/gji/ggv026 |
| [8] |
XIAO L, WARD S N, WANG J. Tsunami squares approach to landslide-generated waves: application to Gongjiafang Land-slide, Three Gorges Reservoir, China[J]. Pure and Applied Geophysics, 2015, 172(12): 3639. DOI:10.1007/s00024-015-1045-6 |
| [9] |
XIAO L, WANG J, WARD S N, et al. Numerical modeling of the June 24, 2015, Hongyanzi Landslide generated impulse waves in Three Gorges Reservoir, China[J]. Landslides, 2018, 15(12): 2385. DOI:10.1007/s10346-018-1057-2 |
| [10] |
DAY S, LANES P, SILVER E, et al. Submarine landslide deposits of the historical lateral collapse of Ritter Island, Papua-New Guinea[J]. Marine And Petroleum Geology, 2015, 67(1): 419. DOI:10.1016/j.marpetgeo.2015.05.017 |
| [11] |
MARAMAI A, GRAZIANI L, ALESSIO G, et al. Near-and far-field survey report of the 30 December 2002 Stromboli(Southern Italy) Tsunami[J]. Marine Geology, 2005, 215(1/2): 93. DOI:10.1016/j.margeo.2004.11.009 |
| [12] |
CASALBORE D, BOSMAN A, CHIOCCI F L. Study of recent small-scale landslides in geologically active marine areas through repeated multibeam surveys: examples from the Southern Italy[J]. Springer Netherlands, 2012, 31(1): 573. DOI:10.1007/978-94-007-2162-3_51 |
| [13] |
PAKOKSUNG K, SUPPASRⅡ A, IMAMURA F, et al. Simulation of the submarine landslide tsunami on 28 September 2018 in Palu Bay, Sulawesi island, Indonesia, using a two-layer model[J]. Pure and Applied Geophysics, 2019, 176(8): 3323. DOI:10.1007/s00024-019-02235-y |
| [14] |
ZHAO L, MAO J, BAI X, et al. Finite element simulation of impulse wave generated by landslides using a three-phase model and the conservative level set method[J]. Landslides, 2016, 13(1): 85. DOI:10.1007/s10346-014-0552-3 |
| [15] |
KARAHAN M, ERSOY H, AKGUN A. A 3D numerical simulation-based methodology for assessment of landslide-generated impulse waves: a case study of the Tersun Dam reservoir(NE Turkey)[J]. Landslides, 2020, 17(12): 2777. DOI:10.1007/s10346-020-01440-4 |
| [16] |
黄波林, 殷跃平, 李滨, 等. 柱状危岩体坐落式崩塌产生涌浪的简化数值模型与校验[J]. 岩土力学, 2021, 42(8): 2269. HUANG Bolin, YIN Yueping, LI Bin, et al. Simplified numerical model and verification for the impulse wave generated by situated collapse of a dangerous columnar rock mass[J]. Rock and Soil Mechanics, 2021, 42(8): 2269. DOI:10.16285/j.rsm.2020.1639 |
| [17] |
HUANG T, ZHANG H, SHI Y. Numerical simulation of landslide-generated tsunamis in lakes: a case study of the Xiluodu reservoir[J]. Science China Earth Sciences, 2023, 66(2): 393. DOI:10.1007/s11430-022-9989-1 |
| [18] |
荆海晓, 陈国鼎, 李国栋. 水下滑坡产生涌浪波特性的数值模拟研究[J]. 应用力学学报, 2018, 35(3): 503. JING Haixiao, CHEN Guoding, LI Guodong. Numerical study on characteristics of water waves generated by submerged landslide[J]. Chinese Journal of Applied Mechanics, 2018, 35(3): 503. |
| [19] |
司鹏飞, 田利勇. 基于流固耦合方法的滑坡涌浪数值模拟[J]. 水利水电科技进展, 2020, 40(6): 28. SI Pengfei, TIAN Liyong. Numerical simulation of landslide-induced impulse waves based on fluid-solid coupling method[J]. Advances in Science and Technology of Water Resources, 2020, 40(6): 28. |
| [20] |
薛宏程, 马倩, 彭杏瑶, 等. 基于Herschel-Bulkley流变模型的滑坡涌浪数值模拟方法研究[J]. 水利学报, 2023, 54(3): 268. XUE Hongcheng, MA Qian, PENG Xingyao, et al. Numerical simulation of landslide generated impulse waves based on Herschel-Bulkley rheological model[J]. Journal of Hydraulic Engineering, 2023, 54(3): 268. |
| [21] |
王志超, 李大鸣. 基于SPH-DEM流-固耦合算法的滑坡涌浪模拟[J]. 岩土力学, 2017, 38(4): 1226. WANG Zhichao, LI Daming. Fluid-structure coupling algorithm based on SPH-DEM and application to simulate landslide surge[J]. Rock and Soil Mechanics, 2017, 38(4): 1226. DOI:10.16285/j.rsm.2017.04.038 |
| [22] |
李承德. 基于耦合FDEM-CEL方法的可破裂岩体滑坡涌浪数值模拟研究[D]. 南昌: 南昌大学, 2023 LI Chengde. Numerical simulation of landslide surge in fracturable rockmass based on coupled FDEM-CEL method[D]. Nanchang: Nanchang University, 2023 |
| [23] |
肖莉丽. 库岸滑坡涌浪数值模拟研究[D]. 武汉: 中国地质大学, 2016 XIAO Lili. Numerical simulation of landslide generated waves in reservoir[D]. Wuhan: China University of Geosciences, 2016 |
| [24] |
汪洋. 水库库岸滑坡速度及其涌浪灾害研究[D]. 武汉: 中国地质大学, 2007 WANG Yang. The research on speed of the landslide and its surge hazard in reservoir[D]. Wuhan: China University of Geosciences, 2007 |
| [25] |
WANG J, XIAO L, WARD S N. Tsunami Squares modeling of landslide tsunami generation considering the 'Push Ahead' effects in slide/water interactions: theory, experimental validation, and sensitivity analyses[J]. Engineering Geology, 2021, 288: 106141. DOI:10.1016/j.enggeo.2021.106141 |
| [26] |
WANG J, XIAO L, WARD S N, et al. Tsunami Squares modeling of the 2007 Dayantang landslide generated waves considering the effects in slide/water interactions[J]. Engineering Geology, 2021, 284(1/2): 106032. |
| [27] |
肖莉丽, 王佳佳, 李枝强, 等. 考虑滑体-水体相互作用的滑坡涌浪产生过程动力学模型研究[J]. 岩石力学与工程学报, 2022, 41(12): 2404. XIAO Lili, WANG Jiajia, LI Zhiqiang, et al. Research on dynamic models of landslide tsunami generation considering/water interactions[J]. Chinese Journal of Rock Mechanics and Engineering, 2022, 41(12): 2404. DOI:10.13722/j.cnki.jr-me.2021.1244 |
| [28] |
WANG J, WARD S, XIAO L. Tsunami Squares modeling of landslide generated impulsive waves and its application to the 1792 Unzen-Mayuyama Mega-slide in Japan[J]. Engineering Geology, 2019, 256: 121. DOI:10.1016/j.enggeo.2019.04.020 |
| [29] |
KELFOUN K, GIACHETTI T. Landslide-generated tsunamis at Réunion Island[J]. Journal of Geophysical Research Atmosphere, 2010, 115(F4): 12. |