## Abstract

This paper presents the numerical results of impulsive waves generated by landslides of solid block, granular materials and heavy block sinking. An impulse product parameter P is developed and a wide range of effective parameters are studied. The volume-of-fluid (VOF) and overset mesh methods have been used to study landslide-generated tsunamis. Also, a Lagrangian tracking approach coupled with the VOF to simulate the granular movement was developed. The effect of the water reservoir depth, the landslide height, the landslide density and the geometrical parameters on the wave height (elevation) has been investigated using the open-source OpenFOAM software. The results have been presented for dimensionless distances and the normalized geometry of the landslide in the ranges 5–7, and 1–2, respectively. These numbers have been normalized to the height of the landslide (a). According to the results of simulations, the tsunami formation process is divided into three stages, which were analyzed in detail by considering the interactions between the solid and the water reservoir. The Scott Russell wave has the highest impulse product parameter among the impulse wave mechanisms which is 58.6% of the total impulse production. In addition, the duration of the wave propagation has been computed based on the wave height.

## HIGHLIGHTS

The impulse wave generation is modeled.

The importance of design parameters in the reservoir dams are highlighted.

The tsunami formation process is divided into three stages.

The novel numerical method is used to simulation.

Coupled Lagrangian tracking using VOF over a set mesh has been applied.

## INTRODUCTION

Tsunamis (impulsive waves) can be generated by sudden movements of volumes of water induced by landslides, earthquakes, volcanic eruptions, and asteroids impacts. Among these, landslides play a key role in impulsive wave generation. This phenomenon happened in the Lituya Bay in 1958 (Fritz *et al.* 2009), in the Vajont Valley in 1963 (Panizzo *et al.* 2005) and in Papua New Guinea in 1998 (Synolakis *et al.* 2002). Also, this phenomenon can occur in bays, reservoirs, lakes, and islands. In a dam reservoir, the behavior of the impulse wave and the density and height of landslides significantly influence the life of dams and their efficiency (Fritz *et al.* 2003a, 2003b, 2004).

Romano *et al.* (2019) numerically studied tsunamis generated by solid and impermeable landslides and proposed a novel approach to impulse wave generation. Water waves generated by landslides have been widely studied theoretically, numerically, and experimentally (Pelinovsky & Poplavsky 1996; Watts 1998; Liu *et al.* 2005; Panizzo *et al.* 2005; Tinti *et al.* 2005; Daneshfaraz & Kaya 2008; Synolakis *et al.* 2002). Heidarzadeh *et al.* (2020) investigated tsunamis generated by landslides that were caused by volcanic eruptions. They found that the wavelengths of landslide-generated waves are shorter than earthquake-generated waves and they exhibit greater dispersive effects. Mulligan *et al.* (2020) proposed a new method to simulate landslide wave generation. Their technique was an improvement upon the conventional methods that had been used. A comprehensive review of waves generated by landslides was presented by Bullard *et al.* (2019a, 2019b). The causes of wave generation associated with landslides were analyzed by Mulligan & Take (2017). It was found that the main criterion for the wave propagation process is the momentum transfer near the collision zone.

One of the first numerical works in this field was performed by Skvortsov & Bornhold (2007). More recently, Jing *et al.* (2020) developed numerical models to investigate the dispersive effects of water waves generated by an accelerated landslide. Heller (2009) used a combination of numerical, analytical, and numerical methods to study impulse waves. They analyzed landslide strength, movement of the landslide, and the wave run-up (overtopping). It was shown that, with increasing landslide strength, the wave height and run-up elevation vary. Dutykh *et al.* (2012) studied the flow patterns of wave propagation under the impact of a landslide.

In the recent past, scientists have proposed various water treatments (Hasegawa *et al.* 2019; Tom 2021), and explained conservation of water resources (Ali & Ahmad 2020; Sahoo *et al.* 2021), and remediating strategies for environmental (Shahid *et al.* 2019; Ajiboye *et al.* 2020). Nanotechnology has also been proposed considering efficient application (Maiga *et al.* 2020).

There are a few researchers who have investigated the interaction of a landslide with the free water of a dam reservoir using numerical methods. Lindstrøm (2016) examined the inﬂuence of landslide mobility by conducting numerical simulations on porous landslides. Their results showed that the maximum near-ﬁeld amplitude decreased with increasing the landslide porosity. The main objective of that work was to study the characteristics of impulsive waves generated by landslides and heavy blocks.

Shen *et al.* (2021) proposed a risk dynamic mechanism to overcome the flaws of safety control in reservoir dam. It was found that the proposed model can increase the abilities of reservoir dams to operate under emergency situations. Nobarinia *et al.* (2021) proposed a novel model to predict the peak outflow from a breached embankment dam. However, to date, studies have yet to be done to more fully investigate the mechanisms of impulse wave generation in reservoir dams. In this paper, numerical simulation has been performed to assess the effects of impulse wave generation types on the impulse product parameters P.

The rest of the paper is arranged as follows. Section 2 presents details of the mathematical modeling, the dimensionless numbers, and assumptions. Section 3 discusses the details of the geometry, the description of case studies, and details of the numerical technique. After validation of the results, Section 3–5 will present a comprehensive discussion about the results. Finally, engineering advice and design criteria are presented in Section 6.

## GOVERNING EQUATIONS

### RNAS model

*et al.*2005; Ding

*et al.*2007; McDonough 2009; Bedram & Moosavi 2012; Ghaderi

*et al.*2020a, 2020b):where denotes the gradient operator, represents the divergence operation, and is the tensor product operator. The symbols

*p*,

*F*, and

_{σ}*τ*are the pressure, the surface tension, and the deviatoric stress tensor, respectively. The Stokes’ hypothesis and the symmetry of the stress tensor apply. Equation (3) is provided by the volume-of-fluid (VOF) model which governs the volume fraction

*α*of one fluid in the mixture. The VOF method determines

*α*for each cell within the computational domain. Cells having

*α*equal to 1 contain only water and those equal to 0 contain only the gas phase (Kunkelmann & Stephan 2009).

### DPM model

*et al.*1999):

*d*

_{p}is the diameter of the solid particles and it has constant value of 5 cm. The drag force equation is adopted as follows (Li

*et al.*1999; Nazari

*et al.*2021):where C

_{D}denotes the drag coefficient, u

_{p}represents the velocity of particles, and d

_{p}is the particle diameter. In addition, because of the size of the solid particles and the interactions between them, the Brownian force (F

_{Brownian}) is considered. Details of Brownian forces can be found in Nazari

*et al.*(2021). The movement of solid particles on the landslide and their impact with the water surface is coupled by the combination of DPM and VOF, which is a Lagrangian tracking approach.

The physical situation being considered is as follows.

### Landslide model

An inclined embankment of a water body contains a mass of material (shown in Figure 1). The landslide material is triangular in shape with width b and height a. Initially, the mass is above the water surface. It moves under the influence of gravity downwards (and to the right in Figure 1). In the figure, the mass is shown after it has submerged somewhat within the water. The mass slides along the embankment until it reaches an impact location where it stops (depth H). The symbol s represents the streamwise direction of flow. The symbol h represents the initial height of the mass above the impact location. The moving mass initiates a wave (represented by the dashed line drawn a top the free surface location).

*a*, the wide of the landslide

*b*, the water reservoir depth

*H*, and the landslide height

*h*(Ribeiro

*.**et al.*2020; see Figure 1):where

*H** denotes the dimensionless depth,

*h** and

*G*are the dimensionless landslide height, the geometrical ratio and streamwise, respectively.

Wave generation by landslides and turbulence effects were predicted using the Reynolds-averaged form of Equations (1)–(3). This is often called the RANS approach. The *κ* − *ε* turbulence model is used to account for turbulence in the flow. The initial values of the turbulent kinetic energy, *k* and the turbulent dissipation rate, *ε* were 0.735 and 3.835 , respectively. For a complete description of landslide collisions with water, it is necessary to consider physical properties such as the viscosity, the density, and the surface tension coefficient. These parameters are listed in Table 1.

Material . | (kg/m^{3})
. | (N/m) . | (Pa.s × 10^{−3})
. |
---|---|---|---|

Air | 1.21 | – | – |

Water | 977 | 0.7 | 0.0727 |

Landslide | 2,100 | – | – |

Solid particle | 2,100 | – | – |

Material . | (kg/m^{3})
. | (N/m) . | (Pa.s × 10^{−3})
. |
---|---|---|---|

Air | 1.21 | – | – |

Water | 977 | 0.7 | 0.0727 |

Landslide | 2,100 | – | – |

Solid particle | 2,100 | – | – |

## PROBLEM DESCRIPTION AND NUMERICAL TECHNIQUE

The open-source field operation and manipulation (OpenFOAM) CFD software package version 1912 was used to perform the numerical simulations. The OpenFOAM code is written in C ++ and uses the finite-volume discretization method to solve the conservation equations of mass and momentum, along with the equations of state. The base code of the solver is overInterDyFoam (Jasak 2009). This solver can take into account different mesh movement models. The second-order upwind scheme is used to handle the convective terms except the phase-fraction term which is discretized using the Vanleer second-order scheme. The Gauss-linear second-order approach is employed to deal with the diffusion terms. The PISO algorithm is applied to couple the pressure and the velocity components. The under-relaxation factors for the pressure, momentum, and energy equations are 0.3, 0.7, and 1, respectively. In addition, the minimum residuals for pressure, velocity, and phase fraction convergence are 10^{−7}, 10^{−6}, and 10^{−8}, respectively.

Figure 1 presents a schematic view of the domain while Figure 2 shows the computational elements of the overset domain, the background domain, and the porous domain. In order to generate a high-efficiency computational mesh, a block-mesh was utilized. The computational mesh is refined gradually from the outer boundary towards the inner domain by halving the sizes of cells in a sequence of local refinement boxes. It is found that these sizes are large enough to guarantee that the results are not affected by the dimension of the domain.

The overset mesh method is based on the use of two (or more) domains. The outer one (i.e. background domain) allows the motion of one or more inner domains (i.e. floating domains) that contain a solid body. The mutual exchange of information between the two domains is achieved by interpolation. Therefore, the two domains, which overlap each other, can be used to simulate different features of the hydrodynamics problem at hand. In Figure 2 a sketch depicting the features of the method is shown. Contrary to other techniques (e.g. immersed boundary or moving mesh methods), this method offers the great advantage that the resolution around the moving body is extremely accurate (i.e. bodyfitted approach) and remains constant throughout the simulation. Thus, the strength of the overset mesh method lies in its ability to represent complex geometries while maintaining a good quality mesh, especially for large amplitude body motions.

Due to the three-dimensional nature of this study, the computational domain is enclosed by three boundaries. No-slip and the fixed-flux pressure conditions are imposed at the landslide walls, and non-reflective-boundary conditions at the background walls and the far-field boundaries. The landslide is triangular with the dimensions a = 0.2 m and b = 0.5 m. It is initially at rest, at the distance 2 m from the reservoir bottom (Figure 3).

## VALIDATION OF RESULTS

To validate the results of the current study, a case with landslide material density of 2,500 kg/m^{3} is used. Comparison of our data with the numerical and experimental visualizations presented by Liu *et al.* (2005) for a three-dimensional numerical simulation of the wave propagation caused by a landslide is illustrated in Figure 4. The vertical axis *η* refers to the water surface elevation. Figure 4 shows a comparison of the elevation of the water surface between the present work and those of Liu *et al.* (2005). The maximum relative deviation between the water surface elevations for the two cases in Figure 4 is ∼ 2%.

The numerical mesh consists of rectangular cells. These cells are refined gradually from the outer boundary toward the inner domain. A grid-independence test has been carried out to compute the proper number of numerical cells for a convergent simulation. To obtain grid independent results, simulations have been performed using three different mesh topologies and with a landslide material density of 2,500 kg/m^{3}. The three meshes are denoted by A, B, and C. The total number of elements for the respective meshes are: A = 700,000, B = 800,000, and C = 900,000. It was observed that the B and C meshes produced almost identical results along the water surface elevation with an approximate percent error of less than 0.3%. Hence, mesh B was chosen as the preferred mesh to balance accuracy and computational time. A summary of the grid-independence test results is shown in Figure 5. The maximum skewness of the grid is 0.3611, which is suitable for obtaining accurate results.

The numerical predictions were also analyzed with respect to different time-step sizes. Three different values of the time-steps are considered and the water surface elevation was computed, yielding the results presented in Figure 6. The figure shows that results are time-step independent for values of 0.00005 seconds.

## RESULTS AND DISCUSSION

The process of producing a landslide-induced impulse wave involves three stages. These stages are termed as: the formation, propagation, and run-up. This paper focuses on the formation of a single wave. The landslide will move with gravity and impact with the water surface. Each case involves simulations of several positions of landslide geometry. The dimensionless ratios that were studied (h/a) are: 5, 6 and 7. The normalized geometry (a/b) of the landslide takes values of 1, 1.5 and 2.

### Effect of the dimensionless distance

Figure 7 presents the time-elapsed images of the water elevation at the formation stage. Results presented in Figure 7 indicate that larger waves are near the side of the channel. This means that the generated wave starts to move along the free surface. The movement of this wave obliterates the initial shape of the free surface and leads to the propagation of smaller waves in the +x direction. Figure 7 reveals that with the increase in the landslide height (strength), the wave height increased.

### Effects of the normalized geometry of the landslide

Figure 8 shows the formation of the single wave as a function of the normalized dimension of the landslide mass. For smaller values of a/b, there is the creation of greater wave heights. The average elevation of water with normalized landslide mass geometry of 2, 1.5 and 1 and a non-dimensional landslide height of 6 are 1.15, 1.26, and 1.3, respectively. This figure shows that the height of the landslide *a*, has a greater influence on the wave elevation than the wide of the landslide *b*.

In the hills overlooking the dam reservoirs, landslides with different geometry may cause impulse waves. Based on the type of these hills, the landslides with different widths can be separated. The effect of landslide geometry has also been demonstrated by Sun *et al.* (2020). It was found that, with the increase in the normalized geometry, the impulse wave height (wave energy) decreases.

### Effect of landslide mass density

Figure 9 illustrates variation in the wave elevation with respect to time. Results show that the water elevation increases with an increase in the landslide density. The results are to be expected because more dense materials will result in a greater mass (and thus momentum) of the material. The landslide suddenly impacts and transfers its momentum to the water which is initially at rest. Subsequently, surface waves form during the sinking process.

The effect of landslide mass density has also been investigated by Schuster & Highland (2003) and Djukem *et al.* (2020). One of the main features of momentum transferring is the disintegration of the quiescent water surface into smaller droplets at the early stages due to the gravity.

### Effects of the water depth

Water depth is an important parameter for landslide-induced impulse waves. Figure 10 shows the effect of the dam water depth on the wave formation. Results show that the wave elevation increases with a decrease in the dam water depth. The reason for this is related to the effects of the bed. The landslide contact with the bed normally causes an increase in the wave height. It is important to notice that a decrease in the water depth obviously leads to more significant bed effects.

### Granular material

Granular material constituting the landslide mass was also studied. Here, the material is not cohesive and spreads as it slides along the embankment, toward the water. To begin, comparisons between the present work and prior research will be set out. In Figure 11, the results of the present study were compared with the results of Shan & Zhao (2014). The figure shows excellent agreement of water surface elevations at a particular instant in time. The normalized landslide mass geometry (a/b) is 1, and the starting location (h/a) and the density of the mass are 6 and 2,500 kg/m^{3}, respectively. Earthen mounds around dam reservoirs are an example of granular material which may slide and created an impulse wave. The average size of the solid particles is 5 cm.

Figure 12 shows a sequence of contours that reveal the demarcation between air and water regions. Also shown in the images are the corresponding positioning of the granular material. It can be seen that the material does not maintain a triangular shape. Rather, the granular material spreads while it slides. Disturbances in the water can be seen relatively early in the sequence of images with a small disturbance forming far from the sidewall and a larger disturbance at the impact site where the granular material enters the water.

## SCOTT RUSSELL WAVE

Next, attention is turned toward waves generated by a heavy rectangular cross-sectional object sinking vertically into water, again using the VOF and the overset mesh methods. This configuration is known as a solitary Scott Russell wave formation in a long rectangular tank. A combination of VOF and overset mesh was used to clarify details of the solitary wave formation and its propagation. First, we simulated a long rectangular tank with a heavy/dense rectangular block adjacent to the left side, as shown in Figure 13. The generated wave profile of the present work is compared with experimental observations of Ataie-Ashtiani & Shobeyri (2008) (upper case of Figure 13) and Monaghan & Kos (2000) (lower case of Figure 13). Based on the excellent agreement between the present calculations and two prior references, confidence is provided for the current calculations. With this verification completed, various new scenarios will be simulated and discussed. The density of sinking box is 2.5 times the water density and the initial height above the water.

### Effects of the heavy block density

The first issue to consider is the influence of block density on the wave formation. A series of calculations with density ranging from 1,800 to 3,300 kg/m^{3} was performed and the results are shown in Figure 14. It is seen that a lower density block reduces the magnitude of the wave profile. It is well known that single wave profiles in long tanks significantly depend on block density, whereas generally the wave elevation increases with density. The characteristics of the wave is similar to the landslide-induced wave discussed earlier. In Figure 14, wave profiles with different block densities are depicted for the Scott Russell wave case study. The density with 1,800 kg/m^{3} exhibits the smallest wave profile peak of approximately 0.295 m. The wave profile peaks increase with density (see Figure 15).

### Effect of the water depth

Five water depths were investigated with a constant block density of 2,500 kg/m^{3}. The resulting wave profiles are illustrated in Figure 16. During the sinking, the peak of the Scott Russell wave has a maximum velocity of 10.5 m/s, which decreases gradually. Figures 16 and 17 show the effect of the water depth on the wave formation at time 0.7 s. Results show that the water elevation increases with decreasing water depth.

### Effects of the block height

*Δh*and

*g*are the block height and gravitational acceleration, respectively. It was found that blocks with a small height create a minimal solitary wave due to their lower impact momentum. Results of Figures 18 and 19 indicate that with an increase in the block height, the wave profile increases because of the increased mass (and thus momentum) of the block upon impact. An increasing box height causes an increase of the momentum. The wave height increases with increasing momentum. Figure 19 discloses water ascension rates corresponding to the four block heights. Results are displayed at a time of 0.7 s. The crest of the wave increases from 0.29 m to 0.34 m as the box height increases from 0.4 m to 0.55 m. Also, the crest profile is stretched along the vertical direction. Behaviors of the crest profile significantly affect the vortex under the wave.

## THE IMPULSE PRODUCT PARAMETER

The impulse product parameter P applies to landslides and to granular and block impulse waves. Numerical results based on P indicate that the largest waves occur for a block sinking into the water (Figure 20). These large waves are a negative phenomenon, particularly for dam reservoirs.

## CONCLUSION AND RECOMMENDATION

The main focus of the present paper is on the generation of landslide-induced impulse waves and Scott Russell waves. Three different positions have been selected based on the simulations conducted for a VOF and overset mesh using the OpenFOAM C ++ library. Three-dimensional simulations of the impulse wave formation for landslide densities ranging from 2,100 to 2,900 kg/m^{3} and Scott Russell waves for densities ranging from 1,800 to 3,300 kg/m^{3} were carried out. It is found that the formation of water waves corresponding to a landslide density of 2,500 kg/m^{3} is consistent with recent experimental results. The present study opens exciting possibilities for future research relating to the design of dam reservoir. Some conclusions about the different features of the wave generation are summarized as follows.

1. The coupling the overset mesh technique overcomes a drawback of the overset mesh method as far as the modelling of a solid body moving in contact with an impermeable surface. The proposed numerical method can be used for Scott Russell wave generation.

2. The landslide-induced impulse wave process is characterized by three stages. In this study, we focused on the first stage; however, the three stages are:

- Formation of the single wave,

- Propagation of the wave,

- Run-up (overtopping).

3. With the increase in the normalized geometry, the impulse wave height decreases.

4. With the increase in the box height, the impulse wave height increases.

5. With the increase in the water depth, the impulse wave height decreases.

6. With the increase in the heavy block density, the impulse wave height increases.

7. The impulse product parameter of the Scott Russell wave is calculated to be 4.98. This value for the landslide and granular are 3 and 1.89, respectively.

As the landslide density and the dimensionless distance of the landslide increase, the magnitude of the wave increases. Also, with the decrease in water depth, the wave elevation increases.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## REFERENCES

^{®}: a new approach