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

A two-layer shallow flow model with two axes of integration, well-balanced discretization and application to submarine avalanches J.M. Delgado-S´anchez∗ . , F. Bouchut† , E.D. Fern´andez-Nieto∗, A. Mangeney‡ . , G. Narbona-Reina∗ 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. Depthaveraging 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 kind 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. Keywords and phrases. Submarine avalanches; tsunami propagation; bilayer shallow water model; hydrostatic reconstruction; well-balanced finite volume method. ∗Departamento de Matem´atica Aplicada I. ETS Arquitectura, Universidad de Sevilla. Avda. Reina Mercedes 2, 41012, Sevilla, Spain. [email protected], [email protected], [email protected] †Universit´e Paris-Est, Laboratoire d’Analyse et de Math´ematiques Appliqu´ees (UMR 8050), CNRS, UPEM, UPEC, F-77454, Marne-la-Vall´ee, France [email protected] ‡Institut de Physique du Globe de Paris, Seismology team, University Paris-Diderot, Sorbonne Paris Cit´e, 75238, Paris, France. [email protected] 1 1 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 modelling 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 nonhydrostatic 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 2 models, the modelling 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 Figure 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 depthaveraged 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 modelling. 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 modelling [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 3 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 modelling [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 analysed. 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. 4 Figure 1: Sketch of the location of the generated wave. The point x0gives the correct position but local coordinates would give ˜xfor the wave emplacement. 2 Proposed model In this section, we present the derivation of a two-layer model of Savage-Hutter type to study submarine avalanches and generated tsunamis. The domain is made of two layers over a fixed topography. The first one is composed of an homogeneous inviscid fluid and the second one is the granular layer. In what follows, variables corresponding to the fluid layer will be denoted with subindex 1 and those to the granular layer with subindex 2. M Figure 2: Coordinate systems. As discussed above, we consider two different coordinate systems to model the fluid and granular layers. These coordinate systems are chosen according to the physical requirements of the typical flow occurring in each layer. Thus we consider Cartesian coordinates, denoted by ~x = (x, z), for the fluid 5 Figure 3: Notation for the layers at time t layer whereas a particular system of local coordinates is introduced for the granular layer as follows. We define the reference plane determined by a fixed point (¯x, ¯z) and its slope tan θ. For every point M in the grain layer with Cartesian coordinates ~x = (x, z), we introduce the coordinates ~ X= (X, Z), X being the abscissa of the orthogonal projection of Monto the reference plane and Zbeing the distance from Mto the reference plane (see Figure 2). Notice that the local coordinates of the point Mare given by (X/cos θ, Z) when using this notation. The relation between Cartesian coordinates ~x and ~ X coordinates is ~x = (X−Zsin θ, b(X) + Zcos θ), where b(X) is the level of the reference plane: b(X) = ¯z+ tan θ(X−¯x).(1) The fixed topography is defined by ˜ b(X) as the distance, with sign, from the bottom level to the reference plane. Notice that the bottom level is thus given by b(X) + ˜ b(X)cos θ. We will denote by ∇~a and div~a the gradient and divergence operator respect to Cartesian coordinates for ~a =~x and to ~ Xcoordinates for ~a =~ X, whereas b ∂t=∂t|~x and ∂t=∂t|~ X. 2.1 Governing equations We denote the velocity field of the layers by ~ui= (ui, wi) for i= 1,2, where uirepresent the horizontal velocities and withe vertical ones. The densities are denoted by ρiand are assumed to be constant for both layers. The model that governs the physical system is given by the Euler equations, that in Cartesian coordinates read as div~x ~ui= 0,(2) ρib ∂t(~ui) + ρi~ui·∇~x ~ui+ div~x Pi=−ρi∇~x(~g ·~x),(3) 6 for i= 1,2, where ~g = (0, g) and gis the gravity constant. The pressure tensor Piis defined as the negative stress tensor, Pi=pxx ipxz i pzx ipzz i,(4) with pxz i=pzx i. We consider exactly the previous equations for the fluid layer (i= 1). For the grain layer (i= 2), we perform the change of variables ~x 7→ ~ Xwhich enables to reformulate system (2)–(4) in a more suitable form. With the aim of clarity we detail the main operators appearing in the change of variables procedure. Thus, for a given vector field ~ Φ, recall that J∇~x ·~ Φ = ∇~ X·(JA~ Φ),(5) where A−1=∇~ X~x and J= det A−1; for any scalar function f(~x) we have ∇~xf= (∇~x ~ X)T∇~ Xf. where (·)Tdenotes the transpose of the corresponding matrix (see [5, Lemma 2.1]). Since (∇~x ~ X)−1=∇~ X~x =1−sin θ tan θcos θ then, ∇~xf=cos2θ−sin θ sin θcos θcos θ∇~ Xf. (6) Note that the change of variables is then valid when θ6=π/2, because J= 1/cos θ. We also introduce the following notation for the grain layer: U2 W2=cos θsin θ −sin θcos θu2 w2,(7) that is, U2(respectively, W2) is the velocity of the second layer parallel (respectively, perpendicular) to the reference plane (see Figure 3). According to this, the following system of partial differential equations is obtained to describe the granular layer:    ∂XU2cos θ+∂ZW2= 0, ρ2∂tU2+ρ2cos θ∂X(U2 2) + ρ2∂Z(W2U2) + ρ2cos θ∂X(~g ·~x) = −cos θ∂XPxx 2−∂ZPzx 2 ρ2∂tW2+ρ2cos θ∂X(U2W2) + ρ2∂Z(W2 2) + ρ2∂Z(~g ·~x) = −cos θ∂XPxz 2−∂ZPzz 2. ,(8) where the pressure tensor P2is defined by P2=cos θsin θ −sin θcos θP2cos θ−sin θ sin θcos θ=Pxx 2Pxz 2 Pzx 2Pzz 2.(9) We introduce analogously the tensor P1that will be used for the boundary conditions at the interface. We also assume that Pxx i=Pzz ifor i= 1,2. 7 Let us now introduce the kinematics and boundary conditions in order to complete the system. Firstly, we define the thicknesses and boundaries of the physical system depicted in Figure 3. We denote by h1(x, t) the thickness of the fluid layer at the position xalong the horizontal axis at time t, while H2(X, t) is the thickness of the grain layer perpendicular to the reference plane. The bottom is defined by z=B(x), with B(x) = b(X)+˜ b(X)cos θor by Z=˜ b(X) when using ~ Xcoordinates. The interface between the grain and fluid layer is defined by z=I(x, t), with I(x, t) = B(x) + H2(X, t) cos θ=b(X)+(˜ b(X) + H2(X, t)) cos θ, while for ~ Xcoordinates, the interface is Z=˜ b(X) + H2(X, t). Finally the water free surface is defined by the equation z=S(x, t) = I(x, t) + h1(x, t). At the three boundaries, free surface, interface between layers and bottom, the density is discontinuous, so the kinematic conditions are directly deduced from the Rankine-Hugoniot condition for the mass conservation, that is written as b ∂tρ+ div~x(ρ~u) = 0. Writing down this condition through the interfaces we obtain the general boundary condition [ρ],[ρ~u]·~ N= 0,(10) where [·] denotes the jump of a quantity through the interface and ~ Nis the time-space normal to the corresponding interface. In our system, the density ρtakes the value 0 up to the free surface, ρ1in the fluid layer and ρ2in the grain layer. The velocity ~u takes the value ~u1in the fluid layer, ~u2in the grain layer and 0 under the bottom interface. At the free surface, defined by z−S(x, t) = 0, we have ~ N= (−b ∂tS, −∂xS, 1), so from (10) we obtain the kinematic condition b ∂tS+u1∂xS−w1= 0 at z=S(x, t).(11) The normal vector to the bottom, defined by z−B(x) = 0, is ~ N= (0,−∂xB, 1). So the kinematic condition at the bottom level is given from (10) by the usual no-penetration condition, u2∂xB−w2= 0 at z=B(x). In order to use this condition for layer 2, it must be written in ~ Xcoordinates. We observe that z−B=Z−˜ bat the bottom level, so the previous equation can be written as u2∂x(z−B)−w2∂z(z−B) = 0 or, equivalently, u2∂x(Z−˜ b)−w2∂z(Z−˜ b) = 0. Now using (6) and (7) we get W2=U2cos θ ∂X˜ bat Z=˜ b(X).(12) Finally, the interface between the fluid and the grain layers is defined by z=I(x, t), so ~ N= (−b ∂tI, −∂xI, 1). Then, from (10) we obtain the following identity, ρ1(b ∂tI+u1∂xI−w1) = ρ2(b ∂tI+u2∂xI−w2). 8 Since there is no exchange of mass between the fluid and granular layers, the kinematic conditions at the interface are given by setting to zero the left and right hand side of the previous equality. This corresponds to the transport condition of the points at the interface. Then we have b ∂tI+u1∂xI=w1at z=I(x, t),(13) b ∂tI+u2∂xI=w2at z=I(x, t).(14) As in the bottom case, we rewrite the kinematic condition for layer 2 using ~ Xcoordinates. For this, notice that (14) is equivalent to b ∂t(z−I) + u2∂x(z−I) + w2∂z(z−I) = 0 at Z=˜ b(X) + H2(X, t). Note that z−I=Z−(˜ b+H2) at the interface between layers. So, in view of (6) and (7), we conclude ∂tH2+U2cos θ ∂X(˜ b+H2) = W2at Z=˜ b(X) + H2(X, t).(15) Now we set the boundary conditions. We introduce the spacial unitary normal vector to each boundary with positive vertical component, denoted by ~nSfor the water free surface, ~nIfor the interface and ~nBfor the bottom, respectively. Then, at z=S(x, t), we assume the free surface tension condition given by P1~nS= 0 at z=S(x, t).(16) At the interface z=I(x, t), we impose the continuity of the normal component of the tensors and a friction force between the fluid and grain layers denoted by F, to be specified later: (P1~nI)N= (P2~nI)Nat z=I(x, t),(17) (P1~nI)T= (P2~nI)T=Fat z=I(x, t),(18) where the subscript (·)N(respectively, (·)T) denotes the normal (respectively, the tangent) component of a vector. Finally, we consider the friction at the bottom z=B(x) as follows, (P2~nB)T=−~nB·(P2−P1)~nBU2 |U2|µat z=B(x),(19) µbeing the friction coefficient that will be given by a friction law. For instance, the Coulomb friction law gives µ= tan δ0, where δ0is the friction angle of the granular material in layer 2. 2.2 Asymptotic analysis and averaging process As commonly done for the shallow averaged models, we consider a dimensional analysis of the set of equations defined by (2)–(3) for the first layer, and (8) for the second one together with the kinematic conditions (11), (12), (13), (15) and the boundary conditions (16)-(19). We introduce Hand Las the characteristic height and length of the domain. We assume a shallow domain, so ε=H/L is supposed to be small. We will denote with (∗) the non-dimensional variables, then (x, z)=(Lx∗, Hz∗),(X, Z)=(LX∗, HZ∗), t = (L/g)1/2t∗, 9 3 Numerical method In this section the solver considered for the numerical tests is presented. The friction term between layers can be discretized semi-implicity, so we describe the numerical treatment of the other components of the model in what follows. Firstly, let us remark that model (35) can be rewritten as                        ∂t(H1J) + ∂XH1(U1−∂X(H2U2) sin θcos θ)= 0, ∂t(H1U1J) + ∂XH1U1(U1−∂X(H2U2) sin θcos θ)+gH1∂XH1+b+ (˜ b+H2) cos θ= 0, ∂tH2+ cos θ∂X(H2U2)=0, ∂t(H2U2) + cos θ∂XH2U2 2+gH2cos θ∂X(1 −r)(b+ (˜ b+H2) cos θ) + ξr(H1+ (b+ (˜ b+H2) cos θ))= =−gsgn(U2)H2cos θ(1 −r)µ, where ξ= 1. We have introduced this parameter ξbecause a simplification of the model, that considering the rigid lid hypothesis on the momentum equation of the second layer, corresponds to set ξ= 0. This hypothesis can be considered for the case of submarine avalanches because the gradient of the water free surface, that is, ∂X(H1+b+ (˜ b+H2) cos θ) does not have a great influence on the movement of the granular layer. By another way, the advantage of considering this hypothesis is that the eigenvalues of the model are always real when ξ= 0. For example, complex eigenvalues are typical in coupled submarine avalanches models in case of an aerial avalanche entering the fluid. In such a case, some numerical treatment of complex eigenvalues must be considered (see [12]). By setting ξ= 0 we avoid this difficulty and very similar results are obtained. Then, for the numerical tests we consider ξ= 0. By simplicity, for a domain [a, b] we consider a fixed grid with step ∆x, with xi=a+ (i−1/2)∆x, xi+1/2=a+i∆x. Let us denote the vector of unknows by W= (H1,Q1, H2, Q2)T= (H1J, H1JU1, H2, H2U2)T. A HLL-PVM-1U method is considered (see [11]), that can be defined as follows: Wn+1 i=Wn i−∆t ∆x(˜ Fn i+1/2−˜ Fn i−1/2)−∆t ∆xSn i,(42) where Wn+1 idenotes the approximation of the average of the vector solution at volume (xi−1/2, xi+1/2) and time tn+1 =tn+ ∆t. In what follows the definition of the numerical flux functions ˜ Fi+1/2and the approximation of source and non-conservative terms Siis presented. With the purpose of simplicity, we omitt the superindex nin the description. Let us denote by Jian approximation of J; we can set Ji= max 1−sin θ(˜ bi+1 +H2,i+1 −˜ bi−1−H2,i−1)/(2∆x),tol, where tol is a minimum tolerance that allows to consider situations where the model is near the limit of application. It can be related to the maximum slope of the granular free surface. For example, it can be set as tol = tan(85◦). 16 An hydrostatic reconstruction is considered for the fluid and the granular layers (see [4]). It allows us to treat wet/dry fronts and to discretized the Coulomb friction term and to obtain a well-balanced method for solutions at rest for the fluid and the granular layers. Firstly, at a given interface xi+1/2, the following left and right reconstruction are defined. For the fluid layer we define H− 1,i+1/2= max(H1,i −max(∆z1,0),0), H+ 1,i+1/2= max(H1,i+1 −max(−∆z1,0),0), where ∆z1= (˜ bi+1 −˜ bi+Hn 2,i+1 −Hn 2,i) cos θ+ ∆xtan θ. Then, H± 1,i+1/2=H± 1,i+1/2Ji+1 2±1 2and Q± 1,i+1/2=H± 1,i+1/2U1,i+1 2±1 2. Secondly, the reconstruction of the variables associated to the granular layer is done. Let us denote by fi+1/2the approximation of the Coulomb friction term. It is defined as follows fi+1/2=       fcµU∗ 2,i+1/2 |U∗ 2,i+1/2|if |U∗ 2,i+1/2|> fcµ∆t, U∗ 2,i+1/2 ∆totherwise, where fc=gcos θ(1 −r),and U∗ 2,i+1/2=U2,i +U2,i+1 2−fc ∆t ∆x∆z1+ξ r 1−r(H+ 1,i+12 −H− 1,i+1/2). This approximation of the Coulomb friction term is considered in the definition of the hydrostatic reconstruction of H2, defined as follows: H− 2,i+1/2= max(H2,i −max(∆z2,0),0), H+ 2,i+1/2= max(H2,i+1 −max(−∆z2,0),0) where ∆z2= (˜ bi+1 −˜ bi) + ∆x cos θfi+1/2 fc + tan θ. Finally, the solid discharge reconstruction is defined by Q± 2,i+1/2=H± 2,i+1/2U2,i+1 2±1 2 and the reconstructed vector of unknowns will be denoted as W± i+1/2= (H± 1,i+1/2,Q± 1,i+1/2, H± 2,i+1/2, Q± 2,i+1/2)T. By using these definitions of the state reconstructions we consider the following definition of the numerical flux: ˜ Fi+1/2=1 2F(W+ i+1/2) + F(W− i+1/2) −1 2α0,i+1/2(W+ i+1/2−W− i+1/2) + α1,i+1/2(F(W+ i+1/2)−F(W− i+1/2) + Si+1/2), 17 where F(W) =     H1U1 H1U2 1 H2U2cos θ H2U2 2cos θ     and Si+1/2=       0 gH1,i +H1,i+1 2(H+ 1,i+1/2−H− 1,i+1/2) 0 g(1 −r) cos θH2,i +H2,i+1 2(1 −r) cos θ(H+ 2,i+1/2−H− 2,i+1/2) + ξ r(H+ 1,i+1/2−H− 1,i+1/2)        . Moreover, Si=1 2Si−1/2+Si+1/2+1 2Ci−1/2+Ci+1/2, where Ci+1/2=            −H1,i+1 +H1,i 2 Q2,i+1 −Q2,i ∆xcos θsin θ −H1,i+1U1,i+1 +H1,iU1,i 2 Q2,i+1 −Q2,i ∆xcos θsin θ 0 0            . This definition corresponds to the application of the hydrostatic reconstruction to the HLL-PVM-1U method introduced in [11]. Moreover, for the hydrostatic reconstruction, based in a reconstruction in terms of the gradient pressure appearing in the energy of the model, we have used the following equality: ∂Xb+ (H2+˜ b) cos θ+rH1= (1 −r)∂Xb+ (H2+˜ b) cos θ+r∂Xb+ (H2+˜ b) cos θ+H1. Note that the last term is the gradient of the water free surface, defined in terms of the hydrostatic reconstruction of the fluid layer. The first term is discretized by the definition of the hydrostatic reconstruction for the second layer. The first and second components of Ci+1/2are discretizations of the terms H1∂tH2sin θand H1U1∂tH2sin θ, respectively. Finally, the HLL-PVM-1U method is defined by the following two coefficients, α0,i+1/2= (SR,i+1/2|SL,i+1/2|−SL,i+1/2|SR,i+1/2|)/(SR,i+1/2−SL,i+1/2), α1,i+1/2= (|SR,i+1/2|−|SL,i+1/2|)/(SR,i+1/2−SL,i+1/2), where the characteristic velocities SL,i+1/2and SR,i+1/2are defined as an approximation of the bounds of the eigenvalues given by (38), namely, SL,i+1/2=¯ Ui−Ψi Ji , SR,i+1/2=¯ Ui+1 −Ψi+1 Ji+1 , 18 where ¯ Ui:= U1,i +Jn iU2,i cos θ 2and Ψi=r3 2 1 2(U1,i −JiU2,i cos θ)2+gH2,i(cos θ)3J2 i+gH1,i!1/2 . This numerical method satisfies the following well-balanced properties, corresponding to two kind of steady solutions. The first one corresponds to the classical one consisting of water at rest solution without movement of the granular layer, defined by U1= 0, U2= 0, ∂XH1+b+ (H2+˜ b) cos θ= 0,∂Xb+ (H2+˜ b) cos θ≤µ. The second one corresponds to the case in which the fluid is not at rest but the granular layer does not move. That occurs when U2= 0 and ∂Xb+ (H2+˜ b) cos θ+ξr 1−r(H1+b+ (H2+˜ b) cos θ)≤µ. Note that this second solution depends on the rigid lid hypothesis. Without this hypothesis, by setting ξ= 1, the gradient of the water free surface can produce the movement of the granular layer in some cases. Theorem 2. The numerical method defined in this section preserves: (i) exactly steady state solutions corresponding to constant water free surface with granular free surface with a slope smaller than µ, if the initial condition is defined veyfying U0 i,1=U0 2,i = 0,(˜ bi+1 −˜ bi+H0 2,i+1 −H0 2,i) cos θ+ ∆xtan θ+H1,i+1 −H1,i = 0, 1 ∆x(˜ bi+1 −˜ bi+H0 2,i+1 −H0 2,i) cos θ+ ∆xtan θ≤µ. (ii) exactly steady states of the granular layer, even for the case in which the fluid is not at rest, if Un 2= 0,1 ∆x(˜ bi+1 −˜ bi+Hn 2,i+1 −Hn 2,i) cos θ+ ∆xtan θ +ξ r 1−r((˜ bi+1 −˜ bi+Hn 2,i+1 −Hn 2,i) cos θ+ ∆xtan θ+H1,i+1 −H1,i)≤µ. The proof is straigthforward: it is enough to see that under the above conditions we obtain that ˜ Fi+1/2= 0 and Si= 0 in the case (i); for the case (ii). this equality is deduced for the third and fourth components of ˜ Fi+1/2and Si. We cannot prove the positivity of the numerical method, even though the numerical tests have not shown any problem of negative heights. A way to obtain the positivity could be to develop a Suliciu type solver (see [4]). 19 4 Numerical tests The numerical tests have been done with a numerical method based on a hydrostatic reconstruction with a PVM-1U-HLL solver (see [4] and [11]) defined in terms of the bounds of the eigenvalues given by (38). In this section, we present three numerical tests. In the first test, the influence of the choice of the slope of the reference plane is analysed for the cases of dry, submerged and aerial avalanches. In the second test, the system is written in a non dimensional form in order to study the effect of two aspect ratios on the spreading and deposit of the avalanche and on the maximum water wave generated at the free surface. Finally, the interaction of an initially aerial avalanche with the fluid is shown in the third test to check that the model properly deals with wet/dry fronts. 4.1 Test 1: Influence of the slope of the reference plane The goal of this test is to show the flow dynamics obtained for each coordinate system in the cases they are designed for and, subsequently, to quantify the error made in other situations. As initial condition, we consider a granular mass released on an inclined plane of slope α(see Figure 4 for sketch and notation) that will flow down-slope. The initial constant upper surface of the granular layer is denoted by hs0and its vertical thickness is maximum at point x0. We simulate four cases: (A) a dry granular mass without a fluid layer over an inclined plane, (B) a dry granular mass over two planes (inclined/horizontal), (C) a submerged granular mass inside a fluid layer with initial free surface height equal to 1 m, and (D) an initially sub-aerial granular mass that will enter a fluid layer of initial free surface height equal to 0.6 m. For cases (C) and (D) the bottom is defined as two planes, inclined/horizontal, (see Figure 4). The experiments have been carried out for two different slopes: the slope α= 35◦, with x0= 1 m and hs0= 0.96 m; and the slope α= 45◦, with x0= 0.6 m and hs0= 0.9 m. The parameters involved in the model have been chosen as follows: the densities of the fluid and granular layers are ρ1= 1000kg/m3and ρ2= 1900kg/m3, respectively, thus providing r= 0.5263. The friction coefficient between the layers is γ= 0.001 and, at the bottom, the Coulomb friction law is used with a friction coefficient µ= tan 28◦, except in case A in which µ= tan 38◦(respectively, µ= tan 48◦) when α= 35◦(respectively, α= 45◦). For the coordinate system, we test different values for the slope θof the reference plane: θ=α, that corresponds to the slope of the inclined plane, θ= 0◦(i.e., Cartesian coordinates), that corresponds to the slope of the horizontal plane, and an intermediate value θ=α/2. The results are obtained by using ∆x= 0.005 m in the discretization of the domain [-2 m,2 m]. The CFL condition is 0.5. Let us first investigate the effect of the chosen coordinate system on dry granular flows (Figures 4a- 4b) by comparing our numerical results to those of the code SHALTOP, specially designed to accurately describe topography variations in shallow depth-averaged models [8, 9, 33, 32, 10]. This model can be written in a Cartesian coordinate system but the equations are derived by applying the shallow approximation and depth-averaging in the direction normal to the slope. As a result, it describes flows over an arbitrary topography and contains the whole curvature tensor of the terrain. This comparison is also made for the initial submerged test (Figure 4c) by considering the reduced gravity approximation, 20 (a) Case A. Dry granular mass over an inclined plane. (b) Case B. Dry granular mass over two planes, inclined/horizontal. (c) Case C. Submerged granular mass. 0.6 (d) Case D. Aerial granular mass. Figure 4: Test 1. Initial conditions for the studied cases. i.e., by replacing gwith (1 −r)gin SHALTOP as explained in [10]. Figures 5, 6, 7 and 8 show the final deposit for each case. The difference between the coordinate systems is obviously higher for α= 45◦than for α= 35◦. Figures 5, 6, 7 and 8 show the final deposit for each case. The difference between the coordinate systems is obviously higher for α= 45◦than for α= 35◦. Overall, the deposit obtained for θ=αgoes further, followed by the solution obtained with Cartesian coordinates and at last by the intermediate slope solution. Interestingly, the behavior is not the same during the flow on the inclined plane as observed on Figures 6a, 6b, 7a and 7b. Contrary to the deposit, the faster solution is obtained when using the Cartesian coordinates (dashed green line) while the solution for θ=α(dashed blue line) remains backward. In particular, for α= 45◦, the granular front already reached the horizontal bottom at these intermediate times for the Cartesian coordinates solution, as opposed to the local coordinates solutions, whose front positions are still on the inclined plane. The smaller velocity obtained for the local coordinates system compared to the Cartesian coordinates system for flows on inclined planes is due to the difference in the driving force related to pressure gradients and Coulomb friction. Indeed, the terms gH2cos θ∂X(b+ (˜ b+H2)cos θ+rH1) and gsgn(U2)H2cos θ(1 −r)µare not multiplied by cos θin the Cartesian coordinate system. In particular we can see clearly this effect in Figure 5 corresponding to Case A. In this figures we also show the deposits from the local point of view, that is perpendicular to the bottom plane. We can observe that the solution of the proposed model for θ=αcoincides with the solution given by SHALTOP. A small difference can be seen at the front position due to a different numerical treatment of the wet/dry transition. Note also that the plotted solution of SHALTOP is obtained by a longitudinal cut of the 2D provided solution. For the Cartesian model the front position 21 is overestimated and contrarily the solution obtained with θ=α/2 is underestimated, as expected. The results obtained are not surprising regarding the cases θ= 0◦and θ=α. The idea behind the local coordinates is to better capture the evolution of the material when it travels over the inclined plane of slope α, since the computed velocity is parallel to the reference plane. Thus, the solution provided by the model with θ=αwill be suitable for the evolution of the mass along the inclined plane, but inappropriate once the material reaches the horizontal part. Analogously, the model in Cartesian coordinates works properly on the horizontal part of the bottom but does not provide a good approximation for the inclined evolution of the granular mass because the computed velocity is horizontal. If we look at the comparison with the reference solution provided by SHALTOP, we can see that the flow on the inclined plane is well reproduced with θ=αas expected, contrary to θ=α/2 and θ= 0◦. However, the best approximation of the deposit is obtained with θ=α/2. In particular, when α= 35◦, the solution obtained for θ=α/2 is the only one that retains mass at the rear part of the avalanche over the inclined plane (see Figures 6a, 7a, 8). On the contrary, when α= 45◦, none of the solution reproduce the retained mass simulated by SHALTOP (see Figures 6b, 7b), although the closer one is again the solution obtained for θ=α/2. The change in granular flow dynamics obtained with these different coordinate systems has a strong impact on the generated water waves as shown for the submerged case in Figures 7a, 7b and for the aerial case in Figure 8. In the submerged case, the difference is again larger when the slope is steeper. In the aerial case, there is a big difference between the solutions obtained with the Cartesian and the proposed local coordinate system. One of the most interesting conclusion of these results is that the solution obtained for a reference plane with slope θbetween 0◦and the maximum inclination of the bottom gives the most reasonable solution globally (from the inclined to the horizontal part of the bottom) when compared to SHALTOP. In particular, it is able to keep part of the mass on the inclined plane while other solutions (θ= 0◦or θ=α) could not. 22 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 -0.2 0 0.2 0.4 0.6 0.8 1 0.8 1 1.2 0.6 0.7 0.8 0.9 1 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 (a) α= 35◦. -0.2 0 0.2 0.4 0.6 0.8 1 1.2 -0.2 0 0.2 0.4 0.6 0.8 0.6 0.7 0.8 0.9 0.65 0.7 0.75 0.8 0.85 0.9 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 0 0.02 0.04 0.06 0.08 0.1 0.12 (b) α= 45◦ Figure 5: Test 1. Case A: dry granular mass over an inclined plane. Figures on the left show the deposits and zoom of the rear part. In figures on the right, the thickness of the granular mass is measured perpendicular to the bottom for each case. 23 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.1 0.2 0.3 0.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 -0.2 0 0.2 0.4 0.6 0.8 (a) α= 35◦. Left, evolution of the avalanches on the inclined plane (t= 0.4 s). Right, deposits. -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 -1.5 -1 -0.5 0 -0.2 0 0.2 0.4 0.6 0.8 1 (b) α= 45◦. Left, evolution of the avalanches on the inclined plane (t= 0.3 s). Right, deposits. Figure 6: Test 1. Case B: dry granular mass over inclined/horizontal planes. 24 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 0 0.2 0.4 0.6 (a) α= 35◦. Left, evolution of the avalanches on the inclined plane (t= 0.4 s). Right, deposits. -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0 0.2 0.4 0.6 0.8 (b) α= 45◦. Left, evolution of the avalanches on the inclined plane (t= 0.5 s). Right, deposits. Figure 7: Test 1. Case C: submmerged granular mass over inclined/horizontal planes. 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Figure 8: Test 1. Case D: aerial granular mass, α= 35◦. Left, evolution of the avalanches on the inclined plane (t= 0.3 s). Right, deposits. 25 4.3 Test 3: Aerial avalanche In this test we consider the case of an initially sub-aerial granular avalanche that reaches a fluid area. Moreover we added a bump on the bottom to produce non trivial solutions, involving in particular several wet/dry transitions between the avalanche and the fluid. A reference plane with an inclination angle θ= 10ois considered. The computation is performed over the interval [-3 m ,3 m] with a space step ∆X= 3/250 m and a CFL parameter equals to 0.5. The fixed bottom topography is determined by the function ˜ b(X) = 0.1e−X2/0.01, that is a nearly plane bed with a bump in the middle zone. As in Test 1, we have also set the density ratio to r= 0.5263, the friction coefficient between the phases to γ= 0.001 and a Coulomb friction law for the granular layer with a friction coefficient µ= tan 28◦. The initial condition is set to U1(x, 0) = U2(X, 0) = 0 m/s, and h2(X, 0) = (1 m,if 1.3 m < X < 1.7 m 0 m,otherwise , h1(x, 0) = max {0.2−I(x, 0),0}m, where I(x, 0) is the ordinate of the interface between the grain and fluid layer at the position x(see figure 14a). −2−1 0 1 2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 x z Fluid surface Granular surface Bottom (a) Initial condition −2−1 0 1 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 x z (b) Steady solution Figure 14: Test 3. Left: Initial condition. Right: final steady solution. The final steady state is complex as shown on Figure 14b. The granular deposit has a bumpy shape with an area filled with fluid over the granular layer, isolated from the rest of the fluid. Figure 15 shows the time evolution of the system in the early stage. Water waves generated by the landslide at the free surface are observed, firstly, when the avalanche enters the water and, later, when it collides with the bump (Figure 15(b), (c) and (d)). As a consequence of the bump, part the granular layer returns back (see Figure 15(e) and (f)) before it stops on the slope and forms a deposit. 32 −2−1 0 1 2 −2 −1 0 1 2 x z Fluid surface Granular surface Bottom (a) t= 0.1 s −2−1 0 1 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 x z (b) t= 0.3 s −2−1 0 1 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 x z (c) t= 0.5 s −2−1 0 1 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 x z (d) t= 0.6 s −2−1 0 1 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 x z (e) t= 1.1 s −2−1 0 1 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 x z (f) t= 1.5 s Figure 15: Evolution of the system in the early stage (zoom). 33 5 Conclusions A two-layer depth-averaged model has been proposed by considering two different axes of integration in order to describe submarine landslides over complex topographies and the generated perturbation of the water free surface. For the fluid layer, the shallow approximation and depth-averaging are performed in a Cartesian (i.e., horizontal/vertical) coordinate system while for the granular layer a reference plane following the averaged topography is chosen, called local coordinate system. This makes it possible to account for the different characteristics of (i) the fluid flow, mainly horizontal, and (ii) the avalanche flow that is essentially tangent to the topography. Furthermore, it allows to properly locate the water surface perturbation on top of the avalanche. In particular, when considering a model completely written in local coordinates, the location of the perturbation on the water free surface depends artificially on the choice of the angle of the reference plane. The resulting model has been rewritten in terms of only one horizontal system of coordinates. The proposed model verifies an exact energy balance. A well-balanced finite volume method is proposed to solve the equations. It preserves exactly the granular mass at rest, even if the fluid is still moving. This is an additional well-balanced property that has not been usually studied for coupled two-layer models. The classical water at rest solution is also exactly preserved. This method can also be applied for the two-layer Cartesian system. We have simulated dry and submarine avalanches as well as initially sub-aerial avalanches entering the fluid layer and the generated perturbations at the water free surface. These simulations have been performed on (i) an inclined bed and (ii) a two-slope bottom starting with an inclined plane reaching a horizontal plane. We have investigated the influence of the reference frame angle and compared our results to those obtained with a Cartesian coordinate system and with the code SHALTOP designed to simulate dry granular flows over complex topography. We have found that the simulated dynamics and deposits as well as the generated water surface perturbations are very sensitive to the chosen coordinate system. In particular, significant errors are obtained when using a Cartesian coordinate system for flows on an inclined bed, while a sloping coordinate system prevents reproducing the mass spreading on horizontal planes. Overall, the use of an inappropriate coordinate system always overestimates the runout distance of granular flows. For a two-slope bottom, choosing the average slope for the reference frame provides the best compromise. By varying the initial geometry of the granular mass released on an inclined plane as well as its depth within the water layer, we quantified the error made when using Cartesian coordinate systems instead of a local coordinate system following the plane. The slope angle of the plane has been also varied. By using non-dimensional variables, we proved that the model can be written only in terms of three aspect ratios, corresponding to ratios between the characteristic heights of the fluid, the avalanche, the bottom, and the characteristic length of the initial profile of the granular layer. Our results show that even for moderate slopes (up to 30◦), strong relative errors on the avalanche dynamics and deposits (up to 60%) and on the generated water waves (up to 120%) are made when using Cartesian coordinate systems for both the landslide and water motion as done in, e.g., [15], instead of the two-axis of integration as proposed here. In a forthcoming work we will investigate the improvement of the submarine avalanches modelling by including some of the effects commented in the introduction, as the dispersive effects for the water layer in this framework of two-axis of integration. Submarine avalanches produce a high vertical acceleration on the fluid, so non-hydrostatic effects are necessary to introduce. Dispersive effects are also necessary to 34 compare with experimental results, that usually provide the time series of the water surface evolution. Furthermore, since the submarine avalanche can generate a tsunami, we can also use this model to follow its propagation because in areas where there is no granular mass the model reduces to a one layer shallow model. The extension to the two-dimensional case is mandatory if we would like to apply the proposed model to realistic submarine avalanches. In this sense, it is easy to carry out a 2D extension of the model where the inclination of the reference plane follows the x-axis, being parallel to the y-axis. Nevertheless, to chose a reference plane that fits optimally the bathymetry can implies a rotation as well as tilting. In this case, the extension of the model entrains a thorough development that must be carefully done. Acknowledgements This research has been partially supported by the Spanish Government under grant MTM 2015-70490- C2-2-R with the participation of FEDER, by the ANR contract ANR-11-BS01-0016 LANDQUAKES, the USPC PEGES project and the ERC contract ERC-CG-2013-PE10-617472 SLIDEQUAKES. Appendix A Proof of Theorem 1 In this appendix we present the proof of Theorem 1. (i) For clarity in the exposition we recall here the system to be studied for a general friction force F. The specific expression will be solved through the search of a dissipative energy balance, as exposed below. The system reads: ∂t(H1J) + ∂X(H1(U1+∂tH2sin θ)) = 0,(43a) ∂t(H1U1J) + ∂X(H1U1(U1+∂tH2sin θ)) = −gH1∂X(H1+I)−J F ρ1 ,(43b) ∂tH2+ cos θ∂X(H2U2)=0,(43c) (43d)∂t(H2U2) + cos θ∂X(H2U2 2) = −gH2cos θ∂X(rH1+I)−gsgn(U2)H2cos θ(1 −r)tanδ0+F ρ2 , where we use the notation of the interface, I(X, t) = b+ (˜ b+H2) cos θ. We multiply the momentum equation for the layer 1 in (43b) by U1and we use the mass equation (43a). The left hand side gives: U1JH1∂t(U1) + U1∂t(JH1) + H1(U1+∂tH2sin θ)∂X(U1) + U1∂X(H1(U1+∂tH2sin θ)) =JH1∂tU2 1 2+H1(U1+∂tH2sin θ)∂XU2 1 2 =∂t(JH1 U2 1 2) + ∂XH1(U1+∂tH2sin θ)U2 1 2. 35 Now we develop the first term on the right hand side: −gH1∂X(H1+I)U1=−∂X(gH1(H1+I)U1) + g(H1+I)∂X(H1U1) =−∂X(gH1(H1+I)U1)−∂X(g(H1+I)H1∂tH2sin θ) +gH1∂tH2sin θ∂X(H1+I)−gH1∂t(JH1)−gI∂t(JH1) =−∂X(gH1(H1+I)(U1+∂tH2sin θ)) −∂t(gIJH1)−gH2 1∂tJ −gJ∂t(H2 1 2) + gJH1∂tI+gH1∂tH2sin θ∂X(H1+I) =−∂X(gH1(H1+I)(U1+∂tH2sin θ)) −∂t(gIJH1)−∂t(gJH2 1 2)−gH2 1∂tJ +gH2 1 2∂tJ+gJH1∂tH2cos θ+g∂tH2sin θ(∂X(H2 1 2) + H1∂XI) =−∂X(gH1(H1+I)(U1+∂tH2sin θ)) −∂t(gJH1(I+H1 2)) + ∂X(gH2 1 2∂tH2sin θ) −gH2 1 2∂tJ+gJH1∂tH2cos θ−gH2 1 2∂X(∂tH2sin θ) + g∂tH2sin θH1∂XI =−∂tgJH1(I+H1 2)−∂XgH1(H1 2+I)(U1+∂tH2sin θ)) + gH2 1 2U1 +gJH1∂tH2cos θ+g∂tH2sin θH1∂XI | {z } R1 . Notice that using the definition of J, we have ∂tJ=−∂X(∂tH2sin θ), that we used to simplify the last equality. Now we focus on the last term; since ∂Xb= tan θ, it follows: R1=gH1∂tH2cos θ(J+ tan θ∂XI) = gH1∂tH2cos θ(1 + tan2θ) = gH1∂tH2 1 cos θ. Finally, the energy balance for the first layer reads: (44) ∂tJH1U2 1 2+g(H1 2+I)+∂XH1(U1+∂tH2sin θ)U2 1 2+g(H1 2+I)+gH2 1 2U1 =gH1∂tH2 1 cos θ−JU1 F ρ1 . Similarly, to obtain the energy balance for the layer 2, we multiply (43d) by U2and use (43c). The left hand side gives: U2H2∂t(U2) + U2∂t(H2) + H2cos θU2∂X(U2) + U2∂X(H2cos θU2) =H2∂t U2 2 2+H2cos θU2∂X U2 2 2 =∂tH2 U2 2 2+∂XH2cos θU2 U2 2 2. 36 The first term on the right hand side is: −gH2cos θ∂X(rH1+I)U2=−∂X(gH2cos θ(rH1+I)U2) + g(rH1+I)∂X(H2cos θU2) =−∂X(gH2cos θ(rH1+I)U2)−g(rH1+I)∂tH2 =−∂X(gH2cos θ(rH1+I)U2)−grH1∂tH2−g(b+ (˜ b+H2) cos θ)∂tH2 =−∂X(gH2cos θ(rH1+I)U2)−∂tgH2b+˜ bcos θ+H2 2cos θ −grH1∂tH2. Then the energy balance for the layer 2 reads: (45) ∂tH2 U2 2 2+gH2b+˜ bcos θ+H2 2cos θ+∂XH2cos θU2U2 2 2+g(rH1+I) =−grH1∂tH2−g|U2|H2cos θ(1 −r)tanδ0+F ρ2 U2. The total energy balance of the system is found by the combination rcos θ·(44) + (45): ∂trcos θJH1U2 1 2+gH1 2+gI+H2 U2 2 2+gH2b+˜ bcos θ+H2 2cos θ +∂Xrcos θH1(U1+∂tH2sin θ)U2 1 2+g(H1 2+I)+rcos θgH2 1 2U1 +H2cos θU2U2 2 2+g(rH1+I)=−g|U2|H2cos θ(1 −r)tanδ0−(JU1cos θ−U2)F ρ2. So, if we define F=ρ1β(JU1cos θ−U2),(46) we find the inequality given in the Theorem 1. (ii) To prove the second part of Theorem 1 it is enough to apply the result presented in [34] to the characteristic polynomial of the transport matrix of the model. It is possible to rewrite the model is written under the structure ∂tW+A(W)∂XW−∂X(D(W)∂XW) = b, where Ais the transport matrix, Wis the vector of unknowns of the mode, matrix Dis defined by the coefficients related to the second order derivatives, by using that ∂tH2=−cos θ∂X(H2U2) at the first and second equations of the model, and vector bcontains the definition of friction terms. The result presented in [34] states that if all the roots of a polynomial anxn+an−1xn−1+. . . +a0 of degree nare real, then they lie in the range bounded by −an−1 nan±(n−1) a2 n−1 n2a2 n−2an−2 n(n−1)an1/2 . 37 References [1] S. Abadie, D. Morichon, S. Grilli, and S. Glockner. Numerical simulation of waves generated by landslides using a multiple-fluid navier–stokes model. Coastal Engineering, 57(9):779 – 794, 2010. [2] S. M. Abadie, J. C. Harris, S. T. Grilli, and R. Fabre. Numerical modeling of tsunami waves generated by the flank collapse of the cumbre vieja volcano (la palma, canary islands): Tsunami source and near field effects. Journal of Geophysical Research: Oceans, 117(C5), 2012. [3] B. Ataie-Ashtiani and S. Yavari-Ramshe. Numerical simulation of wave generated by landslide incidents in dam reservoirs. Landslides, 8(4):417–432, Dec 2011. [4] F. Bouchut. Nonlinear stability of finite volume methods for hyperbolic conservation laws. Frontiers in Mathematics, 49(1):101 – 140, 2004. [5] F. Bouchut, E.D. Fern´andez-Nieto, A. Mangeney, and P.Y. Lagr´ee. On new erosion models of savage-hutter type for avalanches. Acta Mechanica, 199(1):181–208, 2008. [6] F. Bouchut, E.D. Fern´andez-Nieto, A. Mangeney, and G. Narbona-Reina. A two-phase shallow debris flow model with energy balance. ESAIM: Mathematical Modelling and Numerical Analysis, 49(1):101 – 140, January 2015. [7] F. Bouchut, E.D. Fern´andez-Nieto, A. Mangeney, and G. Narbona-Reina. A two-phase two-layer model for fluidized granular flows with dilatancy effects. Journal of Fluid Mechanics, 801:166–221, 2016. [8] F. Bouchut, A. Mangeney, B. Perthame, and J.P. Vilotte. A new model of saint venant and savagehutter type for gravity driven shallow water flows. Comptes Rendus Mathematique, 336(6):531 – 536, 2003. [9] F. Bouchut and M. Westdickenberg. Gravity driven shallow water models for arbitrary topography. Commun. Math. Sci., 2(3):359–389, 09 2004. [10] Morgane Brunet, Laurent Moretti, Anne Le Friant, Anne Mangeney, Enrique Domingo Fern´andez Nieto, and Francois Bouchut. Numerical simulation of the 30–45 ka debris avalanche flow of montagne pel´ee volcano, martinique: from volcano flank collapse to submarine emplacement. Natural Hazards, 87(2):1189–1222, Jun 2017. [11] M.J. Castro-D´ıaz and E.D. Fern´andez-Nieto. A class of computationally fast first order finite volume solvers: Pvm methods. SIAM Journal on Scientific Computing, 34(4):A2173–A2196, 2012. [12] M.J. Castro-D´ıaz, E.D. Fern´andez-Nieto, J.M. Gonz´alez-Vida, and C. Par´es-Madro˜nal. Numerical treatment of the loss of hyperbolicity of the two-layer shallow-water system. Journal of Scientific Computing, 48(1):16–40, Jul 2011. [13] D. Dutykh and H. Kalisch. Boussinesq modeling of surface waves due to underwater landslides. Nonlinear Processes in Geophysics, 20:267–285, April 2013. 32 pages, 16 Figures, 68 references. Other author’s papers can be downloaded at http://www.denys-dutykh.com/. 38 [14] E.D. Fern´andez-Nieto, F. Bouchut, D. Bresch, M.J. Castro D´ıaz, and A. Mangeney. A new savagehutter type model for submarine avalanches and generated tsunami. Journal of Computational Physics, 227(16):7720 – 7754, 2008. [15] D.L. George, R.M. Iverson, and C.M. Cannon. New methodology for computing tsunami generation by subaerial landslides: Application to the 2015 tyndall glacier landslide, alaska. Geophysical Research Letters, 44(14):7276–7284, 2017. [16] S.T. Grilli, M. Shelby, O. Kimmoun, G. Dupont, D. Nicolsky, G. Ma, J.T. Kirby, and F. Shi. Modeling coastal tsunami hazard from submarine mass failures: effect of slide rheology, experimental validation, and case studies off the us east coast. Natural Hazards, 86(1):353–391, Mar 2017. [17] S.T. Grilli and P. Watts. Tsunami generation by submarine mass failure. i: Modeling, experimental validation, and sensitivity analyses. Journal of Waterway, Port, Coastal, and Ocean Engineering, 131(6):283–297, 2005. [18] P.H. Heinrich, A. Piatanesi, and H. H´ebert. Numerical modelling of tsunami generation and propagation from submarine slumps: the 1998 papua new guinea event. Geophysical Journal International, 145(1):97–111, 2001. [19] J. Horrillo, A. Wood, G.-B Kim, and A. Parambath. A simplified 3-d navier-stokes numerical model for landslide-tsunami: Application to the gulf of mexico. Journal of Geophysical Research: Oceans, 118(12):6934–6950, 2013. [20] K. Hutter and I. Luca. Two-layer debris mixture flows on arbitrary terrain with mass exchange at the base and the interface. Continuum Mechanics and Thermodynamics, 24(4):525–558, Nov 2012. [21] R.M. Iverson. Elementary theory of bed-sediment entrainment by debris flows and avalanches. Journal of Geophysical Research: Earth Surface, 117(F3), 2012. [22] R.M. Iverson and R.P. Denlinger. Flow of variably fluidized granular masses across threedimensional terrain: 1. coulomb mixture theory. Journal of Geophysical Research: Solid Earth, 106(B1):537–552, 2001. [23] R.M. Iverson and D.L. George. A depth-averaged debris-flow model that includes the effects of evolving dilatancy. i. physical basis. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 470(2170), 2014. [24] L. Jiang and P. Leblond. The coupling of a submarine slide and the surface waves which it generates. Journal of Geophysical Resources, 97:731 – 744, 1992. [25] G.S. Khakimzyanov, O. I. Gusev, S.A. Beizel, L.B. Chubarov, and N.Y. Shokina. Simulation of tsunami waves generated by submarine landslides in the black Sea. Russ. J. Numer. Anal. Math. Modelling, 30(4):227–237, 2015. 39 [26] James T. Kirby, Fengyan Shi, Dmitry Nicolsky, and Shubhra Misra. The 27 april 1975 kitimat, british columbia, submarine landslide tsunami: a comparison of modeling approaches. Landslides, 13(6):1421–1434, Dec 2016. [27] I. Luca, K. Hutter, C.Y. Kuo, and Y.C. Tai. Two-layer models for shallow avalanche flows over arbitrary variable topography. International Journal of Advances in Engineering Sciences and Applied Mathematics, 1(2):99–121, Sep 2009. [28] I. Luca, K. Hutter, Y.C. Tai, and C.Y. Kuo. A hierarchy of avalanche models on arbitrary topography. Acta Mechanica, 205(1):121–149, Jun 2009. [29] G. Ma, J.T. Kirby, and F. Shi. Numerical simulation of tsunami waves generated by deformable submarine landslides. Ocean Modelling, 69:146 – 165, 2013. [30] Gangfeng Ma, James T. Kirby, Tian-Jian Hsu, and Fengyan Shi. A two-layer granular landslide model for tsunami wave generation: Theory and computation. Ocean Modelling, 93:40 – 55, 2015. [31] Gangfeng Ma, Fengyan Shi, and James T. Kirby. Shock-capturing non-hydrostatic model for fully dispersive surface wave processes. Ocean Modelling, 43-44:22 – 35, 2012. [32] A. Mangeney, F. Bouchut, N. Thomas, J. P. Vilotte, and M. O. Bristeau. Numerical modeling of self-channeling granular flows and of their levee-channel deposits. Journal of Geophysical Research: Earth Surface, 112(F2), 2007. [33] A. Mangeney-Castelnau, F. Bouchut, J. P. Vilotte, E. Lajeunesse, A. Aubertin, and M. Pirulli. On the use of saint venant equations to simulate the spreading of a granular mass. Journal of Geophysical Research: Solid Earth, 110(B9), 2005. [34] R.W.D. Nickalls. A new bound for polynomials when all the roots are real. The Mathematical Gazette, 95(534):520–526, 2011. [35] S.B. Savage and K. Hutter. The dynamics of avalanches of granular materials from initiation to runout. part i: Analysis. Acta Mechanica, 86(1):201–223, Mar 1991. [36] D. R. Tappin, P. Watts, and S. T. Grilli. The papua new guinea tsunami of 17 july 1998: anatomy of a catastrophic event. Natural Hazards and Earth System Sciences, 8(2):243–266, 2008. [37] S. Viroulet, C. Kharif, O. Kimmoun, and D. C´ebron. Shallow water waves generated by subaerial solid landslides. Geophysical Journal International, 193(2):747–762, 02 2013. [38] P. Watts, S. T. Grilli, J. T. Kirby, G. J. Fryer, and D. R. Tappin. Landslide tsunami case studies using a boussinesq model and a fully nonlinear tsunami generation model. Natural Hazards and Earth System Sciences, 3(5):391–402, 2003. [39] S. Yavari-Ramshe and B. Ataie-Ashtiani. Numerical modeling of subaerial and submarine landslidegenerated tsunami waves—recent advances and future challenges. Landslides, 13(6):1325–1368, Dec 2016. 40 View publication statsView publication stats