‪Impulse wave generation: a comparison of landslides of block and granular masses by coupled Lagrangian tracking using VOF over a set mesh

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 influence of landslide mobility by conducting numerical simulations on porous landslides. Their results showed that the maximum near-field 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.

The physical situation being considered is as follows.

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).

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.

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⁻⁷, 10⁻⁶, and 10⁻⁸, 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).

To validate the results of the current study, a case with landslide material density of 2,500 kg/m³ 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³. 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.

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.

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.

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.

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.

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 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³, 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.

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.

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³ 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³ exhibits the smallest wave profile peak of approximately 0.295 m. The wave profile peaks increase with density (see Figure 15).

Five water depths were investigated with a constant block density of 2,500 kg/m³. 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.

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.

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³ and Scott Russell waves for densities ranging from 1,800 to 3,300 kg/m³ were carried out. It is found that the formation of water waves corresponding to a landslide density of 2,500 kg/m³ 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.

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