A two-layer shallow flow model with two axes of integration, well-balanced discretization and application to submarine avalanches

Abstract

We propose a two-layer model with two different axes of integration and a well-balanced finite volume method. The purpose is to study submarine avalanches and generated tsunamis by a depth-averaged model with different averaged directions for the fluid and the granular layers. Two-layer shallow depth-averaged models usually consider either Cartesian or local coordinates for both layers. However, the motion characteristics of the granular layer and the water wave are different: the granular flow velocity is mainly oriented downslope while water motion related to tsunami wave propagation is mostly horizontal. As a result, the shallow approximation and depth-averaging have to be imposed (i) in the direction normal to the topography for the granular flow and (ii) in the vertical direction for the water layer. To deal with this problem, we define a reference plane related to topography variations and use the associated local coordinates to derive the granular layer equations whereas Cartesian coordinates are used for the fluid layer. Depth-averaging is done orthogonally to that reference plane for the granular layer equations and in the vertical direction for the fluid layer equations. Then, a finite volume method is defined based on an extension of the hydrostatic reconstruction. The proposed method is exactly well-balanced for two kinds of stationary solutions: the classical one, when both water and granular masses are at rest; the second one, when only the granular mass is at rest. Several tests are presented to get insight into the sensitivity of the granular flow, deposit and generated water waves to the choice of the coordinate systems. Our results show that even for moderate slopes (up to 30°), strong relative errors on the avalanche dynamics and deposit (up to 60%) and on the generated water waves (up to 120%) are made when using Cartesian coordinates for both layers instead of an appropriate local coordinate system as proposed here.

Introduction

In this work we deal with the dynamics of sub-aerial and submarine landslides or avalanches entering the sea or a lake. As it is well known, these avalanches generate a wave at the water free surface and, depending on their intensity, may cause a tsunami. The study and modeling of these natural events are thus a key issue to assess the associated hazards. While dry avalanches have been extensively studied, submarine granular flow models and generated tsunami are still extremely challenging due to the more complex processes at work such as the fluid/grain interaction involving pore pressure effects or friction at the interface between the avalanche and the fluid layers. Derivation of such models for real avalanches over complex topography is even more tricky. This requires to make approximations for granular and fluid flows to get reasonable computational cost.

Many models have been developed to describe avalanches and generated waves as summarized in [39]. The difficulty lies in the different dynamics of the granular and fluid (i.e. water) masses. On one hand, to correctly capture the generated wave and its latter propagation, it is essential to consider dispersive effects in the fluid model as it has been widely observed (e.g. [38], [36], [2], [13], [25]). On the other hand, for the granular mass the difficulty mainly lies in defining the appropriate rheology including pore pressure effects [16], [15], [7]. At last, the interaction between the two materials plays an important role. As pointed out in [16] and precedent authors' works, the tsunami generation depends on the initial acceleration, volume and submergence of the landslide while the granular flow dynamics itself depends on the water flow on top of it.

Essentially two ways of simulating landslides generated tsunamis have been proposed.

The first family of models uses two independent models for the landslide and the tsunami wave. First, the granular flow is simulated using a landslide model. Then, the landslide free surface is imposed as a boundary condition in the water model to compute the wave generation and propagation. Some of these models are presented in [2], [17], [38], [31] where Navier-Stokes or Boussinesq equations are used for the water model. The main advantage of these models is to accurately describe the water wave including non-hydrostatic and thus dispersive effects. Other attempts consider a Boussinesq or a non-hydrostatic shallow water model to include dispersive effects, following the same idea and results (see e.g. [3], [13], [25], [37]).

Nevertheless, using a prescribed avalanche shape as a boundary condition in the water model prevents taking into account the interaction between both materials. It is generally argued that this interaction can be neglected due to the small size of the water waves compared to the large depth of the sea [24], [13]. This is however not true for landslides occurring near the shore where the water motion can affect the avalanche dynamics. As a result, even if these models well simulate the generated tsunami wave, they are not able to correctly capture the landslide dynamics and deposit.

The second family of models consider the interaction between landslide and water in a coupled model. Although there are several models considering Navier-Stokes equations for both layers as in [1], [19], [29], the most usual approach is to use shallow depth-averaged models for the landslide evolution (e.g. [22], [14], [18], [30], [26]). These models are mainly based on the pioneer work of Savage and Hutter (SH) [35] that is a shallow model for dry avalanches composed of granular material whose rheology is based on the Mohr-Coulomb theory. One of the advantages of shallow depth-averaged models for the landslide part is that their computational cost is much lower than Navier-Stokes models while they still describe the interaction between water and granular mass when two layers are considered. In such models, the modeling of the water layer is made using either depth-averaging with or without dispersive effects or with the Navier-Stokes equations as in [30], [26].

The accurate simulation of the granular mass motion at the generation stage is a key point to correctly approximate the initial wave at the free surface, and then, the tsunami evolution. In most depth-integrated models, the same coordinate (Cartesian or local) system is used for both layers, even though it is appropriate only for one [15]. As discussed above, shallow depth-averaged models for water motion should be derived in Cartesian coordinates while the avalanche flow should be described based on a coordinate system related to topography variations. The avalanche generates a wave at the water free surface just above it (i.e. at the same longitudinal position). When using topography-related coordinate system for both the avalanche and water layers, the model generates a perturbation at the surface at a wrong place (see sketch in Fig. 1). The steeper the slope is, the wronger the approximated starting location of the generated wave will be, thus inducing a wrong tsunami propagation. Some works have introduced different coordinate systems for each material, as [18] and [30]. In [18] a bilayer depth-averaged model for submarine avalanches is presented with Cartesian coordinates for the water layer and local coordinates for the avalanche. The drawback is that the landslide shape is considered as boundary condition for the water layer equations so that no interaction between the fluid and granular layers are taken into account. As discussed above, this simplification is justified in their paper by the fact that the action of the water wave on the avalanche is not so important for deep avalanches and because the water waves amplitude is small compared to the water depth [24]. In [31], the authors present a 3D shock-capturing method for the Navier-Stokes equations, named NHWAVE. This solver is based on σ-coordinates applied to the incompressible Navier-Stokes equations. In [30], this model is extended to landslide-induced tsunamis modeling. They developed a model where the landslide equations are averaged in a slope-oriented coordinate system and a Coulomb mixture rheology is used for the landslide material following [22]. The coupling between water and granular layer is made through kinematic boundary conditions and continuity of the normal stresses while friction between water and granular mass is neglected. The resulting model improves the approximation of the landslide evolution even though the model does not account for a real bottom topography. As the authors pointed out in their paper, this is a key point to explain the lack of precision in the landslide deposit prediction. The interaction between both layers is made by an interpolation procedure at the interface between the two considered grids. This model has been compared to experimental data in [16].

We would like to mention here additional effects that are important for the correct simulation of the whole system, such as the effect of the interstitial fluid in the granular mass, the way the interaction between the water and the granular mass is tackled and how the real bottom topography is taken into account. Since the avalanche is composed of a granular material filled with water, submarine avalanches models should account for the influence of the interstitial fluid on the flow as done for debris flows modeling [21], [23], [15]. Other important characteristics of this geophysical problem is the interaction between water and granular material, that is mainly defined by a drag force at the interface and by pressure gradient terms. Finally, bottom variations have to be accurately considered as they mainly drive landslide trajectory and velocity, thus playing a key role on the generated tsunami. Note that in shallow water depth-averaged models, appropriate for tsunami propagation, only the horizontal velocity is kept in the model whereas the flow velocity for landslides is mainly aligned along the slope. The shallow approximation and depth-averaging should be done in the vertical direction for water wave models and in the direction perpendicular to the topography for avalanche models. Hence, we propose to use a Cartesian coordinate system for the water model and a coordinate system related to bottom variations for the avalanche. A depth-averaged Savage-Hutter model over a general bottom curvature was introduced in [5] for dry granular avalanches. This model also takes into account erosion processes that makes the basal coordinate system (local coordinates) time dependent. Thanks to the curvature terms added to the system it is possible to prove an energy balance. This approximation is also used in [28] together with an analysis of different rheologies, the extension to a two dimensional system in [27] and to debris flows problems in [20]. The extension of the general bottom curvature to submarine avalanches is developed in [14] by taking into account a two-phase granular mixture without including erosion effects. Numerical results based on this model are compared with laboratory experiment in [3] and to a Boussinesq-type model considering a rigid slide. In a weaker attempt to account for the real topography, tilted coordinates instead of Cartesian coordinates are commonly used in landslide modeling [21], [23], [6], [7].

The first objective of the present work is to develop an averaged bilayer model for submarine avalanches and tsunamis that takes into account a double coordinate system: Cartesian coordinates for the water layer and a new coordinate system related to the fixed bottom topography for the granular mass. This new system is based on tilted coordinates on a reference plane whose direction follows the main topography variation. This allows us to describe easily the relation between the coordinate systems, providing an exact coupling procedure. Both, the granular and water layers, are assumed to be immiscible and their interaction is described through a drag force and a pressure gradient. It can be viewed as an alternative to landslide Savage-Hutter models proposed in [14], [18] and [30] in order to better reproduce the water-granular dynamics and interaction. The second objective is to propose a well-balanced finite volume method for the proposed model that preserves exactly two kind of stationary solutions. Note that the proposed model degenerates to the two-layer shallow water system in Cartesian coordinates by considering a horizontal reference plane. Then, the numerical method can also be applied to classical bilayer shallow water systems. We do not focus only on the water and material at rest solution; the proposed numerical method also preserves exactly granular masses at rest even if the fluid is still moving. This property is reached by a generalization of the hydrostatic reconstruction for the proposed two-layer model.

The paper is organized as follows: Section 2 presents the coordinate systems, derives the proposed model and shows some of its properties. A numerical method is described in Section 3 to discretize the proposed model. Section 4 describes three numerical tests. The first one is designed to study the influence of the choice of the slope of the reference plane for dry, submerged and aerial initial landslides. In particular, the difference with the case in which only Cartesian coordinates are used is shown. In the second test, the system is written in a non dimensional form in terms of the characteristic height and length of the submarine avalanche and the water and avalanche depths. The influence of two aspect ratios on the spreading of the avalanche and on the maximum water wave at the free surface is analyzed. The last test is designed to highlight the robustness of the numerical method in the presence of wet/dry areas by simulating the collapse of an initially aerial avalanche entering the fluid layer. In Appendix A the proof of the energy balance associated to the model and a bound of the associated eigenvalues is presented.