Numerical Study of an Impulse Wave Generated by a Sliding Mass

An impulse or tsunami wave is a large wave generated by the displacement of a big volume of water. The triggering cause can be an impulsive geophysical event as earthquakes, subaerial/ submerged landslides, rock falls, and snow avalanches.

On the one side, the knowledge on tsunami generation from sudden impulsive events, as underwater earthquakes, is nowadays considered satisfactory. This was possible thanks to computational advances that led to better modelling of tectonic source mechanisms and coastal inundation models, and, additionally, from large-scale laboratory experiments, analytical results and post-tsunami field surveys [1]. In such case, because of the deformation speed of the marine bed with respect to the propagation of the wave, it is possible to directly relate the deformation of the soil with the wave characteristics.

On the other side, the tsunamigenesis from submarine–subaerial landslides is still considered a not fully understood phenomenon. Indeed, in this case the wave generation process can not be considered impulsive, and the addition of source terms to the equation of motion is needed for the correct modelling [2]. Moreover, the triggering causes may results highly diverse and affected by several factors, as the form and deformation of the sliding mass. Despite many advances have been done in this direction, many researchers still agree that more work is needed to determine how landslide flow processes that involve strong and gradual deformation of subaerial and submarine land are related to tsunami generation [3].

The case of an impulse wave generated by a subaerial landslide is typical of Alpine environments, favored by the presence of steep valleys and large basins at high altitude. A remarkable example is the Vajont dam disaster, occurred in October 1963 in northern Italy: during initial filling, a massive landslide from the surrounding Monte Toc caused a mega-tsunami in the artificial lake. Consequently, 50 million cubic meters of water over-topped the dam and poured into the valley, leading to the complete destruction of several villages and

towns, and to 1,917 deaths.

On the other hand, one of the clearest and disastrous examples of tsunami generated by submarine landslide is the one occurred in Papua New Guinea in 1998. The triggering cause was an earthquake of magnitude VII (Mercalli intensity scale), which, despite not leading to significant direct damage, caused a 2-m vertical drop in the pacific plate along 25 miles of shore. This unleashed a huge tsunami, with 10 m average height and a run-up of 7–8 m, which devastated a 30 miles shoreline and led to a total of around 2200 officialvictims.

The run-up is the height reached by the water above the normal sea level, and it is the most important factor to take into account when evaluating the consequences led by tsunamis. Indeed, except for the largest ones, tsunami waves does not break but they come onshore as a rapid rising turbulent surge of water. Consequently, they appear like an endless tide with a maximum run-up that, historically, has reached up to 30 m a.s.l. [4].

The Papua New Guinea event brought to light for the first time the dangerousness of this type of phenomenon, demonstrating how minor earthquakes can lead to disastrous consequences if they produce submarines slips [5].

The prevention of damages caused by impulse waves and the design of countermeasures requires a detailed knowledge of the characteristics of the waves that can be generated in a particular basin. On the one hand, the experimental study by means of small-scale prototypes allows a good understanding of the phenomenon. However, the detailed simulation of largescale events and the variability of key parameters are more easily taken into account through the numerical simulation. Relevant parameters that may affect the outcomes are the shore slope, landslide mass and shape, basin depth and configuration.

The numerical simulation of this type of phenomena requires the development of a robust platform that takes into account the mutual interaction of liquid, solid and gaseous phases. In this paper we focus on fixed mesh methods, which proved to be the most efficient and adaptable to a wide range of situations. Some examples are those of Liu et al. [2], which uses a large-eddy-simulation (LES) approach and volume of fluid (VoF) algorithms to track the interfaces between phases. Yang et al. [6] introduces the immersed boundary method (IBM) method for the solid tracking while the interface between fluids is captured by a level-set (LS) method. More recently, Sanders et al. [7] and Calderer et al. [8] also solved the interaction between solid and multiphase fluids by means of Fluid Structure Interaction (FSI) methods.

In this work, we propose a numerical platform to simulate the effects of a landslide phenomenon on a water basin. In particular, we perform a direct numerical simulation (DNS) of a two-phase flow on a fixed grid framework, thus, resolving accurately all the length scales involved in the problem. A free surface approximation [9, 10] is adopted in order to simulate the movement of the water independently from the air effect, thus, obtaining a more stable and efficient calculation. The interface tracking is carried-out by means of a conservative LS method [11], while the interaction between solid and liquid is solved by means of direct forcing IBM [12]. An adaptive mesh refinement strategy (AMR) [13] is adopted in order to dynamically improve the mesh definition in the interface region and in zones where the basic mesh size is not sufficient to correctly solve the turbulent scales or the interfacial characteristic lengths. After presenting the numerical framework in Section 2, we propose two simple validation tests in Section 3. Finally, the main test case involving the slip of a sliding object inside a water basin is presented in detail in Section 4.

In this section, a general description of the numerical features involved in the simulation of impulse waves generated by a sliding mass is given. The simulations reported in this work were carried out on the TermoFluids computational fluid dynamics (CFD) software package [14] where all the necessary techniques were implemented and tested.

The movement of the interface is tracked by means of an interface-capturing method, characterized by a scalar level-set function, $\phi(x, t)$, that embeds the moving interface on a fixed grid. $\phi(x, t)$ is a regularized distance function which varies between 0 and 1 in the transition region between the phases. It is advected by means of a transport equation, which in the case of divergence-free velocity field is written as

In detail, this work uses the conservative level-set (CLS) scheme implemented and verified in [11]. The CLS method performs a level-set reinitialization at every time step, according to the following equation

where the artificial diffusion term has the role to avoid discontinuities and to keep the interface thickness constant, while the compressive flux term acts as an interface sharpener.

The level-set function is used to define the position of each phase in the domain, and, consequently, to assign the fluid properties, as density and viscosity. As an example, the density distribution is evaluated as

being r_l and r_g the density of liquid and gaseous phase, respectively. The velocity field, u(x,t), needed for the advection of the level-set function is obtained from the solution of momentum equations. A multiphase flow system is globally governed by continuity and Navier-Stokes (NS) equations, which in a domain occupied by two incompressible fluids separated by an interface $\Gamma$, correspond to

where S includes all the source terms, as surface tension and gravity force. These two equations are discretized over an unstructured grid according to the finite-volume collocated scheme. The discrete solution is obtained by applying a fractional step projection method together with an explicit time integration scheme. In particular, the fluid motion is solved according to a single-phase scheme aimed at optimizing the solution of free-surface problems on 3-D unstructured meshes. The scheme, presented in detail in previous works [9, 10], is based on a novel treatment of the interface for the deactivation of the light phase, allowing the optimization of the classic two-phase model for the cases in which the influence of the lighter phase is negligible. The deactivation is performed by directly imposing the appropriate pressure at the surface boundary and without the need of extending the velocity field to the interface-zone.

The interaction between the two-phase flow and a solid moving with prescribed velocity is solved by means of a second-order direct forcing immersed boundary method (IBM), described in detail by Favre et al. [12]. The IBM method introduces a particular treatment of the Navier-Stokes equations at the interface between solid and fluids. In particular, we consider the following NS equation in discrete form

where RHS includes convective, diffusive and source terms and $\vec{f}$ is an additional source term used to consider the effects of the solid motion on the fluid. $u^{\mathrm{P}}$ is the predicted velocity introduced in the Fractional Step Method for the explicit advancement of the Navier-Stokes equations [15]. The forcing term, $\vec{f}$, is evaluated as

where $U_s$ is directly prescribed in the nodes interiors to the solid body or must follow a particular second-order interpolation [16] for the forcing points that coincide with the interface between solid and fluids.

In case of using Cartesian meshes, the adaptive mesh refinement strategy (AMR) is used to reduce consistently the grid spacing in the vicinity of the interface, thus, allowing an accurate representation of interface phenomena, both at the interface between liquid and gas and between fluids and solid. In particular, the quad/octree algorithm proposed by Antepara et al. [13] and extended to multiphase flows in previous works [17, 18] is here used. The strategy includes a particular treatment of the mass fluxes which allows the exact conservation of mass, momentum and kinetic energy (key element for the correct resolution of turbulent flows [19]). The refinement process is automatically activated by position sensors, which refine the mesh in the regions in which the interface will move. Supplementary refinement is activated in zones with high vorticity, in order to improve the solution of the smallest convective scales in turbulent problems. An example of a second level mesh refinement on a 2D mesh (quad refinement) is shown in Fig. 1, referring to the water entry validation case proposed in the next section: a parent cell refines into 16 child cells in the proximity of the interfaces.

The numerical method is validated by simulating the entry of a 2-D disk of diameter D = 1 m within a initially quiet liquid surface. The disk falls with velocity u = 1 m/s, and the initial distance from its centre to the water surface is set to d = 1.5 m. The domain consists of a rectangle of sides 10 x 10 m, to which Neumann boundary conditions are applied. The liquid density and viscosity are set to $\rho=1000 \mathrm{~kg} / \mathrm{m}^3$ and $u=1 \times 10^{-3}$ Pa s, respectively, and Earth gravity is assumed. The employment of the free surface method implies the negligibility of the gas influence on the flow. Additionally, the effect of the surface tension is considered irrelevant.

The basic mesh has a regular Cartesian grid that accounts for 1 x 104 elements. However, the local definition is increased in the vicinity of the liquid-gas and fluid-solid interfaces thanks to a second level adaptive mesh refinement, as shown in the screenshot of Fig. 1. Consequently, throughout the simulation, the mesh reaches an approximate average of 1 x 105 cells, which dynamically adapt to the position of the free surface.

The profiles of the free surface at different instants of dimensionless time, $t^*=t u / h$,, are presented in Fig. 2. The results are compared with those of Calderer et al. [8], which performed the same test, showing a good agreement.

As a further demonstration, we propose the 3D case of the sphere entry in water. In this case, the falling object is a sphere of diameter D = 0.3 m, and the domain is a rectangular box of sides 2 x 2 x 1 m, half-filled with water. The sphere centre is initially located at d= 0.3 m on the water surface and it falls with velocity u = 0.7 m/s. The Cartesian based mesh, that counts 2.56 x 105 elements, is refined according to the octree criterion and reaches an average of 1.5 x 106 cells. A second level mesh refinement is used; hence, each parent cell refines into up to 64 child cells in the interface region. The screenshots presented in Fig. 3 show the evolution of the flow, until the total immersion of the sphere and the following realistic propagation of circular waves.

The main objective of this work consists in the validation of our model for the case of a landslide sliding within a water basin. The benchmark for this experiment was developed by Liu et al. [2], while additional details about set-up and results are presented in Ref. [20]. In an attempt to evaluate accurately the run-up and run-down generated by three dimensional sliding masses, they performed several experiments, both in laboratory and by means of numerical simulations. In particular, we will consider the case of the subaerial wedge-shaped landslide. Experimental data for free-surface fluctuations and run-up/run-down at various locations are taken as reference to validate our numerical model.

As shown in Fig. 4 (a), the domain consists of a rectangular box of sides 6.6 x 3.7 x 3.7 m, cut by a 1:2 slope. The mass is represented by a wedge of sides 0.91 x 0.455 x 0.61 m, which follows the same slope. The object is initially raised by D = 0.45 m above the water level , hw = 3m . In a general case, the equation that describes the movement of a solid body that slips inside a water basin has been proposed by Pelinovsky [21]. For simplicity, in this work we use a simplified version of the profile reported by Wu [20] to prescribe the motion of the body (subaerial case), which was obtained directly from experimental measurements.

The complete time history profile is reported in Fig. 4(c), showing both experimental and simplified curves. Figure 4(b) depicts the position of the sensors (S4 to S13) used to monitor the oscillations of the free surface, while run-up gauges 2 and 3 are used to measure the run-up and run-down in two different points of the coast. Figure 4(d) shows the computational mesh used for the simulation. It consists of an unstructured tetrahedral mesh, with increased cell density in proximity of the free-surface, in order to achieve a better precision in the measurements of surface fluctuations. It accounts for nearly 1.8 x 106 cells.

Figure 5 shows the results obtained in the measurement of run-up against simulation time. Run-up 2 demonstrates a good agreement with experimental data [2], with a maximum run-up slightly over-predicted at t = 2.7 s ($R_{n u m}=0.13 \mathrm{~m}$ versus $R_{\exp }=0.1$). The same trend is reflected in the run-up 3 plot. However, in this case a short delay of the numerical solution can be noticed after a few seconds of simulation.

The time measurement of relative free-surface fluctuations, h, are shown in Fig. 6. For simplicity only S4 to 8 and S10 probes are reported. In this case, in addition to present work and experimental data, we add to the comparison the numerical results by Yang and Stern [6] which performed the same experiment to validate their two-phase numerical model. Also in this case, the measurements satisfactorily follow those of reference. It is possible to notice that present work particularly fits to Ref. [6] in some time frames in which they both deviate from the experimental solution. This is clear, for example, in free-surface peaks of S5 and S7. This can lead to the conclusion that part of the mechanical behaviour of the apparatus used in the laboratory experiment lead to spurious effects on the liquid unpredictable by numerical simulations.

Finally, in Fig. 7 the comparison between pictures taken from the flow simulated numerically by Ref. [2] and present work are reported for different time instants. Despite the approximated form of the sliding wedge velocity profile adopted in this work, a good overall agreement is found. The model correctly represents the flow from the initial immersion of the object until its stop (t = 2.7 s) and the subsequent propagation of waves on the free surface and the coast.

In this paper we demonstrate the capability to integrate the level-set aided free surface model with the Immersed Boundary scheme to simulate the interaction of a body with prescribed motion and the surrounding fluids. The method can be used on unstructured or tetrahedral meshes dynamically refined by AMR. The platform is altogether suited to the simulation of land sliding into water pools. In fact, the waves induced on the free surface faithfully reproduce those measured through laboratory experiments, allowing the evaluation of the dangerous run-up effect of tsunami waves. Consequently, the tool allows a comprehensive study of this important phenomenon or similar ones, by taking into account several of the many variables that may influence them.

In future work, our intention is to totally couple the movement of solid and fluid by means of fluid-structure interaction algorithms, in order to obtain a more complete and robust platform.