Numerical Modeling of Tsunamis Generated by Subaerial, Partially Submerged, and Submarine Landslides

Abstract

Using the existing two-dimensional experimental data and Open-source Fields Operation and Manipulation (OpenFOAM) software, this study performs a comprehensive comparative analysis of three types of landslide-generated tsunamis (subaerial, partially submerged, and submarine). The primary objective was to assess whether numerical simulations can accurately reproduce the experimental results of each type and to compare the predictive equations of the tsunami amplitudes derived from experimental and simulated data. The mesh size and dynamic viscosity parameters were initially optimized for a specific partially submerged landslide tsunami scenario and then applied across a broader range of experimental scenarios. Most of the simulated wave amplitudes remained within the 50% error margin, although significant discrepancies were observed between landslide types. When focusing on the crest amplitude of the first wave, the simulations of subaerial landslides least deviated from the experimental data, with a mean absolute percentage error of approximately 20%, versus approximately 40% for the partially submerged and submarine landslides. The predictive equations derived from the simulations closely matched those from the experimental data, confirming that OpenFOAM can effectively capture complex landslide–tsunami dynamics. Nonetheless, variations in the coefficients related to slope angles highlight the need for further calibration to enhance the simulation fidelity.

1. Introduction

Tsunamis can be generated not only by seismic-induced ground deformation but also by large-scale landslides. Indeed, landslide-generated tsunamis have caused large changes in water surface elevation and devastated human communities throughout recorded history. Landslide-generated tsunami events are categorized into three primary types based on the relative positions of the landslide concerning the initial water surface: subaerial, partially submerged, and fully submerged (submarine).

Subaerial landslides, which occur above the water surface, are sources of severe tsunami events [1,2,3,4,5,6]. For instance, a tsunami generated by a subaerial landslide at Lituya Bay, Alaska, USA, in 1958 produced the largest recorded tsunami run-up (524 m) worldwide [1,2]. Another massive subaerial landslide tsunami at the Vajont Dam in Italy caused approximately 2000 fatalities in the downstream village [4]. Partially submerged landslides, induced by displacement of a section of land above and below the water surface, have also generated large tsunamis [7,8,9,10]. For instance, during the 2010 Haiti earthquake, approximately 400 m of the coastline collapsed into the sea, generating tsunami waves up to 3 m high [7]. More recently, the 2018 flank collapse of the Anak Krakatau volcano in the Sunda Strait, Indonesia, involved the movement of volcanic material above and below the water surface, producing tsunamis that devastated coastal communities on Java and Sumatra [10]. Submarine landslides, which occur entirely underwater, have also generated catastrophic tsunamis [11,12,13]. Examples are the 1929 Grand Banks event off the coast of Newfoundland, Canada, and the 1998 Papua New Guinea earthquake. The latter event, with tsunami run-ups exceeding 15 m at Sissano Spit, resulted in over 2200 fatalities [12,13]. The 2018 Palu earthquake in Indonesia potentially involved all three types of landslide-generated tsunamis, causing approximately 4000 fatalities in Palu City and Donggala Regency [14,15,16,17,18,19].

To date, the characteristics of landslide tsunamis have been experimentally studied in two-dimensional flumes and three-dimensional wave basins. Subaerial and submarine landslide tsunamis, in particular, have been extensively investigated in physical landslide models using solid blocks [20,21,22,23,24] or deformable materials [25,26,27,28,29,30,31,32,33,34]. Some studies have utilized solid-block and deformable materials to clarify their impacts on generated tsunamis [34,35]. Although fewer in number, physical experiments have been conducted on partially submerged landslide tsunamis [32,34,36,37,38,39]. All such experimental studies have provided valuable insights into the dynamics of landslide tsunamis and have often provided empirical equations for predicting the maximum water surface elevation of a landslide tsunami [21,22,25,26,27,29,30,32,34]. The dynamics of landslide-generated tsunamis have also been explored through analytical approaches [24,40]. For instance, Liapidevskii et al. [40] developed an analytical model for shallow turbidity flows, providing valuable insights into the behavior of the failing mass and the subsequent wave propagation.

Other researchers have developed numerical models of landslide tsunamis. In the early 1990s, water and landslide movements were modeled using nonlinear shallow-water equation models (i.e., nondispersive depth-integrated models) [41,42,43,44,45,46]. For instance, a nonlinear two-layer model, with one layer simulating the behavior of sea water, and the other simulating the behavior of landslides, as proposed by [43], has been utilized to simulate submarine landslide tsunamis in studies by [44,45,46]. However, although nonlinear shallow water equation models have accurately simulated the tsunami behaviors of coseismic tsunamis [47], landslide tsunamis reportedly fall into the intermediate water-wave regime, where dispersion effects are non-negligible [28,48]. Consequently, landslide tsunamis, particularly their propagation stages, have been better modeled with nonlinear and dispersive Boussinesq equations [28,48,49,50,51,52]. For instance, Dutykh and Kalisch [51] applied Boussinesq modeling to examine surface waves resulting by submarine landslides and demonstrated the effectiveness of this approach in capturing the fundamental physics of wave propagation. However, the complex nature of landslide and water movement during the generation stage of landslide tsunamis is most appropriately captured by three-dimensional models based on the Navier–Stokes equations. The efficacy of these models on landslide tsunamis has been confirmed in comparisons with experimental data of physical models [48,53,54,55,56,57,58,59]. For instance, Kim et al. [48] simulated the behavior of a deformable landslide and the resulting tsunami using the three-dimensional TSUNAMI3D model. They confirmed that the model adequately reproduces the results of physical experiments [60]. Urakami and Yoneyama [57] developed a three-dimensional fluid–rigid body coupled model that replicates the solid landslide tsunami experiments in [36]. They also demonstrated good agreement with the experimental results; in particular, their model captured the time histories of water level changes and maximum water levels. Sabeti and Heidarzadeh [58] verified that solid-block submarine landslide simulations using the FLOW3D software captured well the results of their own experiments.

Open-source Fields Operation and Manipulation (OpenFOAM) is a freely available open-source tool for computational fluid dynamics, which has been actively applied to various ocean and coastal wave problems [61,62]. Owing to its accuracy and flexibility, OpenFOAM has been adopted in several recent studies of landslide tsunamis [63,64,65,66,67,68,69,70]. For instance, Romano et al. [63] demonstrated that OpenFOAM using the overset mesh technique can successfully reproduce a solid-body submarine landslide tsunami. Similarly, Huang and Chan [64] showed that OpenFOAM can accurately simulate tsunamis generated by subaerial and submarine landslides. Sabeti et al. [65] demonstrated that both FLOW3D and OpenFOAM can accurately reproduce the experimental results of subaerial rigid-block landslide tsunamis. OpenFOAM has also effectively simulated tsunamis caused by deformable landslides and can accurately model the complex dynamics of such events [66,67,68,69,70]. Additionally, Rauter et al. [69] used OpenFOAM to model both small-scale laboratory experiments and large-scale catastrophic events, confirming its robust applicability to subaerial landslide tsunami simulations.

Although OpenFOAM and numerous other three-dimensional numerical simulation models have proven applicable to landslide tsunamis, most studies have focused on subaerial or submarine landslides; only a few studies [66] have tested numerical models against the results of partially submerged landslide tsunamis, which have also been recognized as considerable hazards. For instance, the 2010 Haiti Earthquake, the 2018 Sunda Strait Tsunami, and the 2018 Palu Tsunami were triggered by partially submerged landslides. Some existing studies have highlighted the risks associated with these tsunamis in other areas, such as La Yesca reservoir, Mexico [37], and Black Lake, Canada [71]. Therefore, it is crucial to ascertain whether the three-dimensional models, generally considered most suitable for simulating landslide tsunamis, are equally effective for partially submerged landslide tsunamis. Additionally, despite the need to compare the accuracies of these models across the three types of landslide tsunamis—subaerial, submarine, and partially submerged—a comprehensive comparative analysis of all three types (to the authors’ knowledge) is currently lacking. Furthermore, experimental studies in this field usually apply multivariate regression to collected datasets, constructing predictive equations for maximum water levels [21,22,25,26,27,29,30,32,34]. Given that numerical simulations are generally intended to incorporate a broad range of new case studies to augment experimental data [34], understanding how regression equations derived from numerical datasets differ from those based solely on experimental results is essential. However, no studies have assessed these differences in the context of landslide tsunamis.

Applying OpenFOAM, the present study attempts to reproduce a wide range of experimental scenarios in [32], encompassing subaerial, submarine, and partially submerged landslide tsunamis. While numerous studies have numerically simulated landslide-generated tsunamis [41,42,43,44,45,46,48,49,50,51,52,53,54,55,56,57,58,59,63,64,65,66,67,68,69,70], to the authors’ knowledge, no previous studies have systematically compared the accuracy of simulations based on different landslide types (i.e., subaerial, submarine, and partially submerged landslides). Additionally, most existing studies have validated simulation accuracy by primarily comparing time histories or focusing on maximum water levels. However, the accuracy of simulating other important parameters, such as wave celerity, periods, and wavelengths, has not been thoroughly investigated. Thus, the primary objective of the present study is to elucidate the accuracy differences among simulations of the three types of landslide tsunamis through a comprehensive comparative analysis. The analysis focuses on a range of key parameters, including the maximum and minimum water levels, wave periods, propagation speeds, and wavelengths. In addition, the authors develop and compare predictive regression equations derived from experimental and numerical datasets of partially submerged tsunamis, focusing on the discrepancies in the coefficients of these equations. This is particularly significant in understanding how regression equations derived from numerical datasets differ from those based solely on experimental results. Through this comprehensive approach, the authors seek to clarify the applicability of OpenFOAM to diverse landslide tsunami scenarios.

2. Methodology

2.1. Numerical Model: OpenFOAM

The numerical analysis was performed in OpenFOAM. Among the various available solvers, the “multiphaseInterFoam” solver was selected because it can simulate incompressible multiphase flow. The governing equations of the solver are the continuity equation (Equation (1)) and the Navier–Stokes equation (Equation (2)), which are discretized using the finite volume method and then solved with the PIMPLE algorithm, a combination of Pressure Implicit with Splitting of Operators and the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE):

$$\nabla ·\mathbf{u}=0,$$

(1)

$${\displaystyle \frac{\partial \rho \mathbf{u}}{\partial t}}+\nabla ·\left(\rho \mathbf{u}\mathbf{u}\right)=-\nabla p^{*}+\nabla ·\left(2\left(\mu +{\mu }_t\right)\mathbf{D}\right)+\rho \mathbf{g},$$

(2)

where $\mathbf{u}$ is the velocity vector, $t$ denotes time, $p^{*}$ is the pseudo-dynamic pressure, $\mathbf{D}$ is the strain rate tensor, $\mathbf{g}$ is the gravitational acceleration, and $\rho$, $\mu$, ${\mu }_t$ are the local fluid density, molecular dynamic viscosity, and eddy viscosity, respectively. The behavior of multiple phases (water, air, and landslide) was modeled using the following volumetric phase fraction equation:

$${\displaystyle \frac{\partial {\alpha }_i}{\partial t}}+\nabla ·\left({\alpha }_i\mathbf{u}\right)+\sum _j\nabla ·\left({\alpha }_i{\alpha }_j{\mathbf{u}}_{r,ij}\right)=0,$$

(3)

where ${\mathbf{u}}_{r,ij}$ is the relative velocity between the phases, and ${\alpha }_i$ is a phase indicator function denoting the proportion of a grid cell occupied by fluid i. Specifically, ${\alpha }_i$ equals 1 when the cell is completely filled with fluid i and 0 when it is completely occupied by other fluids. The local fluid density ($\rho$) and molecular dynamic viscosity ($\mu$) at any point ($\mathbf{x}$) and time ($t$) are obtained by summing the weighted contributions from the densities ${\rho }_i$ and viscosities ${\mu }_i$:

$$\rho =\sum _i{\alpha }_i\left(\mathbf{x},t\right){\rho }_i,$$

(4)

$$\mu =\sum _i{\alpha }_i\left(\mathbf{x},t\right){\mu }_i.$$

(5)

To effectively solve the above equations, it is necessary to appropriately obtain the eddy viscosity (${\mu }_t$) using a turbulence model. Among the range of available models (e.g., the k−ε model) in OpenFOAM, the authors selected the shear stress transport (k−ω–SST) model. This model blends the strengths of both the k−ω model, which performs well near walls, and the k−ε model, which is effective in regions farther away from the wall. The 2003 version of the k−ω–SST model [72], used in the present study, includes several enhancements, such as a production limiter that constrains the turbulent kinetic energy production to prevent over-prediction. Compared to the k−ε model, the 2003 k−ω–SST model demonstrates improved stability and accuracy in predicting flows near the wall, making it particularly well suited for simulating complex fluid flows. This selection was further supported by its successful application in similar landslide tsunami simulations in OpenFOAM [70].

2.2. Benchmark Simulation

The numerical approach was validated on the results of previous two-dimensional experiments [32] as a benchmark. The experiment was conducted in a wave flume of length 14.5 m (not fully utilized), width 0.4 m, and depth 0.8 m (Figure 1). The slope section was constructed from an acrylic board and the remaining horizontal section featured a smooth concrete bottom. During these experiments, deformable landslides composed of glass beads were triggered by the sudden release of a gate positioned above the slope.

The experimental variables were the diameter of the glass beads (d), the slope angle (θ), initial water depth (h), gate height (h_i), initial submergence depth (h_s) for partially submerged and submarine landslides, vertical drop distance (h_a) for subaerial and partially submerged landslides, and landslide mass (m) (each parameter is geometrically defined in Figure 2). The authors of [32] varied these parameters to simulate all three types of landslides, thereby creating 61 subaerial, 103 partially submerged, and 97 submarine landslide tsunami scenarios. Among these scenarios, the present study focuses on those generated by 20 mm diameter glass beads and slope angles ranging from 30° to 60°, which include 27 subaerial, 42 partially submerged, and 27 submarine cases. The experimental results of [32] demonstrated that differences in diameter of glass beads (e.g., 4 mm vs. 20 mm) had negligible effects on the time histories of water levels, justifying the use of only the 20 mm diameter cases in this study. Additionally, slope angles of 8° and 15° in [32] were excluded due to the relatively small water level changes observed in these cases, which resulted in higher measurement uncertainties. Given that the primary objective of this study is to compare the accuracy of simulations for subaerial, submarine, and partially submerged landslide-generated tsunamis, the selected 81 cases with steeper slope angles were considered sufficient to achieve this goal. All cases were simulated in OpenFOAM. The ranges of the experimental parameters imposed by the selection are listed in

Table 1. Ranges of dimensional and nondimensional parameters used in the present study, extracted from the benchmark experiments [32].

Parameters	Subaerial	Partially Submerged	Submarine
θ (°)	30, 45, 60	30, 45, 60	30, 45, 60
m (kg)	5, 10, 15	5, 10, 15	5, 10, 15
h (m)	0.20, 0.25, 0.30, 0.35, 0.40	0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55	0.40, 0.45, 0.50, 0.55, 0.60, 0.65
hi (m)	0.40	0.20, 0.40	0.20, 0.40
ha (m)	0–0.20	0.01–0.23	–
hs (m)	–	0.05–0.25	0.02–0.27
M	0.08–0.94	0.04–0.60	0.03–0.23
A	0–1.00	0.02–0.93	–
S	–	0.11–0.56	0.03–0.67
I	1.00–2.00	0.44–0.89	0.33–0.73
X	2.50–30.00	1.67–24.00	1.54–15.00

. To encapsulate the dynamics across different scenarios, several nondimensional variables—relative landslide mass, M = m/(ρ_wbh²); relative vertical drop distance, A = h_a/h; relative submergence depth, S = h_s/h; relative initial position of landslide, I = h_i/h; and relative streamwise distance, X = x/h, where ρ_w denotes the density of water and b denotes the flume width—were also defined and are summarized in Table 1.

Figure 1 shows the placement of wave gauges (WGs) in [32] and the present study. The accuracy of the simulated results was assessed on the time histories of the water level change obtained from these four WGs.

2.3. Numerical Settings

To replicate the experimental setup, the present study implemented a two-dimensional numerical wave flume with a length of 9.0 m, a width of 0.02 m, and a height of 0.80 m (Figure 3). Cartesian meshes were initially generated using the “blockMesh” utility and the slope area was then refined with triangular elements generated by “snappyHexMesh”. The optimization of the mesh size is detailed in Section 2.4. Impermeable boundary conditions were applied to the slope and bottom of the flume.

The density and dynamic viscosity were set to 1000 kg/m³ and 1.00 × 10⁻⁶ m²/s, respectively, for water and 1.00 kg/m³ and 1.48 × 10⁻⁵ m²/s, respectively, for air. The landslide material (glass beads in the benchmark experiments) was modeled with a density of 2500 kg/m³ to match the experimental conditions. To capture the porosity of the glass beads (35%), the voids in the landslide initially located above the water surface were filled with air, while those in the submerged region were filled with water. For simplification purposes, the landslide was modeled as a Newtonian fluid. Although a more sophisticated landslide rheology is ideal for a detailed representation [69], simplifying the landslide as a Newtonian fluid maintains an acceptable accuracy. However, as highlighted elsewhere [61,69], the dynamic viscosity must be accurately calibrated when modeling landslides under Newtonian assumptions. To ensure that the simulation accurately reproduces the landslide tsunamis, a parameter study on dynamic viscosity was conducted as detailed in Section 2.4.

2.4. Sensitivity Studies

The present simulation was optimized by varying the mesh size and dynamic viscosity parameters and selecting the best values. Three mesh sizes—0.02, 0.01, and 0.005 m—were uniformly applied throughout the simulated area. The dynamic viscosity was varied as 1.00 × 10⁻², 1.00 × 10⁻³, and 1.00 × 10⁻⁴ m²/s. Combining these variations, nine distinct simulation setups were tested. As a case study, the optimization process was conducted on a partially submerged landslide tsunami with the following parameter settings: θ = 45°, m = 10 kg, h = 0.50 m, h_i = 0.40 m, h_a = 0.08 m, and h_s = 0.10 m.

Figure 4 illustrates the time histories of water surface elevations derived from the sensitivity studies. The results are organized into a matrix of 12 panels: three columns corresponding to the mesh sizes (0.02, 0.01, and 0.005 m from left to right, each collecting the results of the three dynamic viscosities) and four rows corresponding to the wave gauges (WG1, WG2, WG3, and WG4 from top to bottom).

Overall, increasing the distance from the generation (slope) zone (i.e., from WG1 to WG4) improved the alignment between the experimental and simulated water surface elevations under all numerical conditions, as the generation zone is more affected by the movement of landslides, which would not be perfectly captured with the assumption of Newtonian fluid. Nevertheless, when the viscosity was set to 1.0 × 10⁻² m²/s, the model consistently underestimated the observed magnitude of the first wave, even at WGs further from the generation zone. In contrast, after lowering the viscosity to 1.0 × 10⁻⁴ m²/s, the water levels (particularly at WG1) tended to exceed the experimentally observed levels. Moreover, the deviations from the experimental results at WG1 became more pronounced in the second and subsequent waves, suggesting that the accuracy of the landslide modeling decreases as the landslide mass moves further down the slope.

When varying the mesh size at a fixed dynamic viscosity, a slightly greater deviation from the experimental results was observed on the 0.002 m mesh than on the 0.001 m and 0.0005 m meshes. Based on these findings, the dynamic viscosity and mesh size were set to 1.0 × 10⁻³ m²/s and 0.001 m, respectively. These settings were considered as the best balance between computational accuracy and efficiency.

Figure 5 illustrates the simulated tsunami generation process under the optimized settings of a partially submerged landslide tsunami. At the beginning of the simulation, the landslide pushed against the water body, initiating a first wave (Figure 5a) that propagated from the slope (Figure 5b). As shown in the subsequent panels, the landslide continually moved downward and a second wave was generated on the slope (Figure 5c). The second wave propagated offshore with a relatively high velocity near the water surface (Figure 5d).

3. Results

3.1. Wave Amplitudes

Figure 6 displays the simulated wave amplitudes (relative to the experimental results) in the 96 cases analyzed after optimizing the simulation parameters in the sensitivity studies. The first and second columns display the relative maximum and minimum amplitudes of the first wave (a_c1/h and a_t1/h, respectively), and the right column displays the relative maximum amplitude of the second wave (a_c2/h). The plots correspond the results obtained at the four WGs [all 96 cases (i.e., a total of 384 plots) are shown in the top row, and those at slope angles of 45°, 30°, and 60° are demonstrated in the second, third, and fourth rows, respectively].

In the partially submerged and submarine landslide tsunami cases, the simulation tended to overestimate the observed maximum amplitude of the first wave (a_c1) (Figure 6a). In the subaerial landslide tsunami cases, some of the simulated amplitudes were overestimated, while others were underestimated. Specifically, in subaerial landslides with a slope angle of 45°, the data points locate close to the identity line, indicating relatively high accuracy, but at slope angles of 30° and 60°, the simulated values are generally below and above the experimental values, respectively. The simulated results of the partially submerged and submarine landslide tsunamis typically exceeded the experimental values; moreover, the discrepancy increased with increasing slope angle. To fairly compare the accuracies of the landslide tsunamis across types, the authors calculated the mean absolute percentage errors (MAPEs) in each case (

Table 2. Comparison of mean absolute percentage errors (MAPEs) in simulations of subaerial, partially submerged, and submarine landslide tsunamis at different wave amplitudes (%).

	Subaerial	Partially Submerged	Submarine
ac1 for all	21.94	43.15	42.47
ac1 for 45°	16.33	29.61	35.84
ac1 for 30°	16.89	25.69	37.02
at1 for 60°	41.01	67.61	69.75
at1 for all	20.64	19.62	8.89
at1 for 45°	15.42	6.97	9.00
at1 for 30°	23.09	11.30	5.14
at1 for 60°	31.24	37.48	14.56
ac2 for all	27.97	21.39	17.48
ac2 for 45°	16.81	12.88	16.80
ac2 for 30°	33.03	16.27	20.15
ac2 for 60°	50.80	33.09	15.10

). The MAPEs of a_c1 in the subaerial, partially submerged, and submarine landslide tsunamis for all angles were 21.9%, 43.1%, and 42.6%, respectively. When focusing on the differences among the slope angles, the MAPEs at slope angles of 30° and 45° are approximately half that of the MAPEs at 60°, regardless of the landslide type.

The MAPE was lower for the minimum wave amplitude of the first wave (a_t1) (Figure 6b) than for a_c1. Specifically, the MAPEs in the a_t1 values of partially submerged and submarine landslide tsunamis (all slope angles) were 19.6% and 8.9%, respectively, significantly lower than those of a_c1. In addition, most data points of these landslide types closely aligned with the identity line at slope angles of 45° and 30°, indicating higher accuracy at these angles than at 60°. Interestingly, simulations of the subaerial landslide tsunamis yielded the smallest MAPE for a_c1 and the largest MAPE for a_t1. Simulations of the subaerial landslide tsunamis substantially underestimated the a_t1 values at 45° and 30° and slightly overestimated them at 60°. It should also be noted that the MAPE for a_t1 was lowest at a slope angle of 45° in simulations of subaerial and partially submerged landslide types, indicating that more accurate simulations for a_t1 were achieved at 45° than at other angles.

The simulated maximum amplitudes of the second wave (a_c2) (Figure 6c) in the submarine landslide cases were relatively well aligned with the experimental results (MAPE for all slope angles = 17.5%). In the subaerial and partially submerged landslide tsunami cases, the MAPE of a_c2 exceeded that of a_t1 for all slope angles. Simulations of these landslide types tended to underestimate and overestimate the experimental results at slope angles of 30° and 60°, respectively. In fact, the MAPE in these landslide types was minimized for a_c2 at 45°, indicating that at this slope angle, the simulations most accurately simulated the second wave in subaerial and partially submerged landslide tsunamis.

Table 3. Comparison of MAPEs at different WGs for subaerial, partially submerged, and submarine landslide tsunamis (%).

	Subaerial	Partially Submerged	Submarine
ac1 for all at WG1	23.04	34.27	27.67
ac1 for all at WG2	19.96	38.88	34.42
ac1 for all at WG3	20.57	46.49	47.30
ac1 for all at WG4	24.18	52.97	60.48
ac1 for 45° at WG1	19.38	21.58	23.26
ac1 for 45° at WG2	13.74	26.40	30.08
ac1 for 45° at WG3	15.04	32.09	42.72
ac1 for 45° at WG4	17.16	38.37	47.28
ac1 for 30° at WG1	13.96	20.66	27.62
ac1 for 30° at WG2	16.43	24.84	28.19
ac1 for 30° at WG3	17.83	25.69	37.02
ac1 for 30° at WG4	19.33	31.56	69.75
ac1 for 60° at WG1	41.27	55.48	40.09
ac1 for 60° at WG2	39.07	60.13	56.58
ac1 for 60° at WG3	37.14	73.90	77.10
ac1 for 60° at WG4	46.57	80.95	105.24

summarizes the MAPEs of a_c1 calculated at each WG. As shown, the MAPEs tend to increase with the distance from the generation zone, indicating that the simulation accuracy decreases as the waves propagate offshore. This reduction in accuracy is more pronounced for partially submerged and submarine landslide tsunamis compared to subaerial landslides. For instance, the MAPE (for all slope angles) at WG1 is 23.0% and at WG4 is 24.2% for subaerial tsunamis, while for partially submerged tsunamis, the MAPE increases from 34.3% at WG1 to 53.0% at WG4, and for submarine tsunamis, from 27.7% to 60.5%. Additionally, the accuracy is shown to decrease as the slope angles become steeper, regardless of the landslide type.

3.2. Wave Periods, Celerity, and Wavelength

Figure 7 compares the first wave characteristics of the simulated and experimental results across the considered scenarios. The three columns of this figure display (from left to right) the tsunami wave period T_c1, the propagation speed (celerity) C_c1, and wavelength L_c1. Here, the T_c1 was obtained by measuring the time difference between the arrivals of the first and second wave crests at each WG. The propagation speed C_c1 of the tsunami wave was calculated as the quotient of the distance between adjacent WGs and the passage time of the first wave crest between these gauges. The L_c1 was calculated as the product of the propagation speed between adjacent WGs and the average of the wave periods measured at these gauges. The top row of Figure 7 displays the aggregated results of all slope angles, while the second, third, and fourth rows are the results at slope angles of 45°, 30°, and 60°, respectively. The MAPEs of these measurements are detailed in

Table 4. Comparison of MAPEs of the simulated wave periods, celerity, and wavelength (%).

	Subaerial	Partially Submerged	Submarine
Tc1 for all	8.92	4.48	4.32
Tc1 for 45°	7.73	3.61	3.86
Tc1 for 30°	13.62	5.58	4.33
Tc1 for 60°	7.21	4.68	5.58
Cc1 for all	5.33	6.69	9.63
Cc1 for 45°	5.01	6.53	8.86
Cc1 for 30°	4.72	5.66	9.25
Cc1 for 60°	6.71	7.49	12.39
Lc1 for all	10.55	8.42	10.89
Lc1 for 45°	9.04	7.17	9.88
Lc1 for 30°	15.04	6.95	10.59
Lc1 for 60°	9.84	10.59	14.20

.

The T_c1 plots in the partially submerged and submarine landslide tsunami cases predominantly aligned near the identity line (Figure 7a), indicating that the wave period was accurately simulated by OpenFOAM. In fact, the MAPE ranged from 3.6% to 4.7% in simulations of partially submerged landslides and from 3.9% to 5.6% in simulations of submarine landslides. In both instances, the wave periods gave lower MAPEs than the wave amplitudes (Table 2). The wave periods were also simulated more accurately than the wave amplitudes in the subaerial landslide cases, but the simulated results tended to be underestimated, particularly at milder slope angles. The MAPEs in the T_c1 for subaerial landslide tsunami cases were 7.7% at 45°, 13.6% at 30°, and 7.2% at 60°. The simulated propagation speeds (C_c1) generally aligned well with the experimental results, with MAPEs of 5.3% for subaerial, 6.7% for partially submerged, and 9.6% for submarine landslides for all slope angles. Interestingly, the accuracy increase (i.e., MAPE decrease) of C_c1 with tsunami type opposed that of T_c1, being highest in subaerial landslides, followed by partially submerged and then submarine landslides. In this study, the wavelength (L_c1) was calculated by multiplying the wave period by the propagation speed. As the MAPEs in the wave period were smallest in the submarine cases and the MAPEs in celerity were smallest in the subaerial cases, the MAPE in the wavelength was minimized (8.4%) in partially submerged landslides. The MAPE is expected to be smallest at the 45° slope angle because the parameters were optimized at this angle. Contrary to this expectation, the MAPEs were minimized at various slope angles and were inconsistent across landslide types.

3.3. Predictive Equations

This subsection establishes two sets of predictive equations for the normalized maximum amplitude of the first wave (a_c1/h): one derived from experimental data alone [32] and the other derived from simulation data alone. Based on the results recorded at the four WGs across all partially submerged landslide cases, these equations incorporate the nondimensional relative landslide mass (M), relative submergence depth (S), relative initial position of the landslide (I), tanθ, and relative streamwise distance (X). The relative vertical drop distance (A) was excluded because it can be estimated from M and S, so it is not a completely independent variable. The predictive equations of a_c1/h were developed through multivariate regression analyses of their distinct data sets. The equation derived from the experimental data [32] is given by

$${\displaystyle \frac{a_{\mathrm{c}1}}{h}}=0.519M^{1.262}{S}^{-0.372}{I}^{0.682}{\left(\tan \theta \right)}^{-0.004}{X}^{-0.423}.$$

(6)

Similarly, the equation derived from the simulation data is

$${\displaystyle \frac{a_{\mathrm{c}1}}{h}}=0.717M^{1.245}{S}^{-0.275}{I}^{0.743}{\left(\tan \theta \right)}^{0.276}{X}^{-0.359}.$$

(7)

Both equations describe the influence of the specified nondimensional parameters on the maximum wave amplitude, highlighting the inherent differences between the two datasets. Figure 8a compares the experimentally measured a_c1/h values with those predicted by Equation (6), while Figure 8b compares the simulated a_c1/h against those predicted by Equation (7). The correlation coefficients (R²) of both equations are notably high (0.96 and 0.99 for the experimental and simulation data, respectively), confirming a strong relationship between a_c1/h and the nondimensional parameters of partially submerged landslides.

Comparing the coefficients in each term of the predictive equations, one observes a notable discrepancy only in the tanθ term, with values of −0.004 in Equation (6) and 0.276 in Equation (7). This finding suggests that the simulation cannot accurately replicate the effects of the slope angle but sufficiently captures the M, S, I, and X values. This notable difference in the coefficient of tanθ explains the higher estimates of the predictive equation derived from the simulation results than of the equation derived from the experimental data (Figure 8c).

4. Discussion

The present study began with a sensitivity analysis to optimize the mesh size and dynamic viscosity in simulations of partially submerged landslide tsunamis. For the specific case of the optimization, the parameters were set to θ = 45°, m = 10 kg, h = 0.50 m, h_i = 0.40 m, h_a = 0.08 m, and h_s = 0.10 m and the mesh size was varied as 0.01, 0.02, and 0.005 m. As the analysis results negligibly differed on the 0.01 m and 0.005 m meshes, a mesh size of 0.01 m was selected to reduce the computational load. This same mesh size was selected in a previous analysis of landslide tsunamis [68], but some studies adopted smaller [67,70] or larger [65] mesh sizes. Three values of dynamic viscosity (1.00 × 10⁻², 1.00 × 10⁻³, and 1.00 × 10⁻⁴ m²/s) were also considered in the sensitivity analysis. The viscosity, which optimized the results (1.00 × 10⁻³ m²/s) and was adopted in subsequent simulations, contrasts with the viscosities used in similar studies assuming Newtonian fluids for landslides. Specifically, an analysis of subaerial landslide tsunamis [59,67], which was dominated by interactions between the landslide material and air, adopted higher viscosities, whereas a study of submarine landslide tsunamis by the same authors [59,67], in which the interactions occurred entirely underwater, adopted lower viscosities. As partially submerged landslides fall between subaerial and submarine landslides, the chosen viscosity of the present study—positioned between the optimal values of subaerial and submarine landslides identified in [59,67]—is deemed appropriate.

Optimizing the parameters under each experimental condition can be time-consuming and is sometimes infeasible, particularly when the experimental data are scarce or difficult to obtain. Therefore, it is necessary to determine the extent to which parameters optimized under a specific set of conditions can be applied to other scenarios. In the present study, the parameters optimized in the partially submerged landslide case were intentionally applied across all scenarios to evaluate their general applicability. Although the optimization was designed to minimize the effort of computing wave amplitudes, most of the simulated values remained within the 50% error margin (Figure 6). Therefore, even with a streamlined optimization approach, OpenFOAM simulations can estimate various landslide scenarios, provided that an error margin of up to 50% is acceptable. This finding demonstrates that a single set of optimized parameters can be applied across different conditions, particularly when extensive parameter optimization for each scenario is infeasible. Essentially, modeling the generation zone of landslide tsunamis is more complex than modeling the propagation stage, as it requires capturing the intricate behavior of landslides. Although a 50% error margin must be considered, the results suggest that OpenFOAM can be effectively applied to simulate the generation process of landslide tsunamis. While simulating the entire landslide tsunami event with OpenFOAM may pose challenges due to computational costs, the generated data could still serve as boundary conditions for other models that simulate tsunami propagation and inundation in real-world scenarios, as demonstrated in studies such as [8,28]. This strengthens the case for using OpenFOAM in predicting and mitigating the risks associated with landslide-generated tsunamis.

When focusing on the results of subaerial cases, the MAPEs of a_c1, a_t1, and a_c2 in the subaerial landslide tsunami cases were 21.9%, 20.6%, and 28.0%, respectively. As the accuracy of the simulations in these cases largely depends on slope angle, the model must be calibrated at the targeted slope angle before applying it to other conditions such as initial water depth and landslide mass. Although the simulations were optimized for partially submerged landslide tsunamis, a_c1 was simulated with relatively low accuracy in these cases (MAPE = 45%), highlighting the challenges of accurately modeling the crest amplitude of the first wave of partially submerged landslide tsunamis, but the a_t1 and a_c2 were simulated with considerably higher accuracies (MAPEs ~20%). Additionally, the accuracy of a_c1 was shown to decrease with wave propagation, especially for partially submerged and submarine landslide tsunamis. As noted in the benchmark experiments by Takabatake et al. [32], the generated waves in this study fall within the intermediate water depth regime. This means that wave dispersion plays a significant role, causing longer-period waves in the generated spectrum to travel faster than shorter-period waves as they propagate offshore. The reduced accuracy at more distant wave gauges suggests that the simulations may not fully capture the energy of the longer-period wave components, especially when the landslide is initially submerged. This discrepancy could be attributed to the assumption of Newtonian fluid dynamics for the landslide. In the partially submerged cases, a_t1 and a_c2 were more accurately reproduced at a slope angle of 45° than at 30° and 60°, emphasizing the need for calibration at the target slope angle for these cases as well. The submarine landslide tsunamis yielded a higher MAPE for a_c1 and lower MAPEs for a_t1 and a_c2, but the simulation accuracies of these cases were negligibly influenced by slope angle. The limited reproducibility of a_c1 in the partially submerged and submarine landslide tsunamis might also be attributable to the optimized mesh size. As the average a_c1 values at the WGs of the partially submerged and submarine landslide tsunamis were relatively small (0.017 and 0.057 m, respectively, versus 0.027 m for subaerial landslide tsunamis), they might be inadequately captured with the adopted mesh size (0.001 m).

Other wave characteristics, namely, the wave period, celerity, and wavelength, were more accurately simulated than the wave amplitudes. As the generated waves fall within the intermediate water depth regime, the effects of wave dispersion cannot be ignored. In this study, the focus for the wave periods was on the first wave, which is primarily influenced by the longer-period components of the generated wave spectrum. The lower MAPEs for wave periods compared to those for wave amplitudes suggest that while the numerical model effectively capture the timing and periodicity of these longer-period components, it struggles to accurately reproduce the wave energy (i.e., wave amplitudes) associated with these components. The lower variability in the MAPEs for propagation speeds can be attributed to the relatively smaller wave amplitudes of the generated waves. When wave amplitudes are smaller, the nonlinear effects that influence propagation speed are minimized, allowing the speed to be predominantly governed by the wave period and water depth. Since the wave periods of the first waves are well reproduced in the simulations, it is likely that the propagation speeds also show good agreement across different landslide types, despite the higher MAPEs observed for wave amplitudes. Notably, these variables were not clearly influenced by slope angle, indicating that they are not predominantly influenced by landslide movement. Thus, the modeling settings, including the assumption of Newtonian fluid dynamics, would not need to be overly sensitive when the purpose of the simulations is to obtain these specific properties.

Numerical simulations of landslide experiments can create new datasets for deriving predictive equations of wave amplitudes [33]. Accordingly, the authors compared the predictive equations derived from the experimental and simulation data. The coefficients of the terms in the two predictive equations were generally consistent, confirming that the simulations accurately reproduce the relationship between the generated wave amplitudes and nondimensional parameters. However, the slope-angle coefficients significantly differed between the two equations. Therefore, when developing predictive equations through numerical simulations, the results must first be optimized with respect to the specific slope angle.

A notable limitation of this study is the treatment of landslides as Newtonian fluids. This assumption simplifies the modeling process by requiring only the adjustment of viscosity to represent the landslide dynamics, thereby reducing the uncertainties associated with parameter settings. The decision to adopt a Newtonian fluid assumption in this study was primarily driven by the need to balance computational efficiency with reasonable accuracy in capturing the key characteristics of landslide-generated tsunamis. Indeed, similarly to previous studies [59,67], the time histories of water surface elevation were relatively simulated by fine-tuning the viscosity values. However, while the Newtonian fluid assumption is relatively effective for simulating the initial landslide movement and the generation of the first wave, it may oversimplify the complex rheological behavior of real-world landslides. In particular, the second and subsequent waves near the generation zone were less accurately simulated under this assumption. This suggests that although the first wave is reasonably well-reproduced near the generation zone, the assumption becomes less reliable as the landslide progresses further down the slope, impacting the accuracy of the subsequent waves. Additionally, while the accuracy of the crest amplitude of the first wave remained relatively high close to the generation zone, it decreased as the wave propagated offshore, particularly for partially submerged and submarine landslide tsunamis. This may be because the longer wave components, which become more prominent as the waves propagate offshore, were not fully captured by the simulations, especially when the landslide begins to move underwater. To address these limitations, a more advanced rheological model is required to accurately simulate the landslide’s movement and its influence on the generated waves. For instance, Romano et al. [70] demonstrated that a Coulomb viscoplastic non-Newtonian rheological model can precisely reproduce the time histories of landslide movements. Other aspects of the numerical modeling, including the boundary conditions and turbulence models, could also be optimized to enhance the simulation accuracy. Refining the mesh size near the water surface might further improve the results, especially when the generated tsunami heights are small. While the present study employed the fixed mesh, alternative strategies such as Adaptive Mesh Refinement (AMR) [73] and moving mesh techniques [74] could offer improvements in both accuracy and computational efficiency. These methods dynamically adjust mesh resolution, enabling more accurate modeling of critical areas such as the water surface. Future research could explore the application of these techniques to enhance the simulation of landslide-generated tsunamis.

It is also important to note that the benchmark experiments [32] were conducted using glass beads with a constant diameter of 20 mm. In contrast, real landslides typically consist of materials with varying particle sizes, which can introduce more complexity into the landslide’s behavior. This variability would require more advanced modeling techniques to better simulate real-world landslide dynamics. To further validate the applicability of OpenFOAM simulations, future studies should explore a broader range of landslide materials with more realistic particle distributions. Finally, it must be acknowledged that the present findings are derived from two-dimensional simulations, which may not completely capture the complexities of three-dimensional, real-world scenarios. In future research, this gap should be bridged by three-dimensional simulations that provide more comprehensive insights into landslide dynamics.

5. Conclusions

Applying OpenFOAM software, the authors conducted a comprehensive comparative analysis of three types of landslide-generated tsunamis: subaerial, partially submerged, and submarine. The primary objectives were to explore whether numerical simulations can accurately reproduce the experimental results of each type, as documented in previous two-dimensional experiments, and to compare the predictive equations of tsunami amplitudes derived from experimental and simulated data. Initially, the mesh size and dynamic viscosity of the landslide were optimized for a specific partially submerged landslide tsunami case. These optimized settings were subsequently applied in OpenFOAM simulations over a broad range of experimental scenarios: 27 subaerial, 42 partially submerged, and 27 submarine landslide tsunamis.

Most of the simulated wave amplitudes fell within the 50% error margin of the experimental results, indicating that OpenFOAM with a single set of optimized parameters can reliably estimate various landslide scenarios when a 50% error margin is acceptable. However, significant discrepancies in the simulation accuracies were noted among the three types of landslides. The deviations between the simulated and experimentally measured crest amplitudes of the first wave were higher in the partially submerged and submarine landslide tsunami cases (MAPEs ~40%) than in the subaerial tsunami cases (MAPE ~20%). In contrast, the trough amplitudes of the first wave and the crest amplitudes of the second wave were more accurately simulated in the partially submerged and submarine cases than in the subaerial case. The slope angle notably affected simulations of subaerial and partially submerged landslide tsunamis, emphasizing the importance of calibrating the model at the targeted slope angle to improve the simulation accuracy. Nonetheless, the wave periods, celerity, and wavelength were simulated with relatively high accuracy, with MAPEs typically less than 10%.

The predictive equations developed from the numerical results closely matched those derived from experimental data, indicating that OpenFOAM can effectively replicate the complex dynamics of landslide tsunamis. However, the coefficients related to slope angles differed between the two equations, suggest that further calibration focused on slope angle would enhance the fidelity of the simulations. These findings are promising but the present simulations based on two-dimensional experiments cannot fully capture the complexities of three-dimensional, real-world scenarios. In addition, the present study modeled the landslide as a Newtonian fluid for simplicity and efficiency. A sophisticated rheology that more accurately models the behavior of landslides will be considered in future research.