Numerical Simulations of Tsunami Wave Generation by Submarine Landslides: Validation and Sensitivity Analysis to Landslide Parameters

Landslides are the second most frequent tsunami source worldwide and caused several destructive events such as the catastrophe of July 1998 in Papua New Guinea, which resulted in more than 2,200 deaths (Synolakis et al. 2002; Satake and Tanioka 2003; Lynett et al. 2003; Synolakis 2003; Tappin et al. 2014; Harbitz et al. 2014; Grilli et al. 2017). Recent large landslide tsunamis, such as the 2018 Palu (Indonesia) and the 2018 Anak Krakatau events, have highlighted the devastating social and economic impacts of landslide tsunamis (Takagi et al. 2019; Grilli et al. 2019). Unlike tectonic tsunamis whose generation, propagation, and coastal amplification are well understood and modeled, the existing knowledge on landslide tsunamis is relatively limited.

Researches on landslide tsunamis have been focused on experimental, analytical, and numerical studies. Analytical solutions are generally available for only simple cases and thus unable to account for more realistic landslide events (Lo and Liu 2017). Laboratory physical modeling is a straightforward way to study landslide tsunamis, and the result of experimental data can be used to validate numerical models (Watts et al. 2003; Grilli and Watts 2005; Liu et al. 2005; Enet and Grilli 2007; Fritz et al. 2009). However, laboratory works are expensive and time-consuming. By contrast, validated numerical simulations, in principle, allow for replication of actual events at a reasonable time and cost in addition to the high flexibility of the numerical models relative to laboratory studies (Lynett and Liu 2005; Liu et al. 2005; Løvholt et al. 2005; Horrillo et al. 2013; Grilli et al. 2017; Heidarzadeh et al. 2020).

Existing numerical models for landslide tsunamis have benefited from various forms of mathematical formulations; examples are shallow-water theory (Harbitz 1992), fully nonlinear potential flow (Grilli and Watts 2005; Grilli et al. 2002), Boussinesq equations (Lynett and Liu 2002; Watts et al. 2003; Fuhrman and Madsen 2009; Zhou and Teng 2010), nonhydrostatic wave equations (Ma et al. 2012) and Navier–Stokes equations (Heinrich 1992; Liu et al. 2005; Yuk et al. 2006; Abadie et al. 2010; Montagna et al. 2011; Rabinovich et al. 1999). Heinrich (1992) developed a shallow-water numerical model to clarify the efficiency of deep-water slumps in generating the tsunami. Monaghan and Kos (2000) used a 2D numerical particle model using the Smoothed Particle Hydrodynamics to clarify the wave formation and the boxlike landslide dynamics. Lynett and Liu (2002) developed a number of dimensionless relationships (e.g., location of secondary runup peak due to edge waves) by analyzing the runup of submerged landslides. Grilli and Watts (2005) presented empirical equations for landslide travel distance and initial wave amplitude and wave surface. Liu et al. (2005) found that, for the submerged cases of the solid sliding mass, the runup decreases as the submergence increases asymptotically and approaching zero as the submergence tends to infinity. Grilli et al. (2017) studied the effects of material rheology on the evolution of landslide tsunamis. Ruffini et al. (2019) presented a numerical model based on the nonhydrostatic shallow water equations to quantify the effect of the water body geometry on far-field landslide tsunami propagation.

In this work, we apply the Computational Fluid Dynamics (CFD) package FLOW3D-Hydro (version 1.0, update 1) (Flow Science 2021) to simulate and validate waves generated by solid-block submarine landslides with the aim of studying the sensitivity of the model to various landslide parameters. FLOW3D has been used to study landslide-generated waves in the past (Choi et al. 2007; Parambath 2010; Kim 2012; Horrillo et al. 2013; Kim et al. 2020). A motivation for this study is the lack of sufficient information about the effects of landslide parameters on the dynamics of the waves; specifically, limited information is published on two parameters of travel distance (S) and submergence depth (d). In particular, there are a few equations relating maximum initial landslide wave amplitude with S and d (Ataie-Ashtiani and Najafi-Jilani 2008; Takabatake et al. 2020). The innovation of this research is a comprehensive sensitivity analysis of landslide parameters as well as the validation of the model using physical experiments that would make the model suitable for future applications. We compare the results against previous experimental data to validate the model. The experimental data used for comparison are those described in Sabeti and Heidarzadeh (2021) and have been obtained through physical laboratory experiments. Sensitivity analyses of three landslide parameters, namely slide volume (V), S, and d, on the absolute value of the maximum initial landslide amplitude (η_max) are performed (Fig. 1) and the wave velocity field around the source region are investigated. It is noted that the maximum initial negative amplitude is usually larger than the positive one for submarine landslide tsunamis and thus here we target the absolute value of the maximum initial negative amplitude, η_max (Fig. 1). This validated landslide numerical model can be used for various studies on landslide-generated waves, including developing predictive equations for landslide parameters (e.g., η_max, S). As experimental works are costly and time-consuming, validated landslide models such as the one developed in this study can be a useful tool for engineers and scientists working on tackling landslide tsunami hazards.

The numerical simulations in this work have been performed using the CFD package FLOW3D-Hydro (version 1.0, update 1), which has been widely used in industry and academia for modeling Fluid Mechanics problems. FLOW3D-Hydro was applied to simulate fluid flow algorithms at the Los Alamos National Laboratory in the 1960s and 1970s (Harlow and Welch 1965). The solver is based on Finite Difference and Finite Volume formulations in the Eulerian and Lagrangian frameworks. The flow of a Newtonian, incompressible fluid with free surface and density (ρ) in a bounded domain is governed by the conservation of mass and momentum equations as follows:where u = velocity vector; t = time; P = pressure; ν = kinematic viscosity; and g = gravitational acceleration acting in the z-direction. FLOW3D-Hydro solves the transient Navier–Stokes equations [Eq. (2)] using Fractional Area/Volume Obstacle Representation (FAVOR) and volume of fraction (Hirt and Nichols 1981) method. FAVOR defines solid boundaries within the Eulerian grid and determines the fraction of areas and volume in partially blocked volume to compute flows corresponding to those boundaries. Thus, boundaries and obstacles are defined independently of grid generation. The Volume of Fluid (VOF) method has been applied in FLOW3D-Hydro to track the free surface motion. To compute the time evolution of water surface, FLOW3D-Hydro applies the following equation for the fraction of fluid function by preserving the step-function nature of the distribution (Hirt and Nichols 1981):where u = velocity vector; t = time; and F = fraction of fluid function.

A two-dimensional flow domain representing the experiment setup described by Sabeti and Heidarzadeh (2021) was defined within the FLOW3D-Hydro [Fig. 2(a)]. The physical experiments were performed using the 4.0-m long wave flume at the Brunel University London and applying a solid concrete block as the sliding mass (Fig. 2). In our numerical model, the entire flow domain is 0.26 m wide, 0.50 m deep, and 1.8 m long. The fluid inside the flume was specified as water with a density of 1,000 kg/m³ at 20°C. The submergence depth (d) and the water depth (h) were 0.080 and 0.375 m, respectively [Fig. 2(a)]. Three solid concrete blocks with the shape of a prism, having variable volumes, are used in our simulations for landslide generation to provide a range for slide volumes (Table 1). These blocks are labeled as small, medium, and large with geometrical parameters listed in Table 1. The specific gravity (γ_s) for these blocks is 2,600 kg/m³ corresponding to the laboratory’s actual blocks. The slope angle in the setup is 45°. A nested grid comprising of two mesh planes with different spatial resolutions has been used to solve the equations (Fig. 3). The finer mesh with the grid size of 0.0010 m is applied in an area of 0.6 m (x-direction) × 0.5 m (z-direction) around the landslide generation. In particular, it covers the impact area of the landslide with the free surface enabling a detailed reproduction of the landslide generation phase (e.g., η_max). The coarser grid, with a size of 0.00175 m, was applied to the other parts of the computational domain (Fig. 3). We examined the model’s sensitivity to grid ratio (size of the larger grid to the smaller one). Spurious reflections at the interface can occur if the ratio of the size of the larger grid to the smaller one is too high (e.g., more than 4). Flow Science (2021) recommends that the grid ratio to be less than or approximately two; it is 1.75 in our simulations. The total number of computational cells is about 9.2 × 10⁶. The grid also covers the initial air space above the water to accommodate the block motion and surface waves. We acknowledge that our model is a simple approximation of real-world landslide events; however, it is adequate for this study as our scope is sensitivity analysis of landslide parameters.

All boundary surfaces of the coarser mesh, except the top surface, were defined as walls; the top surface was designated as a symmetrical boundary. The top, front, and back of the finer mesh area were defined as symmetry, and the other surfaces were of wall type with no-slip conditions around the walls. The symmetry boundary condition indicates that the conditions outside of the boundary are identical to those inside the boundary. As the free movement of waves from the front and back boundaries of the finer mesh are allowed into the coarser mesh, we designated them as symmetry boundaries. In addition, the top boundaries of both fine and coarse meshes are of symmetry types because they represent free water surfaces. Among the available options within the FLOW3D-Hydro for simulating turbulent flows [i.e., Large Eddy Simulation (LES); k–ɛ model and Renormalization Group (RNG)], the RNG model (Yakhot and Orszag 1986) is used in this work. The RNG model is selected because of its wider applicability than the standard k–ɛ model and for its accuracy in modeling turbulent flows (Choi et al. 2007). In the simulations presented herein, the landslide movement has been reproduced by the coupled motion object. Rather than dictating the movement of the slide through the prescribed velocity, the coupled motion object employs the fluid governing equations and gravitational and control forces as well as momentums to model the slide motions (Wei 2005). The friction coefficient in the coupled motion is designated as 0.40–0.50 based on Coulombic friction measurements in the laboratory to calibrate the model (Sabeti and Heidarzadeh 2021). This friction is influenced by the materials of the solid block and the incline. To match the time-dependent intervals of the outputs with the actual wave gauge’s sampling rate, the output interval is set to 0.02 s in the numerical model. Two wave gauges located around the submerged slide are defined in the model to measure the free surface wave fluctuation; these gauges’ locations are fitted in the numerical setup corresponding to the actual locations in the physical experiments (Figs. 2 and 3).

The vertical amplitude of the waves in our experiments is relatively lower than the horizontal wavelength. However, such small amplitudes of the waves can be captured using our mesh with moderate cell densities in the vertical direction using the TruVOF (Volume of Fluid) approach which is unique to FLOW3D-Hydro (Flow Science 2021). TruVOF is a split Lagrangian method that typically produces lower cumulative volume error than the alternative methods. At each time step, the changing areas and volume fractions for each cell are calculated/updated and, the time step size is automatically adjusted to be able to capture the dynamic free surface correctly as it changes within a cell. To guarantee stability and convergence of simulations, FLOW3D-Hydro applies the Courant Number (following equation) condition and ensures that it stays sufficiently below one:where C = Courant Number; Δt = time step; Δx = grid size; and U = flow velocity. FLOW3D-Hydro applies a dynamic time step (Δt) for simulations which means the time step is changing due to various flow conditions. In our simulations, the initial time step was 0.001 s and varied in the range of 0.00035–0.0077 s during the simulations.

Numerical results are compared with previous physical experiments (Sabeti and Heidarzadeh 2021) to examine the model’s capability in the reproduction of real-world measurements and its validation (Oberkampf and Trucano 2002) (Fig. 4). This validation is conducted using the large block (Table 1) and at two wave gauges (WG-1 and WG-2; Figs. 2 and 4). As this study is aimed at investigating the wave characteristics around the source region, it is believed that the two wave gauges around the sliding mass are sufficient. Consistent with the results of the physical model in which the Friction coefficient was in the range of 0.40–0.50, we varied this coefficient in the same range in our numerical model (Fig. 4). Based on the comparison shown in Fig. 4 between modeled wave (i.e., simulations) and the actual one (i.e., observation), our numerical model delivers satisfactory performance. We note that the numerical model was calibrated to produce the best match with laboratory observations by adjusting the friction coefficient. The quality of the match between observations and simulations was measured using the following equation:where ɛ = mismatch error; Obs_i = observation point from physical experiments; and Sim_i = simulation results. Results indicate that the best match is achieved using a friction coefficient of 0.45 as its mismatch errors in the minimum (Table 2); therefore, it has been used for our simulations. The model reproduces the maximum and minimum wave amplitudes as well as the wave period satisfactorily (Fig. 4). There is a small deviation between simulations and observation immediately after the wave trough. Although the model has reproduced such deviation, its magnitude does not necessarily match the observations.

In order to study the sensitivity of the numerical results to the mesh sizes, we varied the mesh sizes by doubling and halving the current meshes and presented the results in Fig. 5. It can be seen that the results from the larger mesh (Fig. 5, Δx1 = 0.0035 and Δx2 = 0.002) are deviating approximately 8% from those of the other two finer meshes. We also see that the waveforms from the two finer meshes are approximately identical (Fig. 5), which demonstrates that the results become stable. Therefore, we are confident that the selected mesh sizes (Fig. 5, mesh sizes 0.00175 and 0.0010 m for coarse and fine meshes, respectively) are suitable for this study. A qualitative comparison of wave snapshots at different times between observations and simulations is presented in Fig. 6, which indicates a good agreement between them. It is useful noting that the model (Fig. 6(a)) reproduces multiple wave troughs and peaks observed in the laboratory (Fig. 6(b)). In general, our numerical model is accurate enough to carry out this research. As we are primarily concerned with the maximum wave amplitude (peak or trough), the numerical model predicts them well. We note that our study is limited to submarine landslides and to solid block failures.

The logarithmic law-of-the-wall is used in FLOW3D-Hydro rather than calculating the viscous sub-layer as a separate region in a turbulent flow. The program calculates the law-of-the-wall velocity in the mesh cells next to an obstacle and applies the momentum equations to give continuity of the velocity profile beyond. The logarithmic layer normally extends in the region 30 ≤ y⁺ ≤ 0.1δ⁺ where y⁺ = (y u_T/v) is the dimensionless wall distance; δ⁺ = (δu_T/v); δ is the thickness of the boundary layer; $u_T=\sqrt{{\tau }_w/\rho }$ is friction velocity; τ_w is wall shear stress; ρ is fluid density; v is the kinematic viscosity; and y is the absolute distance from the wall (Flow Science 2021). For estimating y⁺, we need to manually estimate the shear stress τ_w. In case manual estimation of τ_w is not possible, multiple simulations can be performed to iterate toward a best fit. In general, y⁺ should exceed 30 and remain below several 100’s (e.g., 300–500). Where the wall shear stress (τ_w) is unknown, y⁺ cannot be determined a priori, and iterative meshing of the simulation may be required to find a suitable size. In our simulations, y⁺ ranges from 30 to 163 which guarantees that our meshing process are accurate (Flow Science 2021).

In total, we performed 177 simulations applying our validated model by using three different sliding block sizes (small, medium, and large; Table 1). In our simulations, we considered main landslide parameters that control wave amplitudes namely: submergence depth (d), landslide volume (V), water depth (h), and travel distance (S) (Fig. 1). The other important parameter is the slope angle (θ); however, we used a fixed slope angle of 45° in our research to limit the scope of this work. The simulations are grouped into two main classes: (1) Those with varying submergence depth (d) and constant travel distance (S); and (2) Those with varying travel distance and constant d.

A sequence of three velocity vector plots at different times during landslide motion is shown in Figs. 12 and 13 for the large and small blocks, respectively. The direction of the vector points in the flow direction and the length of the vectors indicate the magnitude of the velocity (Figs. 12 and 13). A complex flow pattern involving outward, inward, and upward flows of water is visible throughout this sequence. A combination of inward and upward flows forms a circulating flow pattern or eddies around the top side of the block (Figs. 12 and 13). A previous study by Fritz et al. (2004) on subaerial landslides showed that such a circulating flow is not generated by subaerial mass movements. The sequence of large block starts at t = 0.060 s after initial landslide motion [Fig. 12(a)] and continues with a time step of 0.2 s covering a time span of approximately 5 s (Fig. 12). The high-water pressure in front of the block pushes the water up towards the surface to form a crest [Fig. 12(c)]. At the beginning of the slide motion, wave amplitude increases in elevation at t = 0.060 s [Fig. 12(a)]. This increase continues until the maximum initial wave amplitude occurs at t = 0.42 s [Fig. 12(c)]. We observe that the largest particle velocity in the wave field occurs below a developed wave trough when the block reaches the end of the slope [v = 0.88 m/s; Fig. 12(c)] as previously reported by Fritz et al. (2009). Compared to the maximum landslide velocity (i.e., the velocity of the solid block) of 0.87 m/s which was measured in our laboratory experiments (Sabeti and Heidarzadeh 2021), the maximum water particle velocity (0.88 m/s) is approximately of the same magnitude. The smaller block (Fig. 13) results in lower maximum water particle velocity (v = 0.688 m/s) and smaller maximum initial landslide amplitude in comparison to the larger block [Figs. 12(c) and 13(c)].

We applied the model FLOW3D-Hydro to solid-block landslide-generated waves in order to investigate the effects of various landslide parameters on the dynamics of the waves. The motivations for this study were to provide a robust and flexible tool for studying landslide tsunamis and to generate additional data on the effects of slide volume (V), travel distance (S), and submergence depth (d) on maximum initial landslide amplitude (η_max). We limited the study to solid-block slides with three sizes of small, medium and large blocks and to a single beach slope (45°). By changing d in the range of 0.04–0.193 m, we observed an inverse relationship between maximum initial landslide amplitude and submergence depth and that the maximum negative initial amplitude varies more significantly than the maximum positive amplitude of the landslide-generated waves. It was also observed that the dominant wave period generated by the large block was 0.7 s whereas it was 1.1 s for the small block. The latter unexpected observation could be attributed to the relatively lower velocity of the sliding mass for the case of the smaller block. For experiments with varying travel distance (S = 0.165 − 0.342 m), an inverse relationship between S and η_max was seen: the shorter the travel distance, the larger the wave amplitude. Velocity vector plots revealed that a circulating flow pattern (eddies) is formed around the top side of the sliding blocks. The FLOW3D-Hydro model validated for the case of a solid-block submarine landslide can be used by engineers and scientists working on tackling landslide tsunami hazards, given the fact that experimental works are generally costly and time-consuming. We add that our model is most valid for the generation phase of landslide-generated waves, and additional studies are required to investigate the propagation and runup of the waves.

All data used in this research are provided in the body of the article.

The authors are sincerely grateful to Mr. Mark Keating (Flow Science Inc.) for his invaluable guidance and support during the modeling process. We gratefully thank Professor Vijaykumar G. Panchang (Texas A & M University, USA) for his constructive comments on this manuscript before submission. The license for using the software FLOW3D-Hydro was granted by Flow Science Inc. to Brunel University London. The authors are sincerely grateful to three anonymous reviewers and the editor for their constructive comments. This study is funded by the Royal Society (the United Kingdom) Grant No. CHL\R1\180173.