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

: Submarine landslide-generated waves have been responsible for signi ﬁ cant damage to coastal communities worldwide. Despite this, the existing knowledge on the mechanism of the phenomenon is limited that can be partly attributed to the absence of adequate validated numerical models. In this study, we applied and validated a Computational Fluid Dynamics numerical model (FLOW3D-Hydro) and used it for investigating the sensitivity of landslide-generated waves to the variations of different parameters. This study is limited to solid-block submarine landslides moving downward at a ﬁ xed slope angle of 45°. We conducted 177 simulations applying the validated model by using three different sliding block sizes (small, medium, and large). The experiments revealed an inverse exponential relationship between maximum initial landslide amplitude and both initial submergence depth and travel distance. We 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; this unexpected result could be attributed to the relatively lower velocity of the sliding mass for the case of the smaller block. DOI: 10.1061/(ASCE)WW.1943-5460.0000694 . This work is made available under the terms of the Creative Commons Attribution 4.0 International license, https://creativecommons.org/licenses/by/4.0/.


Introduction
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 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 deepwater 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 solidblock 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.

Data, Methods, and Model Validation
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 3 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 3 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 6 . 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 Fig. 1. Landslide tsunami parameters used in this study. d, submergence depth; h, water depth; S, travel distance; B, landslide length; T, landslide thickness; θ, slope angle; and η max , absolute value of the maximum initial landslide amplitude. Other parameters are t, travel time; and t*, end time of landslide movement. The travel distance is measured from the toe of the concrete block to the end of the slope. SWL = still water level. 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  (Table 1) is used for this experiment.
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  Table 1 for block dimensions. In these simulations, the travel distance (S) is kept constant. 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 = τ w /ρ √ 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).

Results and Discussions
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.

Simulations with Varying Submergence Depth (d ) and Constant Travel Distance (S)
A total number of 156 simulations were conducted to examine the effect of submergence depth on η max ; each of the three blocks was tested for 52 submergence depths (Fig. 7). The submergence depth (d ) was varied in the range of 0.04-0.193 m while travel distance was kept constant. In this case, the water depth (h) was varied to produce varying submergence depth while the travel distance was constant. It is evident that the water depth also affects the results, but its effects are predicted to be smaller than those of the submergence depth as previously shown by Watts et al. (2005). As expected, an inverse correlation is established between (a) (b) Fig. 8. Spectral analyses of the landslide-generated waves from our simulations for the case of (a) large; and (b) small sliding blocks. d represents submergence depth. See Table 1 for block dimensions. Fig. 9. Correlation between dimensionless maximum wave amplitude (η max /B) and dimensionless submergence depth (d/B) for simulations with varying submergence depth (d) and constant travel distance. The black asterisks are experimental data from Watts (1997). See Table 1 for block dimensions.
maximum initial wave amplitude (η max ) and submergence depth (Fig. 9). The numerous simulations presented in Fig. 7 reveal that the maximum negative initial amplitude varies more significantly than the maximum positive amplitude. As an example, for the case of the large block [ Fig. 7(a)], the maximum negative initial amplitude is in the range of 0.005-0.019 m while it varies from 0.004 to 0.009 m for the maximum positive initial amplitude. The same pattern is seen for the other two blocks [Figs. 7(b and c)]. Comparison of the waveforms from the three sliding blocks also demonstrates that η max increases by an increase in the volume of the blocks, although this observation was reported previously (Watts 1998;Murty 2003). For instance, η max of the large black is 2.4 times larger than that of the small block. An important observation from Fig. 7 is that the period of the landslide-generated wave (i.e., the time duration of a full cycle of the wave) increases as the landslide volume decreases that appears to be unexpected. Spectral analyses of the waveforms reveal the dominant wave periods of 0.7 and 1.1 s for the large and small blocks, respectively (Fig. 8). This could be attributed to the relatively lower velocity of the sliding mass for the case of smaller blocks. In general, the larger the size of the landslide, the longer the period of the wave that is generated, given all other parameters are kept the same. However, this is not holding in our study of landslide-generated waves; most likely due to the influence of landslide velocity.
To investigate the effects of submergence depth (d ) on η max , we plotted dimensionless submergence depth (d/B) versus dimensionless maximum initial amplitude (η max /B) which is shown in Fig. 9. Here, B is slide length. This is inspired by previous researches by  Table 1 for block dimensions. In these simulations, the submergence depth, d, is kept constant. Watts et al. (2003) and Grilli and Watts (2005), who showed that η max is strongly influenced by d. It is noted that here we are interested in understanding the relationship between the two parameters rather than developing an equation. We observe a nonlinear inverse relationship between the two parameters as follows based on a power regression analysis and using the nonlinear least squares technique (MathWorks 2020): Relationship (6) is valid for the parameter range of θ = 45°and d/B ∈ [0.20, 1.82]. In Fig. 9, we also plotted η max /B versus d/B reported from the experimental work of Watts (1997) (WTS-97) that shows our numerical results are generally consistent with the experimental works although they are slightly deviated. The quality of regression analysis, commonly reported as R 2 , for Relationship (6) is 0.96. It is noted that Relationship (6) may not be used as a predictive equation for η max ; rather it is meant to show the relationship between η max and d.

Simulations with Varying Travel Distance (S) and Constant Submergence Depth (d )
To investigate the effect of travel distance (S) on maximum initial wave amplitude (η max ), the simulated waveforms for varying S and constant d are plotted in Fig. 10. For each block, seven different travel distances were examined while the submergence depth for all simulations was kept constant at d = 0.08 m. Overall, 21 simulations were conducted. Results show that η max is strongly influenced by travel distance with an inverse relationship between the two parameters: the shorter the travel distance, the larger the wave amplitude. For instance, η max of the medium block was increased from 0.005 to 0.008 m by reducing travel distance from 0.308 to 0.223 m [ Fig. 10(b)]. The initial wave produces a free surface with a relatively large negative phase followed by a smaller elevation phase. Because each block had a different length (B) (Table 1), the range of travel distance for each case was different. As expected, the larger solid block volumes produce larger maximum initial wave amplitude (Wiegel 1955). Fig. 10 reveals that the range of η max is varied from 0.004 to 0.006 m for the small block and it is within domain 0.011-0.016 m for the large block.
Regression analysis using nonlinear least square method (MathWorks 2020) was applied to the simulated results in order to find a relationship between the travel distance (S) and η max (Fig. 11). Parameters S and η max were made dimensionless by landslide length (B) and submergence depth (d ). We added a few data points from the physical experiments of Najafi-Jilani and Ataie-Ashtiani (2008) (JIL-08) to our simulation-based graphs (Fig. 11); however, these experimental data are added only for comparison and are not used in our regression analyses. A similar trend can be seen between our simulations and JIL-08 experiments although there are some deviations. The following relationships are established between η max /B and S/B as well as  Fig. 11. (a) Correlation between dimensionless maximum initial wave amplitude (η max /B) and dimensionless travel distance (S/B) for simulations with varying travel distance (S) and constant submergence depth; and (b) same as (a) but for the correlation between η max /d and S/d. The black asterisks are experimental data from Najafi-Jilani and Ataie-Ashtiani (2008). between η max /d and S/d: Relationship (7)

Analysis of Velocity Field around the Source Region
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  Table 1. Velocity magnitude scale is given in this figure.
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)].

Conclusions
We applied the model FLOW3D-Hydro to solid-block landslidegenerated 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.

Data Availability Statement
All data used in this research are provided in the body of the article.  Fig. 13. Velocity field at different times around the source region as the solid block is moving down the slope for the small block at different times of (a) 0.020 s; (b) 0.370 s; and (c) 0.720 s. Dimensions of the block are given in Table 1. Velocity magnitude scale is given in this figure.