VOF/Navier-Stokes numerical modeling of surface waves generated by subaerial landslides

La generation de vagues par des glissements de terrain est etudiee a l'aide d'un modele base sur les equations de Navier-Stokes et un suivi d'interface VOF. L'originalite de l'etude repose sur la prise en compte implicite du couplage glissement/fluide. Le modele est valide dans un premier temps, dans le cas de generation d'un soliton par la chute d'un bloc rigide a l'extremite d'un canal. Une etude numerique de l'influence du caractere deformable du glissement est ensuite proposee. Ce travail met en evidence l'importance et la complexite du role de cette deformation sur les caracteristiques des vagues generees. Il convient donc de prendre en compte de maniere plus fine la rheologie du glissement dans le processus de prediction des tsunamis.


I n INtRODUCtION
Whether subaerial or underwater, landslides may generate surface waves, large enough to cause significant hazards for coastal populations and infrastructures. In the ocean, underwater landslides represent the second most important source of tsunami generation that, sometimes, may create more devastating tsunamis than co-seismic tsunami sources of moderate strength [1]. Aerial and subaerial landslides frequently occur in mountain lakes and fjords, whose sides may become unstable. Although less frequent in the ocean, large subaerial landslides may be induced on volcanic islands susceptible to flank collapse, and cause potentially damaging tsunamis for nearby, or possibly far distant, coastal populations. Indeed, in such cases, slide volumes involved may become so large as to generate major, even catastrophic, transoceanic tsunamis that could potentially cause destruction along far distant coasts. A recently studied case in this respect is the potential flank collapse of the Cumbre Vieja volcano on the island of La Palma, in the Canaries [2].
Laboratory experiments have been performed to study submarine and subaerial landslides (e.g., [3][4][5][6][7]), particularly in relation to case studies (see, [3] for a more detailed discussion of work to date), but these may be both lengthy and costly to carry out. By contrast, numerical modeling, if properly validated with laboratory experiments, may be a more flexible and efficient tool (e.g., [1], [6] and [9]). Moreover, numerical modeling more easily provides flow variables for any point of space and, hence, is better suited to a detailed study of physical processes.
Earlier modeling work of so-called tsunami landslides may be classified into two groups. In the first one, slide kinematics is a priori specified in the model, typically based on a semi-empirical equation, e.g., describing the center of mass motion for solid landslides (e.g., [9][10][11][12][13][14]). This method has more often been applied to underwater landslides, for which induced free surface waves are usually initially less complex, although some authors have also applied it to subaerial landslides ( [6], [15]). In the second group of methods, a fully coupled computation of both the slide and induced fluid motion is performed. Monaghan and Kos [16] thus use a Smooth Particle Hydrodynamics (SPH) method to study shallow water wave generation, due to the impact of a falling rigid block. Panizzo and Dalrymple [17] also used a SPH model to study underwater landslide generated waves. Models based on a direct solution of Navier-Stokes (NS) equations, and featuring a free surface tracking algorithm, have also been used. Heinrich [18]

Risques Littoraux Majeurs
et al. [19] thus proposed two-dimensional (2D) NS simulations, with a Volume Of Fluids (VOF) type free surface tracking, of both underwater and subaerial landslide tsunamis. Gisler et al. [20] recently proposed a modeling of the La Palma case study based on a compressible NS model, and Quecedo et al. [21] similarly modeled the Lituya Bay case study with a multi-fluid NS model, in which both air and water motion were simulated. These various NS model results are all quite promising, but the slide/fluid motion coupling still needs to be clarified, as it has not yet been fully studied. Also, a detailed analysis of physical processes involved has rarely been carried out.
In this work, we present results of simulations of surface waves generated by subaerial landslides, obtained with the NS model Aquilon (developed in the Trefle laboratory, UMR 8508). In section 2, we briefly summarize the equations and numerical methods for the model. In section 3, we present results of a validation case for the 2D generation of a solitary wave by a falling block in shallow water. Finally, in section 4, the full potential of our model is illustrated by simulating the fully coupled case of waves generated by a deformable slide moving down a slope.

II n NUMeRICaL MODeL
The computational domain is divided into three fluid subdomains : water, air and slide. In our approach, the latter is also treated as a fluid, whose properties and constitutive law can be adjusted depending on water content and nature. Pierson and Costa [22] for instance indicate that, for water content greater than 50 %, slides behave as Newtonian fluids. For smaller water contents, however, slides behave as viscoplastic fluids, with various constitutive laws depending on whether behaves more as a granular or a debris flow. The three fluids in the subdomains are similarly governed by incompressible NS equations : where U denotes the velocity vector and P the total pressure. The various fluids that co-exist in the computational domain are represented in Eq. (2) by their space-varying density ρ and viscosity μ. Equations (1) and (2) are discretized on a fixed structured grid, and solved using an augmented Lagrangian formulation [23]. The linear algebraic system obtained for each time step is solved using the iterative BiCGSTAB algorithm [24]. Since turbulence is not simulated in our model, results in fluid areas with strong shear stresses are to be taken with caution. Turbulence effects will be investigated in future work, which should be straightforward considering that different turbulence models are already implemented in Aquilon.
The motion of interfaces is represented by advection equations, expressed for a "volume fraction function" quantifying the fraction of each fluid present in interface cells, speci-fying that these interfaces are moving with the flow. Here, two such equations can be expressed for the two interfaces between the three fluids : where C water (resp. C slide ) is equal to 1 in water (resp. in the slide) and 0 everywhere else. Equations (3) and (4) are either solved using an explicit Total Variation Decreasing Lax Wendroff scheme of LeVeque [25] like in section IV, or with the PLIC-VOF reconstruction algorithm [26] like in section III. Lubin et al. [23] found that both methods lead to the same accuracy in the case of wave breaking simulation.
Note, for interface cells, Eq. (2) is expressed using an equivalent density and viscosity, obtained by weighed average of individual fluid properties, based on their volume fraction in the cell. The whole model, composed of equations (1) to (4), is called a "one-fluid" model [27] because it approximates a three phase flow with a single fluid formulation, in which the physical parameters are space varying. With such models, no conditions are required between the various fluid interfaces. The velocity continuity is implicitly specified in the model formulation itself.

III n VaLIDatION CaSe
The numerical model is first validated by simulating the 2D generation of a quasi-soliton in shallow water of constant depth, by a falling rigid block. This application, known as "Russel's Scott wave generator", has been the object of wellcontrolled laboratory experiments as well as other numerical simulations. Figure 1 shows the basic set-up and parameters for both the experiment and SPH simulations performed by Monaghan and Kos [16]. The experiment took place in a 2D flume, 9 m long and of variable depth D. A 38.2 kg rectangular block (0.4 m tall, 0.3 m long and 0.39 m wide; nearly the same as the inside of the flume) is located just above still water level at initial time, and then released. Experiments were repeated for D = 0.288, 0.210, and 0.116 m; in each case, the block vertical position and free surface elevation at a wave gauge located 1.2 m from the leftward extremity of the flume, were measured as a function of time.
Numerical simulations were performed using a 250 x 130 grid, with uneven cell dimensions in x (starting with Δx min = 6 mm near the block and widening for larger x) and constant cell dimensions in y (Δy = 6 mm). To closely reproduce experiments, in which the block was forced to have a vertical motion, horizontal block velocities are set to zero in the model. The block is represented by a very viscous fluid with μ = 10 4 Pa.s. Caltagirone and Vincent [28] have shown that, with such high viscosities, the local fluid flow within the block behaves as a solid that only experiences a global translation or rotation (i.e. ). Slip boundary conditions are specified for the velocity components along all boundaries. Figure 2 shows numerical results at a few selected times, for D = 0.21 m. As in the experimental set-up, the block was slightly shifted rightward (by 25 mm), so that a very thin water sheet can rise to the left of it during its fall, as can be seen on figure 2. Results show, as expected, that the numerical method allows the block to keep its shape during motion, and hence to stay rigid. The main feature of the water flow is the development of a vortex at the block lower right corner. When time increases, this vortex detaches from the block and advects rightward, while losing strength. Monaghan and Kos [16] also report the occurrence of a vortex at this location, followed in experiments by the development of a small plunging breaker on the free surface, near the block right side. Our simulations, however, cannot resolve such a small flow feature with the selected grid size. Adaptive grid refi-nement should be performed to this effect. Finally, also note on the figure the appearance of two counter-rotating trailing vortices in the air flow, that gradually detach from the upper corners of the block. It should be emphasized that, considering no turbulence model is used, such predicted features should be taken with some caution. Figure 3 shows results for the block fall velocity as a function of depth. Despite the fairly coarse grid used in the simulations, the block velocity agrees well with measurements down to two-third of the water depth. Further down, the simulated block moves slower than in reality. The flow, however, is still well simulated in the model, as evidenced by the good agreement of simulated and experimental elevations of waves at the gauge (Table 1). We conclude that, overall, the fluid/solid kinematic coupling is well represented

Risques Littoraux Majeurs
in our multi-fluid model, and free surface elevations are well simulated by the free surface tracking algorithm.

IV n DeFORMING SLIDe CaSe
We now simulate the generation of waves by a deforming subaerial landslide, moving down a 50 % (26.6 deg.) slope. Similar to the idealized shape used in [9], the slide geometry is assumed to be semi-elliptical, with length L = 1 m and thickness T = 0.2 m. The slide is initially at rest and partly emerged, with its center of mass located at d/L = -0.048 below still water level.
The computational grid is specified in a polar coordinate system, with 200 meshes in the angular direction, between θ = 63.44° and 90°, and 400 meshes in the radial direction, between r = 0.22 and 16 m. The radial step is kept constant to Δr/L = 2.7 10 -2 , between r = 0.22 and 8 m, and then increased exponentially. Slip boundary conditions are specified for the velocity components along all boundaries. The time step is set to 10 -2 s for the first 10 iterations and then adaptively calculated later on, based on a Courant condition with CFL = 0.5.
Simulations were performed for different values of the slide dynamic viscosity : μ s = 5000, 500, 100, 10 and 1 Pa.s, in order to estimate the influence of slide deformation (governed by viscosity) on the generated tsunami characteristics (note, water viscosity is μ = 0.001 Pa.s). Figure 4 shows three snapshots of simulations for cases with decreasing viscosity, μ s = 5000, 100 and 1 Pa.s, which demonstrate the strong effect of viscosity on slide deformation during motion. Velocity vectors, plotted in the water, show how the fluid is being re-circulated from the front to the back of the slide, with the creation of a small trailing vortex. One also sees, in each case, a marked depression of the free surface to the back of the slide, as well as narrow wedges of water, moving up the slope, that will cause tsunami runup. At this particular time, in each case, a first long wave is moving offshore (to the right of each figure), followed by a second somewhat steeper wave (to the left of each slide). Looking more closely, we see that the slide with higher viscosity does behave more as a quasi-rigid solid, closely keeping its front-back symmetrical shape during motion. For the slide with medium viscosity, a bulbous front develops during motion that causes a thickening of the slide, yielding the slower moving slide among these three cases. In the case with lowest viscosity, the slide significantly lengthens, while a large vortex of slide material (here behaving as a heavy viscous fluid) forms at the front, and starts rolling onto itself. Figure 5 shows free surface elevations, generated at almost the same time t = 3.2 to 3.5 s in each case with different slide viscosity. At these times, which correspond to about 1.5 s later than in Fig. 4, the second generated wave has evolved into steeper waves, even breaking in one case (case 3). For cases 1, 4, and 5, the waves are quite similar in height and length, despite the marked difference in slide motion and deformation. For case 3, which corresponds to the slower and thicker slide shown in Fig. 4, the free surface wave steepens up to overturning and plunging breaking; this is the case where the maximum of energy has been transferred from the slide to the water motion, causing the free surface to steepen and break at t = 2.8 s. Note, Case 2 is also found to break at t = 3.5 s, but the volume of water involved in the plunging jet is smaller.
One important point that has not yet been investigated is how much energy is eventually conveyed from the slide to the water motion, once a quasi-stationary stage has been reached; this will be left out for further study. We however calculated the height of the generated tsunami at a larger time t = 4.5 s, and found : H = 23, 19, 8, 21, and 20 cm, for cases 1-5, respectively. Thus, wave height is similar for all cases (with a slightly larger height for the rigid case 1), except for case 3, for which dissipation due to the violent breaking has clearly reduced wave height.

V n CONCLUSIONS
We presented simulations of surface wave generation by subaerial slides of idealized geometry, using a multi-fluid NS model, with a VOF free surface tracking algorithm. The comparison of simulations with experimental results for the generation of a solitary wave by a falling block [16] confirms  The influence of slide deformation on wave generation is studied by performing simulations for various slide viscosity, the largest one corresponding to a nearly rigid slide.
We show that the deforming slide shape strongly influences both slide motion and wave generation. In cases with large slide deformation (thickening) during motion, surface wave breaking may even occur at some late stage of motion. Such results indicate that slide rheology should be included in numerical simulations of landslide tsunamis, in order to improve the prediction of wave characteristics.
An interesting preliminary observation in terms of tsunami hazard assessment is, in our simulations, a rigid slide represents the worst case scenario, yielding the largest offshoremoving tsunami height H. This had also been suggested by Grilli and Watts [9] based on potential flow simulations and using arbitrary slide deformation. Since idealized slide shape and incline geometries were used in the simulations, however, this conclusion may not be general. The present model being general, other slide geometries could be investigated in future work.

VI n aCKNOWLeDGeMeNtS
This work was supported in part by a grant of the "Tsunami Risk And Strategies For the European Region" program, funded by the European Commission, under contract n o 037058.