A new Savage-Hutter type model for submarine avalanches and generated tsunami

A new Savage-Hutter type model for submarine avalanches and generated tsunami E.D. Fern´andez-Nieto, ∗ , F. Bouchut † , D. Bresch ‡ , M.J. Castro D´ıaz § , & A. Mangeney ¶ Abstract In this paper we present a new two-layer model of Savage-Hutter type to study submarine avalanches. A layer composed of fluidized granular material is assumed to flow within an upper layer composed of an inviscid fluid (e. g. water). The model is derived in a system of local coordinates following a non-erodible bottom and takes into account its curvature. We prove that the model verifies an entropy inequality, preserves water at rest for a sediment layer and their solutions can be seen as particular solutions of incompressible Euler equations under hydrostatic assumptions. Buoyancy effects and the centripetal acceleration of the grain movement due to the curvature of the bottom are considered in the definition of the Coulomb term. We propose a two-step Roe type solver to discretize the presented model. It exactly preserves water at rest and no movement of the sediment layer, when its angle is smaller than the angle of repose, and up to second order all stationary solutions. Finally, some numerical tests are performed by simulating submarine and sub-aerial avalanches as well as the generated tsunami. 1 Introduction Recent improvements in seabed and sub-surface mapping techniques as bathymetry measurements and seismic imagery, have revealed a large amount of slide scars and a wide diversity of related deposits on many of the world’s continental margins [e. g. Locat and Mienert, 2003, Vanneste et al., 2006]. Submarine avalanches or landslides are poorly studied compared to their subaerial counterparts. This is however a key issue in geophysics. Indeed, submarine granular flows driven by gravity participate in the evolution of the sea floor and in particular of the continental margins. They also represent a threat to the submarine infrastructures, especially for the oil or port industry as well as to many sea shore inhabitants due to the potential tsunamis that can be triggered by such landslides. In this paper we present a new two-layer Savage-Hutter type model, with application to sub-aerial/submarine avalanches over variable topography and generated tsunami. The first layer is filled with a homogenous inviscid fluid with constant density and the second layer is made of a fluidized granular mass. The two fluids (i. e. water and fluidized debris) are assumed to be immiscible. Important questions are (i) the rheological behavior of the fluidized granular mass on a complex topography and (ii) the interaction between the two layers. ∗Departamento de Matem´atica Aplicada I, Universidad de Sevilla. E.T.S. Arquitectura. Avda, Reina Mercedes, s/n. 41012 Sevilla, Spain ([email protected]) †D´epartement de math´ematiques et applications, CNRS & ´ Ecole normale sup´erieure, 45, rue d’Ulm, 75230 Paris cedex 05, France ([email protected]) ‡Laboratoire de Math´ematiques, UMR 5127 CNRS, Univ. Savoie, 73376 Le Bourget du Lac (France) ([email protected]) §Departamento de An´alisis Matem´atico, Universidad de M´alaga. F. Matem´aticas, Campus Teatinos S/N, Spain ([email protected]) ¶Equipe de Sismologie, IPGP, 4, pl. Jussieu, 75232 Paris cedes 05, France and Institute for Nonlinear Science, University of California San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0402, USA ([email protected].fr) 1 Numerical modeling of sub-aerial debris or snow avalanches has been extensively investigated during this last decade with application to both laboratory experiments dealing with granular flows and geological events (see e. g. [32], [29], [39], [49], [21], [2], [3], [27], [5], [1]). Most of the models devoted to gravitational granular flows describe the behavior of dry granular material following the pioneer work of Savage and Hutter (see [45]): a shallow-water type model (i. e. thin layer approximation for a continuum medium) is derived to describe granular flows over a slopping plane based on Mohr-Coulomb considerations: a Coulomb friction is assumed to reflect the avalanche/bottom interaction and the normal stress tensor is defined by a constitutive law relating the longitudinal and the normal stresses through a proportionality factor K. New Savage-Hutter models over a general bottom have been proposed by Bouchut et al. in [6], that take into account the curvature of the bottom. The authors introduce two new models: the first one is deduced under the hypothesis of small variation of the curvature and the second one deals with a general bottom topography. The new curvature terms introduced in the models are necessary for two reasons: they make it possible to preserve water at rest solutions and to exactly verify an energy inequality. In this paper we consider the first hypothesis, i. e. a small variation of the curvature. The equations are derived in a local coordinate system attached to the non-erodible topography and takes into account its curvature (see [6]), in particular the centripetal acceleration due to the bottom curvature. A generalization to 2D aerial avalanches over surfaces with small lateral curvature has been carried out in [49] and [41]. In [7], Bouchut and Westdickenberg generalize the previous models for small or for general slope variation in two dimensions. The discretization of 2D aerial avalanches can be done for example by finite volume by using kinetic schemes [31], by Roe type finite volume methods [11], or distribution schemes [43]. A two-layer Shallow Water type model with compressible effects has been introduced in [34] by Morales de Luna. He considers an upper compressible and a lower incompressible layer. The model is presented in local coordinates, verifies an entropy dissipation inequality and gives an approximation of the free surface compressible-incompressible Euler equations. In most industrial applications and real debris flows, the fluid which is present in the granular material cannot be neglected. Recent attempts have been developed to describe mixtures of grains and fluids in shallow-water two-phase or mixture models ([26], [40], [37], [42]). Iverson and Denlinger extend the SH model in [26] to study avalanches of fluidized granular masses where the pores between the grains are assumed to be filled with a fluid. In [42], Pudasaini, Wang and Hutter generalized the work [26] for a general channel on local coordinates. In both works a simplified system is considered, assuming that the velocity of the fluid within the pores is equal to the velocity of the grains. The same hypothesis is used here: the fluidized mass is assumed to be a porous medium composed of sand grains, filled with the fluid present in the upper layer (see [26]). The dissipation within the granular medium is modeled by a Coulomb friction law taking into account the buoyancy effects over the sand grains. The other key point concerns the definition of the stress tensor for the fluid and grain phases of the second layer. From the vertical momentum equation and dimensional analysis, the vertical stress tensor of the complete layer can be derived. However, it is necessary to know the stress tensor for each phase (i. e. fluid and solid phase) in order to apply different constitutive relations for the fluid and solid phase separately. Therefore, additional hypotheses have to be introduced (see [26]). Finally, very few models have been proposed to deal with the interaction of a fluidized mass and the surrounding fluid in which the avalanche propagates. One of the outcome of the interaction between water and debris is the generation of water waves and possible tsunami for particular configurations of the coastal topography and of the submarine avalanche. Most of the models dedicated to the simulation of landslide generated tsunamis reduce the trigger mechanism to a vertical motion imposed as boundary condition in the water wave propagation model (see for example [20]). Submarine landslides are actually modeled by partially or totally submerged pistons, rigid bodies entering the water or initial water displacement (see e. g. [44], [35]). More recently, submarine avalanche dynamics has been taken into account using depth-averaged or full Navier-Stokes models 2 describing the rheological behavior by a Coulomb friction law or by viscous dissipation ([22], [23], [30], [19]). A similar attempt has been performed by Heinrich et al. in [24] but without taking into account the effects of the fluid on the landslide dynamics (i. e. the sea-bottom deformation induced by the landslide is used as input data in the tsunami model). As a result, the momentum equation for the fluidized granular material does not contain any coupling terms between the two layers which should appear in the pressure gradient terms. Other systems, named active models, with a dynamic displacement of sea bed are used with a coupling between a shallow-water system and visco-elastic equations, see for instance [16], [17]. The interested reader is referred to [15] and [18] for references around Tsunamis and challenging modeling. In this paper we present a 1D model for submarine avalanches, which is a generalization of the Savage Hutter (SH) 1D model [45] for aerial avalanches and the model proposed by Heinrich et al. in [24]. To discretize the model that we introduce in the paper, we propose a well-balanced finite volume method. Firstly we begin by rewriting the model obtained in local coordinates to Cartesian coordinates. We can write the model as a hyperbolic system with conservative terms, source terms and non-conservative terms. One of the characteristics of the model is that as we consider the variations of the topography, the physical flux function depends on the variable xmeasured along the horizontal coordinate. This dependence of the flux with respect to xmakes difficult the derivation of an exact well-balanced method for water at rest (see [9], [36]). Moreover, for the proposed model the water at at rest solution should be understood as: no movement of the water column and no movement of the sediment layer when the angle of the sediment surface is smaller than the angle of repose. In such situations, it is necessary to discretize properly the source terms due to the variations of the bottom angle and the derivatives of the flux function with respect to the angle. The more specific difficulty related to discretizisation of system comes from the Coulomb friction term. Its discretization is important, to simulate properly the landslides and to preserve the stationary solutions corresponding to water at rest and no movement of the sediment layer. We propose a two-step numerical scheme to treat the Coulomb friction term. In a first step, a discretization of a term that can be interpreted as a redefinition of the Coulomb term for stationary solutions is considered. This term is only introduced in the uncentered component of the numerical scheme. In the second step, a semi-implicit treatment of the Coulomb term at each cell is performed. We proof that the numerical scheme constructed in this way preserves the solutions corresponding to water at rest and no movement of the sediment layer for angles smaller than the angle of repose. The paper is organized as follows: the model is derived in Section 2. Section 3 is devoted to study the model properties. In Section 4 we present a well-balanced finite volume numerical scheme to discretize the model. We proof that the numerical scheme exactly preserves water at rest and no movement of the sediment layer, and up to second order all stationary solutions. Finally, in Section 5 a series of numerical tests are performed, including the simulation of a tsunami generated by the motion of a sediment layer, following [24]. In Appendix A we present the details related to the change of variable used in the incrompressible Euler equations to write the equation in a local coordinate system attached to the bottom topography. 2 Derivation of the model In this section, we present the derivation of a two-layer model of Savage-Hutter type to study submarine avalanches and generated tsunamis. We denote with index 1 the upper layer, composed of a homogeneous inviscid fluid of constant density ρ1. We also consider a grain layer of density ρs, and porosity ψ0(see Figure 1). We consider that the pores in the grain layer are filled with the 3                           θ b x z X=(x,z) Z X η0 η ηS h2 1 h 2 h Water surface ρ ρ1 2 Figure 1: A fluid layer over a grain layer and a non-erodible bottom b fluid of the upper layer. Then, the density of layer 2 composed of the fluidized mass is defined as ρ2= (1 −ψ0)ρs+ψ0ρ1.(1) First, the system of equations describing the dynamics of the two-layer system is presented. Next, a change of variables to local coordinates attached to the bottom (see Appendix) is performed and the boundary and kinematic conditions are set. The final model is derived based on a dimensional analysis and a vertical integration of the equations. Starting system of equations We consider the incompressible Euler equations. The unknowns are ~ Vi=ui vi, i = 1,2, being uiand vi, the horizontal and vertical velocity components of each layer, respectively. Then, the incompressible Euler equations can be written as div~ Vi= 0, i = 1,2,(2) ρi∂t~ Vi+ρi~ Vi∇~ Vi=−divPi+ρi∇(~g ·~ X), i = 1,2,(3) where we denote by Pi,i= 1,2, the pressure tensor of each layer Pi=pi,x x pi,x z pi,z x pi,z z , i = 1,2, (with pi,x z =pi,z x), by ρi,i= 1,2 the densities of each layer, by ~ Xa point in Cartesian coordinates ~ X= (x, z), and ~g = (0,−g). In order to model the evolution of the granular layer using the Euler equations, following [26] we suppose that the velocity of the fluid in the pores of the second layer and the grains are the same and P2can be decomposed as P2=Ps 2+Pf 2, where Ps 2and Pf 2are the pressure tensor of the solid phase (grains) and the fluid phase, 1respectively. 1In a binary mixture model the pressure tensor of the mixture is exactly given by P=Ps+Pf− 2 X α=1 ρα(~ Vα−~ Vb)⊗(~ Vα−~ Vb),where ~ Vb= 2 X α=1 ρα~ Vα P2 α=1 ρα is the barycentric velocity. This reduces to P=Ps+Pfif the fluid and solid velocities are the same. 4 Next, a change of variables is performed to equations (2)-(3). Local variables over a nonerodible bottom defined by z=b(x) are considered. Xdenotes the arc’s length of the bottom and Zis measured orthogonally to the bottom (see Figure 1). In what follows we denote by h1and h2the thickness of the fluid and grain layers, respectively, measured orthogonally to the bottom (see Figure 1), by S=h1+h2the free water surface. The details of this change of variables is given in the Appendix. The change of variables is valid when the local radius of curvature of the bed is smaller than h1+h2. Equations (2)-(3) are re-written in the new variables as                    ∂X(Ui) + ∂Z(J Wi) = 0, i = 1,2, ρi∂t(J Ui) + ρi∂X(U2 i) + ρi∂Z(JWiUi) + ρi∂X(~g ·~ X) = −∂X(Pi XX )−∂Z(JPi ZX )+ +ρiWi(∂X(Uiθ) + ∂Z(JWiθ)) + Pi XZdXθ. i = 1,2, ρi∂t(J Wi) + ρi∂X(UiWi) + ρi∂Z(JW2 i) + ρiJ∂Z(~g ~ X) = −∂X(Pi XZ)− −∂Z(JPi ZZ )−ρiUi(∂X(Uiθ) + ∂Z(JWiθ)) − Pi XX dXθ, i = 1,2, (4) where, we denote by Ui,i= 1,2, the velocity parallel to the bottom and by Wi,i= 1,2, the velocity perpendicular to the bottom, with ireferring to layers 1 and 2. The pressure tensor Piis defined by Pi=cosθsinθ −sinθcosθPicosθ−sinθ sinθcosθ=Pi,XX Pi,XZ Pi,ZX Pi,ZZ . Observe that as pi,xz =pi,xz then Pi,XZ =Pi,ZX . Moreover, let us recall that ρ1is the density of the fluid and that ρ2is defined by (1). θis the angle between the tangent vector of the bottom and the horizontal (see Figure 1), and J= 1−ZdXθ is the Jacobian of the change of variables (note, dXθ=∂Xθ, for a non-erodible bed, see Appendix). Observe that J6= 0 if the local radius of curvature of the bed is smaller than h1+h2. Boundary and kinematic conditions We denote by ηSthe unitary normal vector to the free water surface Z=S(S=h1+h2) with positive vertical component, by ηh2the unitary normal vector to the surface Z=h2and by η0= (0,1) the corresponding unitary normal vector to the bottom (Z= 0). The following kinematic conditions are considered ∂tS+U1∂XS−W1= 0,(5) ∂th2+Ui∂Xh2−Wi= 0, i = 1,2.(6) Finally, the following boundary conditions are imposed: •On Z=S: P1·ηS= 0.(7) •On Z=h2: ηh2·(P1− P2)ηh2= 0 (8) Pi·ηh2−ηh2(ηh2· Piηh2) = fric(U1, U2) 0i= 1,2,(9) where fric(U1, U2) is a friction term between both layers. 5 •On Z= 0: (U, W)·η0= 0 ⇒W= 0,(10) P2η0−η0(η0· P2η0) =  −η0·(P2− P1)η0U0 2 |U0 2|tan(δ0) 0 .(11) Note that the Coulomb friction law in equation (11) takes into account the buoyancy effects due to the fact that the grains are submerged within a fluid layer. Remark 1 Equation (9) assumes no water exchange between the two layers. Nevertheless there is a water exchange between the fluid an the porous avalanche, so equation (9) constitutes a simplification of the problem. This entrainment process has been studied first by Beaves and Joseph in [4]. Equation (6) assumes that the second layer has constant porosity (volume fraction) since ρ2is constant then ψ0=constant. Dimensional analysis Next, a dimensional analysis of the set of equations (4), the kinematic and boundary conditions is performed. The non-dimensional variables ( e.) read: (X, Z, t) = (Le X, H e Z, (L/g)1/2et), (Ui, Wi) = (Lg)1/2(f Ui, εf Wi), i = 1,2, hi=He hi, i = 1,2,(Pi XX ,Pi ZZ) = gH(e Pi XX ,e Pi ZZ ), i = 1,2, Pi XZ =gHµie Pi XZ , i = 1,2, (12) where µ1= 1, µ2= tan(δ0), δ0being the angle of repose in the Coulomb term (see [45]). By L and Hwe denote, respectively, the characteristic lengths tangential and normal to a representative basal direction of the domain. We suppose a shallow domain, so ε=H/L is supposed to be small. Note that the Savage Hutter model has been shown to reproduce experimental granular collapse over horizontal plane for aspect ratio ≤0.5 [32]. Using the above change of variables, the system of equations (4) are re-written as (we omit the tildes): ∂X(Ui) + ∂Z(J Wi) = 0, i = 1,2,(13) J∂t(ρiUi) + ρiUi∂XUi+ρiJWi∂ZUi+ρi∂X(b+Zcosθ+Pi XX ρi )ε= =−µi∂Z(JPi XZ ) + ρiWiεUidXθ+∂XθPi XZ µiε, i = 1,2,(14) ε{J∂t(ρiWi) + ρiUi∂X(Wi) + ρiWi∂Z(Wi) + ∂X(Pi XZ )−∂XθPi XX − Pi ZZ dXθ}+ +ρiJ∂Z(b+ cosθZ) = −J∂Z(Pi ZZ)−ρiU2 idXθ, i = 1,2.(15) The kinematic conditions (5)-(6) are re-written as: ∂tS+U1∂XS−W1= 0, ∂th2+Ui∂Xh2−Wi= 0, i = 1,2.(16) Finally, the boundary conditions (7)-(11) are now given as: 6 •Z=S: On Z=Swe have, ηS= (−ε∂XS, 1)/ϕSwith ϕS=p1 + ε2(∂XS)2, then from (7) we obtain −ε∂XSP1XX +µ1P1ZX = 0,(17) −ε∂XSµ1P1XZ +P1ZZ = 0.(18) •Z=h2: On Z=h2we have, ηh2= (−ε∂Xh2,1)/ϕh2with ϕh2=p1 + ε2(∂Xh2)2, then from (8) and (9) we obtain P1ZZ =P2ZZ +O(ε),(19) −εPi XX ∂Xh2+µiPi XZ =−(ηh2Piηh2)(ε∂Xh2) + fric(U1, U2), i = 1,2,(20) −εµiPi ZX ∂Xh2+Pi ZZ = (ηh2Piηh2)i= 1,2.(21) •Z= 0: On Z= 0, we have η0= (0,1), then from (10) and (11) we obtain W2= 0,(22) µ2P2XZ =−(P2ZZ − P1ZZ )U0 2 |U0 2|tan(δ0).(23) Constitutive laws. We suppose that dXθ=O(ε). Then from (15) we obtain ∂Z(P1ZZ ) = −ρ1cosθ+O(ε),(24) ∂Z(P2ZZ ) = −ρ2cosθ+O(ε).(25) If we integrate (24) from Z > 0 to S, we have, up to order ε, P1ZZ =ρ1(S−Z)cosθ, (26) therefore, P1ZZ (h2) = ρ1h1cosθ. Using this last expression, the relations given in (19) and integrating (25) from Z > 0 to S, we have, up to first order Ps 2ZZ +Pf 2ZZ =P2ZZ =ρ1h1cosθ+ρ2cosθ(h2−Z).(27) The last equation defines the total pressure, P2ZZ , perpendicular to the base. The constitutive relation for both the grains and the fluid, i. e. Ps 2ZZ and Pf 2ZZ , are required to close the model. The same problem appears if we study a grain-fluid mixture aerial avalanche. See for example [26] and [41]. In order to obtain an expression for the normal stress of both phases, they suppose that both are linear in Z. Moreover, they suppose that the component of the stress tensor of the fluid phase normal to the basal surface is proportional to the pressure of a fluid layer, without the solid phase. We adapt this hypothesis to our case, taking account of the fact that the fluidized layer has an upper layer of fluid. Concretely, we suppose Pf 2ZZ (Z) = λ1ρ1h1cosθ+λ2ρ1h2cosθ(h2−Z),(28) where λ1and λ2are two parameters. Moreover, by (27), we have Ps 2ZZ (Z) = ρ1h1cosθ(1 −λ1) + cosθ(h2−Z)(ρ2−λ2ρ1) (29) The study of the stress transition conditions at a singular surface between two mixtures which do not have the same number of components, is a very difficult subject. Hypothesis (28) can be seen as a fist trial, in the context of this paper. Some earlier papers looking at the problem of interfacial transition conditions with different number of constituents have been included in the references. For example in [25] Hutter et al. study the transition conditions, with application to glaciers where the upper layer is ice and the under layer is a sediment-ice mixture (see also [47], [50], [51]). 7 Remark 2 Comparisons with experiments are necessary to define λ1and λ2. Nevertheless we can make some possible choices. The first simplification is to consider λ1=λ2. In this case the component of the stress tensor of the fluid phase normal to the base is proportional to the pressure of a fluid, without the solid phase in the second layer. Nevertheless, we prefer at this moment to retain two different parameters. Because the role of λ1could be different of that of λ2. Observe that if we evaluate (28) and (29) for Z=h2we obtain Pf 2ZZ (h2) = λ1ρ1h1cosθ, Ps 2ZZ (h2) = ρ1h1cosθ(1 −λ1); Note that λ1controls the distribution of the pressure at the interface into the two phases of the second layer. We have imposed continuity of the component of the stress tensor normal to the base across the interface of the first and second layer (equation (8)), and we observe that it is verified independently of the definition of λ1. If we want to include an additional condition, for example the continuity of the pressure of the fluid phase of the second layer with the first layer of fluid, then we obtain that λ1= 1. Depending on the material of the second layer, we can also suppose that the fluid that fills the pores of the second layer is nearly isolated of the fluid of the first layer, in this case we can consider λ1≈0. Independently of the additional hypothesis that we can use to set the distribution of the pressure at the interface between the solid and fluid phase, we have still the parameter λ2at our disposal, in order to impose a similar hypothesis to that introduced by Iverson and Delinger, but only for the second layer. Another possible choise is to fix λ1=λ2=ψ0, where ψ0is the porosity of the second layer. We obtain in this case Ps 2ZZ = (1 −ψ0)(ρ1h1+ρs(h2−Z))cosθ, Pf 2ZZ =ψ0ρ1(h1+h2−Z)cosθ. An interesting property of this choice is that P2ZZ at height Z=h2is proportional to the pressure that is obtained in absence of the fluid phase (with proportional constant (1 −ψ0)). Observe that Ps 2ZZ depends on ρs, the density of the solid phase, and not ρ2, the density of the mixture defined by (1). Analogously, Pf 2ZZ is proportional to the pressure (with proportional constant ψ0) that is obtained in the absence of the solid phase. Moreover, if the porosity of the second layer is zero, then P2=Ps 2, which is automatically deduced from this definition. Finally, the following relations are also considered (see for example [26] and [41]): P1XX =P1ZZ ,Ps 2XX =KPs 2ZZ ,Pf 2XX =Pf 2ZZ , where Kmeasures the anisotropy or normal stress effects in the solid phase. The definition of K can be done in different ways. For example Heinrich et al. in [24] consider K= 1, other definitions of Kcan be found in [26]. The effects related to the definition of Kin numerical modelling of experimental and natural flows is studied in [41] and [38]. Remark 3 The value K= 1 corresponds to isotropic conditions, K6= 1 makes ‘overburden pressures’ different from the normal stresses parallel to the basal surface. In soil mechanics, Kcorresponds to the earth pressure coefficient, see [41]. For non-Neutonian rheology Kmay also be different from unity. Using the previous relations, the following expression for P2XX is derived: P2XX =KPs 2ZZ +Pf 2ZZ = =h1cosθρ1(λ1+K(1 −λ1)) + (h2−Z)cosθ(λ2ρ1+K(ρ2−λ2ρ1)).(30) 8 Now, replacing (26) and (30) in (14), and using the incompressibility equation (13), we obtain up to second order ∂t(ρ1U1) + ρ1∂XU2 1+ρ1∂Z(U1W1) + ρ1∂X(b+Scosθ)ε=−µ1∂Z(P1XZ ),(31) and ∂t(ρ2U2) + ρ2∂XU2 2+ρ2∂Z(U2W2) + ρ2∂Xb+Zcosθ+1 ρ2 [h1cosθρ1(λ1+K(1 −λ1))+ +(h2−Z)cosθ(λ2ρ1+K(ρ2−λ2ρ1))]ε=−µ2∂Z(P2XZ ).(32) Integration process In this section, equations (31), (32) and (13) are depth-averaged in the direction normal to the topography. Let us introduce the following notation: we denote by ¯ Ui,i= 1,2 the velocities of each layer averaged perpendicular to the basal surface: ¯ U1=1 h1ZS h2 U1(X, Z)dZ, ¯ U2=1 h2Zh2 0 U2(X, Z)dZ. We also denote U2 1=1 h1ZS h2 U2 1(X, Z)dZ, U2 2=1 h2Zh2 0 U2 2(X, Z)dZ. Assuming that dXθ=O(ε), then J= 1 −ZdXθis reduced to J= 1 up to second order. Therefore, (13) reduces to ∂X(U1) + ∂Z(W1) = 0,(33) ∂X(U2) + ∂Z(W2) = 0.(34) I.1) If equation (33) is integrated from Z=h2to Z=S, we obtain 0 = ∂X(h1¯ U1)−U1(S)∂XS+W1(S) + U1(h2)∂Xh2−W1(h2). Now, using the kinematic conditions (16), the following equation is derived: ∂th1+∂X(h1¯ U1) = 0. I.2) Analogously, by integrating (34) between Z= 0 and Z=h2we obtain 0 = ∂X(h2¯ U2)−U2(h2)∂Xh2+W2(h2)−W2(0), and, using the kinematic condition (16) and the boundary condition (22), the following equation is derived: ∂th2+∂X(h2¯ U2) = 0. I.3) Let us now proceed with equation (31), integrating it from Z=h2to Z=S. We obtain ρ1∂t(h1¯ U1) + ρ1∂X(h1U2 1)−ρ1U1(S)[∂t(S) + U1(S)∂XS−W1(S)]+ +ρ1U1(h2)[∂th2+U1(h2)∂Xh2−W1(h2)] + ρ1ZS h2∂X(b+Scosθ)dZε= =−µ1(P1XZ (S)− P1XZ (h2)).(35) The expressions of P1XZ (S) and P1XZ (h2) are now derived using the boundary conditions and the constitutive laws. 9 Remark 7 The inequality (51) is just the energy conservation for smooth solutions, supposing Kin = 0,δ0= 0 and (1 −Λ2)∂xθ= 0 (that is, Λ2= 1 or θconstant). Remark 8 Observe that the inequality (54) is independent of λ1,λ2and ∂xθwhen K= 1. For K= 1 (54) reduces to |∂x(b+h2cosθ)| ≤ tan(δ0). As the equation of the interface between the fluid and the sediment material is defined by b+h2cosθ, the previous condition implies that the slope of the interface is smaller than δ0. In the case K6= 1, (54) relates the values of Kto the curvature of the bottom, the parameters λ1,λ2and the ratio between the densities of the fluid and the granular material, r(see 2). Observe that for stationary solutions verifying U2= 0 then ∂XU2= 0. Then, if we consider for this case K= (Kact +Kpas)/2, we have K= 1 + tan2φ. (60) By another way, we can study the profiles verifying the equality in (54). We consider a domain [0, L]and we impose the value of the interface at x=L, then for δ0,K,λ1,λ2and rfixed we have the equation αcos θ∂xh2+ ((α+β 2)∂xcos θ)h2= (1 −r) tan(δ0)−(α+β)∂xb (h2cos θ+b)|x=L=A with α= Λ2−rΛ1,β= 1 −Λ2. The solution is h2(x) = I(x)−I(L) + cos θβ/(2α)(A−b(L))cos θ−1+β/(2α) where I(x) = Zx 0(1 −r) tan(δ0) αcos θβ/(2α)−(1 + β α) sin θcos θ−1−β/(2α)dx. In Figure 2 we present two examples of the profiles that we obtain for two different bottom topographies for different values of K. From K= 1 to K= 2 they correspond to the definition (60) with φfrom 0to 45 degrees. In both examples we have set λ1=λ2=ψ0,ψ0= 0.2,δ0= 28 degrees and r= 0.2. Figure 2(a) corresponds to a bottom with constant slope equal to 15 degrees, where L= 2 meters and the interface at x=Lis A= 2. In this example we observe that the interfaces obtained for K > 1are over the interface corresponding to K= 1. The bottom of Figure 2(b) is defined by b(x) = −ln(cos(x)), moreover L= 1.5and A= 3.5. In this example we observe that by the influence of the curvature of the bottom the interfaces corresoponding to K > 1are under the interface obtained for K= 1. 4 Numerical scheme: rewriting the model In this section we describe the numerical scheme that we propose to discretize model (45). We propose a well-balanced finite volume method that exactly preserves the solutions corresponding to water at rest and no movement of the sediment layer verifying (52), (53), (54); and up to second order all stationary solutions. In Subsection 4.1 we introduce the numerical scheme, and we study its properties. However, before defining the numerical scheme we begin by rewriting the proposed model in Cartesian coordinates. We remark that model (45) is written in local coordinates over a non-erodible bottom. In order to solve the problem of defining a proper mesh for an arbitrary topography, we propose to rewrite 16 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Bottom K=1 K=1.5 and K=2 (a) Bottom with constant slope 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 K=1.01 K=1 K=1.1 K=1.2 K=1.5 and K=2 Bottom (b) Bottom with curvature Figure 2: Stationary interface profiles depending on the values of K the model (45) in Cartesian coordinates. To do this, the following rule is used: for a given function f(X(x)), as ∂x = cosθ ∂X ⇒∂f ∂X = cosθ∂f ∂x .(61) Introducing the notation Hi=hi cosθ, Qi=Hi¯ Ui, i = 1,2, equations (45) can be written as                              ∂tH1+∂x(Q1cosθ) = 0, ∂t(Q1) + ∂x(H1¯ U2 1cosθ+gH2 1 2cos3θ) = −g H1cosθ dxb+ +gH2 1 2sinθcos2θ dxθ−g H1cosθ ∂x(H2cos2θ) + fric(U1, U2) cosθρ1 , ∂tH2+∂x(Q2cosθ) = 0, ∂t(Q2) + ∂xH2¯ U2 2cosθ+gΛ2 H2 2 2cos3θ=−g H2cosθ dxb+ +gH2 2 2sinθcos2θ∂xθ−rΛ1g H2cosθ∂x(H1cos2θ)−fric(U1, U2) cosθρ2 +T cosθ, (62) where Tis defined by If |T | ≥ σc⇒ T =−(g(1 −r)H2cos2θ+H2cosθ¯ U2 2dx(sinθ)) Q2 |Q2|tan δ0,(63) If |T | < σc⇒Q2= 0,(64) where σc=g(1 −r)H2cos2θtan(δ0). 4.1 Well-balanced finite volume method In this subsection we present the finite volume method that we use to discretize system (62). There are several difficulties related to the discretization of this system: As we describe below, we can rewrite (62) under the structure of a hyperbolic system with a conservative term, a non-conservative product and two types of source terms (see equation (65) below), where 17 i) the flux function does not only depend on the vector of unknowns, but also on θ(x); ii) the coupling term is a non-conservative product B(W)∂xW. In general, it is not well defined nor as a distribution and the choice of a family of paths is necessary (see [14]); iii) the source terms G1and G2, are defined as functions of the fixed topography. Their numerical discretization can be treated by following the ideas given in [12] or [13] in the framework of a system of balance laws or by rewriting the system for an extended variable and an extra equation in such a way that the source terms are written in the form of non-conservative products (see [36]); iv) the source term corresponding to the Coulomb term presents a different difficulty. We propose a two step method combining a well-balanced discretization of the Coulomb term and the numerical treatment introduced by Mangeney et al. in [31]. The numerical method constructed in this way is exactly well balanced for the solutions corresponding to water at rest and no movement of the sediment layer given by (52), (53) and (54). We can rewrite model (62) under the form of a hyperbolic system with a conservative product, a non-conservative term and source terms: ∂tW+∂xF(θ, W) = G1(x, W) + G2(x, W) + B(W)∂xW+T, (65) where W=    H1 Q1 H2 Q2    , F(θ, W) =        Q1cosθ Q2 1 H1 cosθ+gH2 1 2cos3θ Q2cosθ Q2 2 H2 cosθ+gΛ2 H2 2 2cos3θ        , G1=    0 −g H1cosθ dxb 0 −g H2cosθ dxb    , G2=      0 −gH1 2(H1 2+ 2H2) cosθ ∂x(cos2θ) 0 −gH2 2(H2 2+ 2 rΛ1H1) cosθ ∂x(cos2θ)      , B(W) =     0 0 0 0 0 0 −gH1cos3θ0 0 0 0 0 −rΛ1gH2cos3θ0 0 0    , T =    0 0 0 T/cosθ    . Note that the terms ∂x(H2cos2θ) and ∂x(H1cos2θ) of the second and fourth equations on (62) contributes to (65) in the definition of the non-conservative term B(W)∂xWand in the definition of G2. Remark 9 Observe also that G2can be written in terms of ∂x(cos3θ), nevertheless we propose to define G2in terms of cosθ∂x(cos2θ), motivated by the discretization that we proposed. The purpose is that we want to obtain an exactly well-balanced numerical scheme for water at rest: ¯ U1=¯ U2= 0, b +H2cos2θ=cst, H1cos2θ=cst, that is defined in terms of cos2θ. Remark 10 System (65) can be written in nonconservative form, ∂tW+A(θ, W)∂xW=S(θ, W) 18 where S(θ, W) = G1+G2−∂θF+T. And where A(θ, W)defines the transport matrix of the system: A(θ, W) =     0cosθ0 0 −U2 1cosθ+gH1cos3θ2U1cosθ gH1cos3θ0 00 0cosθ rΛ1gH2cos3θ0−U2 2cosθ+ Λ2gH2cos3θ2U2cosθ    ,(66) where Ui=Qi/Hirepresents the averaged velocity of the i-th layer, and r=ρ1 ρ2. This matrix is similar to the one obtained in the well-known two-layer Shallow Water system (see [36] for example). Unfortunally no explicit expressions of the eigenvectors of the system can be obtained. The characeristic equation of the system is: λ2−2U1λ+U2 1−gH1cos2θλ2−2U2λ+U2 2−gΛ2H2cos2θ=rΛ1g2H1H2cos4θ. (67) It is not easy to verify the genuinely nonlinear character of the 4 characteristic fields, as the eigenvalues and eigenvectors can not be written explicitely in a simple manner. Nevertheless, this fact is easily proved in the case r= 0 as, in this case, the system reduces to a decoupled system of Shallow Water and Savage-Hutter equations. In this case the eigenvalues are those corresponding to each layer separately. Then, a continuity argument ensures the genuinely nonlinear character of the 4 characteristic fields when ris close to zero. In the case r≈1, in [46] authors gives an approximation of the eiganvalues for the two-layer Shallow Water equations. The case r≈1is the situation arising for two fluid with different densities in many oceanographical flows. In the context of submarine avalanches we suppose two different materials, then we are closer to the case r≈0. Although, following [46], we can also give a approximation of the eiganvalues for r≈1. We obtain: λ± ext ≈U1H1+U2H2 H1+ Λ2H2 ±g(H1+ Λ2H2)1 2cosθ, (68) λ± int ≈U1Λ2H2+U2H1 H1+ Λ2H2 ± g(1 −rΛ1 Λ2 )H1H2Λ2cos2θ (H1+ Λ2H2)h1−(U1−U2)2 g(1 −rΛ1 Λ2)(H1+ Λ2H2)cos2θi!1 2 . (69) From equation (69) we can observe that the internal eigenvalues may become complex. This situation occurs when they verify, approximately, the following inequality: (U1−U2)2 g(1 −rΛ1 Λ2)(H1+ Λ2H2)cos2θ>1.(70) In this case, the system loses its hyperbolic character. These situations are related with the appearance of shear instabilities that may lead, in real flows, to intense mixing of the two layers. While, in practice, this mixture partially dissipates the energy, in numerical experiments these interface disturbances grow and overwhelm the solution. Clearly, we cannot expect to simulate these phenomena with a two-immiscible-layer model. Therefore, the inequality (70) in fact gives the range of validity of a model based on the equations (65). In this work only the case where the matrix A(θ, W)has real eigenvalues is considered, i.e. the system is supposed to be strictly hyperbolic. As the source terms modeling the friction between the two layers are discretized semi-implicitly, they do not appear in the finite volume discretization, therefore, they are supposed to be zero in this section. For the discretization of the system, computing cells Ii= [xi−1/2, xi+1/2] are considered. For simplicity, we suppose that these cells have constant size ∆x. Let us define xi+1 2=i∆xand by xi= (i−1/2)∆x, the center of the cell Ii. Let ∆tbe the constant time step and define tn=n∆t. 19 We denote by Wn ithe approximation of the cell averages of the exact solution provided by the numerical scheme: Wn i∼ =1 ∆xZxi+1/2 xi−1/2 W(x, tn)dx. (71) The source terms G1and G2are discretized following the ideas introduced in [13] and [36]. The discretization of B(W)∂xWfirstly requires to interpret this term as a Borel measure (see [14]), depending on the choice of a family of paths linking given states. Here the family of segments are considered as in [36]. The dependence of the flux function on θ(x), makes it difficult to obtain the desired exact well-balanced property for water at rest (see [9]). Following the same ideas that have been exposed in Remark 9, we propose to consider the flux function F(θ, W ) as a function of cosθand cos2θ. More precisely, F(θ, W) = F(cosθ, cos2θ, W),with F(α, β, W) =            Q1α Q2 1 H1 α+gH2 1 2α β Q2α Q2 2 H2 α+gΛ2 H2 2 2α β            . Finally, as mentioned before, the discretization of the source term T(W) corresponding to the Coulomb friction term is critical to simulate properly the landslides and to preserve the stationary solutions corresponding to water at rest and no movement of the sediment layer verifying (52), (53) and (54). We propose a two-step numerical scheme to treat the Coulomb friction term. Let us suppose that the values Wn iare known. In order to advance in time we proceed as follows: •First Step. We define W∗ i= [H∗ 1,i Q∗ 1,i H∗ 2,i Q∗ 2,i]Tas W∗ i=Wn i−∆t ∆xDFn,+ i−1/2+DFn,− i+1/2,(72) where DFn,± i+1/2=DF± i+1/2(Wn i, Wn i+1) are the generalized Roe flux difference computed using the family of segments. DF± i+1/2(Wi, Wi+1) = 1 2F(cosθi+1/2,(cos2θ)i+1/2, Wi+1) − F(cosθi+1/2,(cos2θ)i+1/2, Wi) + S3,i+1/2(cos2θi+1 −cos2θi) +S4,i+1/2(cosθi+1 −cosθi)−S1,i+1/2(bi+1 −bi) −S2,i+1/2(cos2θi+1 −cos2θi)−Bi+1/2(Wi+1 −Wi)(73) ±Pi+1/2Ai+1/2(Wi+1 −Wi)−S1,i+1/2(bi+1 −bi) +(S3,i+1/2−S2,i+1/2)(cos2θi+1 −cos2θi) + +S4,i+1/2(cosθi+1 −cosθi)−Ti+1/2∆x. The matrices appearing in the definition of the numerical scheme can be written as follows: 20 Ai+1/2=          J1 i+1/2−B1,2 i+1/2 −B2,1 i+1/2J2 i+1/2          , Bi+1/2=          0B1,2 i+1/2 B2,1 i+1/20          ,(74) where J1 i+1/2=0 cosθi+1/2 −¯ U2 1,i+1/2cosθi+1/2+c2 1,i+1/2(cos2θ)i+1/22¯ U1,i+1/2cosθi+1/2, J2 i+1/2=0 cosθi+1/2 −¯ U2 2,i+1/2cosθi+1/2+ Λ2c2 2,i+1/2(cos2θ)i+1/22¯ U2,i+1/2cosθi+1/2, B1,2 i+1/2=0 0 −c2 1,i+1/2(cos2θ)i+1/20, B2,1 i+1/2=0 0 −rΛ1c2 2,i+1/2(cos2θ)i+1/20, S1,i+1/2=    0 −gH1,i+1/2cosθi+1/2 0 −gH2,i+1/2cosθi+1/2    ,(75) S2,i+1/2=       0 −gH1,i+1/2 2(H1,i+1/2 2+ 2H2,i+1/2)cosθi+1/2 0 −gH2,i+1/2 2(H2,i+1/2 2+ 2 rΛ1H1,i+1/2)cosθi+1/2        ,(76) S3,i+1/2=      0 3g 4H2 1,i+1/2cosθi+1/2 0 3g 4Λ2H2 2,i+1/2cosθi+1/2      , S4,i+1/2=    Q1,i+1/2 Q2 1,i+1/2/H1,i+1/2 Q2,i+1/2 Q2 2,i+1/2/H2,i+1/2    ,(77) Ti+1/2=    0 0 0 Ti+1/2/cosθi+1/2    ,(78) where Ti+1/2=   T1,i+1/2+T2,i+1/2if |Q2,i+1/2|>∆tσc,i+1/2 cosθi+1/2 τcrit,i+1/2otherwise (79) with T1,i+1/2=−c2 2,i+1/2cosθi+1/2(1 −r) SGN( ¯ U2,i+1/2) tan(δ0), T2,i+1/2=−H2,i+1/2¯ U2 2,i+1/2 sinθi+1 −sinθi ∆xSGN( ¯ U2,i+1/2) tan(δ0), σc,i+1/2= (1 −r)c2 2,i+1/2cosθi+1/2tan(δ0), 21 τcrit,i+1/2=c2 2,i+1/2cosθi+1/2(Λ2−rΛ1)bi+1 −bi+H2,i+1cos2θi+1 −H2,icos2θi ∆x+ +(1 −Λ2)(bi+1 −bi ∆x+H2,i+1/2 4 cos2θi+1 −cos2θi ∆x).(80) In (73)-(80), we use the definitions ck,1+1/2=qgHk,i+1/2cosθi+1/2,(81) ¯ Uk,i+1/2=pHk,i ¯ Uk,i +pHk,i+1 ¯ Uk,i+1 pHk,i +pHk,i+1 ,(82) and Hk,i+1/2=Hk,i +Hk,i+1 2, k = 1,2,cosθi+1/2=cosθi+ cosθi+1 2, (cos2θ)i+1/2=cos2θi+ cos2θi+1 2, as well as the upwinded matrices Pi+1/2=Ki+1/2(SGN(Di+1/2))K−1 i+1/2.(83) Here, if Di+1/2is a diagonal matrix defined by the eigenvalues of the matrix Ai+1/2,Ki+1/2 is the matrix whose columns are the associated eigenvectors. Let us denote by λj,i+1/2, j= 1, .., 4, the eigenvalues of matrix Ai+1/2, then SGN(Di+1/2) =     sgn(λ1,i+1/2) sgn(λ2,i+1/2) sgn(λ3,i+1/2) sgn(λ4,i+1/2)     •Second step. We define Wn+1 i= [H∗ 1,i Q∗ 1,i H∗ 2,i Qn+1 2,i ]T and Qn+1 2,i =       Q∗ 2,i + (T∗ 1,i +T∗ 2,i) ∆tif |Q∗ 2,i|>σ∗ c,i∆t cosθi , 0 otherwhise, (84) with T∗ 1,i =−(1 −r)(c∗ 2,i−1/2)2+ (c∗ 2,i+1/2)2 2cosθiSGN(Q∗ 2,i) tan(δ0), T∗ 2,i =−H∗ 2,i−1/2+H∗ 2,i+1/2 2(¯ U∗ 2,i)2sinθi+1/2−sinθi−1/2 ∆xSGN(Q∗ 2,i) tan(δ0), in which σ∗ c,i = (1 −r)(c∗ 2,i−1/2)2+ (c∗ 2,i+1/2)2 2cosθitan(δ0).(85) with c∗ 2,i+1/2=rgH∗ 2,i +H∗ 2,i+1 2cosθi+1/2. 22 The definition of Qn+1 2,i proposed in equation (84), is based on the numerical treatment of Coulomb friction term introduced by Mangeney et al. in [31]. Observe that the definition of the Coulomb term, implies that if |T|< σcthen Q2= 0. A way to impose implicitely this definition in the numerical scheme is the one proposed by equation (84). Remark 11 Observe that S3,i+1/2(cos2θi+1 −cos2θi) + S4,i+1/2(cosθi+1 −cosθi) ∆x is a second order approximation of ∂θF dxθ|x=xi+1/2. This numerical scheme could be seen as a predictor-corrector numerical scheme for the Coulomb friction term. In the first step, the term Ti+1/2is only considered in the uncentered part of the numerical scheme. Note that in the definition of Ti+1/2(see (79)) a second order approximation of the Coulomb friction term is considered if |Q2,i+1/2|>∆tσc,i+1/2 cosθi+1/2 . Otherwise, we set Ti+1/2= τcrit,i+1/2, that is also a second order approximation of the value of the Coulomb friction term in order that all terms in the last equation of system (65) are balanced taking into account U2= 0. This relation is critical in order to obtain a well-balanced numerical scheme for the solutions corresponding to water at rest and no movement of the sediment in the model. After this first step, a predicted value Q∗ 2,i is computed and then, following [31], the final value Qn+1 2,i is computed using (84). Concerning the stability requirements, we use the following CFL-condition max{kDi+1/2k∞,0≤i≤M}∆t ∆x≤γ, where 0 < γ ≤1, and Mis the number of cells into which the space domain is decomposed. We have the following result: Theorem 2The previous numerical scheme verifies the following properties: i) The numerical scheme preserves all the stationary solutions satisfying ¯ U1=¯ U2= 0, b+ (H1+H2)cos2θ=cst, and |(Λ2−rΛ1)∂x(b+H2cos2θ) + (1 −Λ2)(∂xb+H2 4∂xcos2θ)| ≤ (1 −r) tan(δ0), such as (Λ2−rΛ1)b(xi+1)−b(xi) + H2(xi+1) cos2(θ(xi+1)) −H2(xi) cos(θ(xi))+ +(1 −Λ2)b(xi+1)−b(xi) + H2(xi+1) + H2(xi) 4(cos2(θ(xi+1)) −cos2(θ(xi)))≤ ≤(1 −r) tan(δ0)∆x. (86) ii) The numerical scheme preserves all stationary solutions up to order 2. 23 PROOF.- Using the definition of the Roe matrix Ai+1/2it is easy to prove that Ai+1/2(Wi+1 −Wi) = F(cosθi+1/2,(cos2θ)i+1/2, Wi+1)− F(cosθi+1/2,(cos2θ)i+1/2, Wi)− −Bi+1/2(Wi+1 −Wi).(87) Observe that the terms that are multiplied with Pi+1/2in (73) are equal to the centered components of DF± i+1/2defined by (88) except the Coulomb friction term Ti+1/2∆x. That is, using (87), we could rewrite the numerical fluxes as DF± i+1/2=1 2R1,i+1/2±Pi+1/2R2,i+1/2, where, R2,i+1/2=R1,i+1/2−Ti+1/2∆x, and R1,i+1/2is defined by R1,i+1/2=F(cosθi+1/2,(cos2θ)i+1/2, Wi+1)− F(cosθi+1/2,(cos2θ)i+1/2, Wi)+ +S3,i+1/2(cos2θi+1 −cos2θi) + S4,i+1/2(cosθi+1 −cosθi)− −S1,i+1/2(bi+1 −bi)−S2,i+1/2(cos2θi+1 −cos2θi)−Bi+1/2(Wi+1 −Wi).(88) Let us prove i). Observe that in this case the first and third components of R1,i+1/2are equal to zero. The second component of R1,i+1/2is equal to [R1,i+1/2]2=gH2 1,i+1 2cosθi+1/2(cos2θ)i+1/2−gH2 1,i 2cosθi+1/2(cos2θ)i+1/2+ +3 4gH2 1,i+1/2cosθi+1/2(cos2θi+1 −cos2θi) + gH1,i+1/2cosθi+1/2(bi+1 −bi)+ +gH1,i+1/2 2(H1,i+1/2 2+ 2H2,i+1/2)cosθi+1/2(cos2θi+1 −cos2θi)+ +gH1,i+1/2cosθi+1/2(cos2θ)i+1/2(H2,i+1 −H2,i). We can write gH2 1,i+1 2cosθi+1/2(cos2θ)i+1/2−gH2 1,i 2cosθi+1/2(cos2θ)i+1/2= =gcosθi+1/2H1,i+1/2(cos2θ)i+1/2(H1,i+1 −H1,i), and then obtain [R1,i+1/2]2=gcosθi+1/2H1,i+1/2(cos2θ)i+1/2(H2,i+1 +H1,i+1 −(H2,i +H1,i))+ +g(H1,i+1/2+H2,i+1/2)(cos2θi+1 −cos2θi) + bi+1 −bi. Thanks to the definition of (cos2θ)i+1/2= (cos2θi+cos2θi+1)/2, we can use in the previous equation the following rule (discrete version of the derivative of a product), a b −c d =a+c 2(b−d) + (a−c)b+d 2,∀a, b, c, d ∈ R,(89) to obtain [R1,i+1/2]2=gcosθi+1/2H1,i+1/2(bi+1 +H1,i+1cos2θi+1 +H2,i+1cos2θi+1)− 24 −(bi+H1,icos2θi+H2,icos2θi). For i) we have constant free surface, b+H1cos2θ+H2cos2θ=cst, then the second component of R1,i+1/2is equal to zero. The fourth component of R1,i+1/2is equal to [R1,i+1/2]4=gcosθi+1/2H2,i+1/2Λ2(bi+1 −bi) + Λ2(cos2θ)i+1/2(H2,i+1 −H2,i)+ +Λ2H2,i+1/2(cos2θi+1−cos2θi)+rΛ1(H1,i+1/2(cos2θi+1−cos2θi)+(cos2θ)i+1/2(H1,i+1−H1,i))+ +(1 −Λ2)gH2 2,i+1/2 4cosθi+1/2(cos2θi+1 −cos2θi) + gH2,i+1/2cosθi+1/2(bi+1 −bi). Using (89) and b+H2cos2θ+H1cos2θ=cst we obtain [R1,i+1/2]4=gcosθi+1/2H2,i+1/2(Λ2−rΛ1)bi+1 −bi+H2,i+1cos2θi+1 −H2,icos2θi+ +(1 −Λ2)gcosθi+1/2H2,i+1/2H2,i+1/2 4(cos2θi+1 −cos2θi) + bi+1 −bi.(90) Finally, we conclude that the three first components of R1,i+1/2are zero. For the fourth component observe that [R1,i+1/2]4exactly coincides with (Tcrit,i+1/2∆x/cosθi+1/2) where Tcrit,i+1/2is defined by (80). Moreover, by hypothesis the given stationary solution verifies Qn 2,i = 0, and by (79) we obtain [Ti+1/2]4=Tcrit,i+1/2 cosθi+1/2 , Then, [R2,i+1/2]4= [R1,i+1/2]4−[Ti+1/2]4∆x= 0,⇒R2,i+1/2= 0. So, DF± i+1/2=R1,i+1/2. Additionally, in the second step we have Hn+1 1,i =Hn 1,i,Qn+1 1,i =Qn 1,i, Hn+1 2,i =Hn 2,i. Moreover, by (72) we deduce Q∗ 2,i =−∆t ∆x R1,i+1/2+R1,i−1/2 2= =−∆t 2 ∆x(Λ2−rΛ1)gH2,i+1/2cosθi+1/2(bi+1 +H2,i+1cos2θi+1)−(bi+H2,icos2θi)+ +(Λ2−Λ1r)gH2,i−1/2cosθi−1/2(bi+H2,icos2θi)−(bi−1+H2,i−1cos2θi−1)+ +(1 −Λ2)gcosθi+1/2H2,i+1/2H2,i+1/2 4(cos2θi+1 −cos2θi) + bi+1 −bi+ +(1 −Λ2)gcosθi−1/2H2,i−1/2H2,i−1/2 4(cos2θi−cos2θi−1) + bi−bi−1 Then, by (86) |Q∗ 2,i| ≤ ∆t(1 −r)gH2,i+1/2cosθi+1/2+H2,i−1/2cosθi−1/2 2tan(δ0) = σ∗ c,i∆t cosθi , where σ∗ c,i is defined by (85). Then, by (84) we obtain Qn+1 2,i = 0. Finally, we conclude that the stationary solutions verifying (52), (53), (54) satisfying (86) are exactly preserved. 25 2 3 4 5 6 7 8 −2 −1.95 −1.9 −1.85 −1.8 −1.75 −1.7 −1.65 −1.6 −1.55 −1.5 Water free surface bottom interfaces for 10, 11.31, 12 and 13 degrees 13 11.31 10 13 12 11.31 10 12 (a) Stationary landslide profiles for δ0equals to 10, 11.31 (angle of the plane), 12 and 13 degrees 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 −2 −1.9 −1.8 −1.7 −1.6 −1.5 −1.4 −1.3 −1.2 Water free surface bottom interfaces for 15, 20, and 25 degrees 25 20 15 25 20 15 (b) Statinary landslide profiles for δ0equal to 15, 20 and 25 degrees Figure 8: Landslide impinging a lake: Comparison of stationary landslide profiles for different values of δ0. 0 0.5 1 1.5 2 2.5 3 x 104 −2500 −2000 −1500 −1000 −500 0 500 Water free surface Sediment layer Bottom topography (a) Initial condition 0 0.5 1 1.5 2 2.5 3 x 104 −2500 −2000 −1500 −1000 −500 0 500 Water free surface Sediment layer Bottom topography (b) Final configuration of the landslide (t= 600s.) Figure 9: Tsunami experiment. (Horizontal x(m), Vertical z(m)) 32 that represents most of the tsunami energy. Over the landslide water is sucked, which creates a large trough splitting into two waves, one wave propagating shorewards and the other offshore. Later, the shoreward-propagating wave is followed by a positive wave, responsible for coastal inundation to which attention will be focused. Note that the second positive wave is of about 10m height when approaching the coastal line (see Figure 10(f)). Finally, Figure 9(b) shows the final position of the landslide. The results that we obtain are similar to those descrived by Heinrich et al. in [24]. Conclusions In this paper we have introduced a new model to study submarine avalanches and generated tsunamis. The presented model is a two-layer shallow water model including a Coulomb friction type term for the grain layer (Savage-Hutter model for the second layer). It is presented in local coordinates, by taking into account the curvature of a non-erodible bottom over which the avalanche and the tsunami flow. Some of the properties of the model are: the rank of the stationary solutions verifies exactly a dissipation entropy inequality, and the solutions of the model are solutions of Euler equations with hydrostatic pressure. The final system can be rewritten as a hyperbolic system with non-conservative products. We compare the model with that proposed by Heinrich et al., which is an uncoupled model, that mixes local coordinates for the evolution of the grain layer with non local coordinates for the evolution of the fluid layer. We see how the momentum equation for the grain layer of the Heinrich model is obtained under the assumption that the water surface is flat (rigid lid assumption). This model does not verify the properties of our model. We also present a numerical solver of finite volume type, based on a Roe method for hyperbolic systems with non-conservative products. With a special treatment of the Coulomb term, splitting the discretization of this term into an upwind explicit treatment and a second step to introduce an implicit centered discretization. This allows us to obtain a well-balanced solver for a wide rank of stationary solutions, when the angle of the slope of the grain layer is smaller than the corresponding angle of repose. The resulting model has been able to simulate sub-aerial and submarine avalanches and the generated tsunami by taking into account the interaction between the flowing mass and the surrounding water. Appendix: Change of variable In this Appendix, we perform the change of variable of Euler equations, from Cartesian coordinates ~ X= (x, z) to local coordinates (X, Z) (see Figure 11). We consider that the coordinate Zgives the position of an interior point ~ Xto the bed, measured in the normal direction to the bed. Thus 0< Z < S(t, x),with S(t, x) = h1(t, X) + h2(t, X).(92) Then, the relation between the Cartesian coordinates ~ Xand the coordinates (x, Z) related to the bed is ~ X=x−Zsinθ(x), b(x) + Zcosθ(x),(93) where (x, b(x)) is a point of the bed. ¯xis the x-Cartesian coordinate of the point (X, 0) (see Figure 11). We also consider a local variable X(x) measuring the arc length along the bed. We will denote by ∇~ Xand div ~ Xthe gradient and divergence operators in Cartesian coordinates. We consider the equations of conservation of mass and momentum in (x, z) Cartesian coordinates as follows ~ V=u v,∇ · ~ V= 0,(94) 33 0 0.5 1 1.5 2 2.5 3 x 104 −50 −40 −30 −20 −10 0 10 20 Water free surface Bottom topography (a) 30 s. 0 0.5 1 1.5 2 2.5 3 x 104 −50 −40 −30 −20 −10 0 10 20 Water free surface Bottom topography (b) 60 s. 0 0.5 1 1.5 2 2.5 3 x 104 −50 −40 −30 −20 −10 0 10 20 Water free surface Bottom topography (c) 90 s 0 0.5 1 1.5 2 2.5 3 x 104 −50 −40 −30 −20 −10 0 10 20 Water free surface Bottom topography (d) 150 s 0 0.5 1 1.5 2 2.5 3 x 104 −50 −40 −30 −20 −10 0 10 20 Water free surface Bottom topography (e) 300 s 0 0.5 1 1.5 2 2.5 3 x 104 −50 −40 −30 −20 −10 0 10 20 Water free surface Bottom topography (f) 480 s Figure 10: Tsunami evolution (zoom). (Horizontal x(m), Vertical z(m)) 34 XX ( )θ h Z h1 2 X b( ) x z x X Water surface Landslide surface (grain−water interface) layer 1 layer 2 Figure 11: Local coordinates. ∂t(ρ~ V) + ρ~ V∇~ X~ V=−∇ · P+ρ∇~ X(~g ·~ X),(95) where ~g = (0,−g), gbeing the constant acceleration of gravity. Moreover by Pwe denote the matrix of constraints, P=px x px z pz x pz z , (with px z =pz x). From (93) we have ∇(X,Z)~ X=Jcosθ−sinθ Jsinθcosθ, J = 1 −ZdXθ. Therefore, ∇~ X(X, Z) = 1 Jcosθsinθ −Jsinθ Jcosθ. The following result will be used in this Appendix: Lemma: Using the classical chain rule we have div(X,Z)(J∇~ X(X, Z)~ V) = Jdiv ~ X~ V , J =det(∇(X,Z)~ X),(96) and ∇~ XP= (∇~ X(X, Z))T∇(X,Z)Por (∇~ XP)T= (∇(X,Z)P)T∇~ X(X, Z).(97)  We will also use the following definitions: U W=cosθsinθ −sinθcosθ~ V . and P=cosθsinθ −sinθcosθPcosθ−sinθ sinθcosθ=PXX PXZ PZX PZZ . As pxz =pxz then PXZ =PZX . 35 Continuity equation (conservation of mass) In an incompressible material the velocity field ~ Vis solenoidal; so (94) holds. Multiplying (94) by Jand using (96) we have 0 = Jdiv ~ X~ V= div(X,Z)J∇~ X(X, Z)~ V= = div(X,Z) cosθsinθ −Jsinθ J cosθ~ V= div(X,Z)U J W . Thus, we have the equation, ∂X(U) + ∂Z(J W ) = 0.(98) This correspond to the first equation of (4) Conservation of momentum In this section, we first find the equation for Uand then for W. Equation for U. We add the first component of equation (95) multiplied by cosθto the second component of equation (95) multiplied by sinθ. Then, we obtain ρ∂tU+ρdiv ~ X(U~ V) + ρ(∇~ X(~g ·~ X))Tcosθ sinθ= −div~ XPcosθ sinθ+ρWdiv ~ X(θ~ V) + (∇~ Xθ)TP−sinθ cosθ. In order to apply the rules (96) and (97) we rewrite this last equation as ρJ∂tU+ρJdiv ~ XUu Uv +Jρ(∇~ X(~g ·~ X))Tcosθ sinθ= −Jdiv ~ XPcosθ sinθ+ρJ Wdiv ~ Xθ~ V+J(∇~ X·θ)TP−sinθ cosθ. Then, using the rules (96) and (97), we obtain ρ∂t(JU) + ρdiv(X,Z)U2 JUW +ρ(∇X,Z (~g ·~ X))T1 0= =−div(X,Z)1 0 0J cosθsinθ −sinθcosθPcosθ sinθ+ +ρWdiv(X,Z)θU JθW + +(∇(X,Z)θ)T1 0 0J cosθsinθ −sinθcosθP−sinθ cosθ. This correspond to the second equation of (4). Equation for W: Now, multiplying the first equation of (95) by (−sinθ) and the second one by (cosθ) and adding them, we obtain ρ∂tW+ρdiv ~ X(W~ V) + ρ(∇~ X(~g ·~ X))T−sinθ cosθ= 36 −div ~ XP−sinθ cosθ−ρUdiv ~ X(θ~ V)−(∇~ Xθ)TPcosθ sinθ. In order to apply the rules (96) and (97), we rewrite this last equation as ρ∂t(JW ) + Jρdiv ~ XWu Wv +Jρ(∇~ X(~g ·~ X))T−sinθ cosθ= −Jdiv ~ XP−sinθ cosθ−ρJ Udiv ~ Xθ~ V−J(∇~ Xθ)TPcosθ sinθ. Then, using the rules (96) and (97), we obtain ρ∂t(JW ) + ρdiv(X,Z)W U JW 2+ρ(∇X,Z (~g ·~ X))T0 J= =−div(X,Z)1 0 0J cosθsinθ −sinθcosθP−sinθ cosθ− −ρUdiv(X,Z)θU JθW − −(∇(X,Z)θ)T1 0 0J cosθsinθ −sinθcosθPcosθ sinθ. This correspond to the third equation of (4). Acknowledgments. We thank Lev Tsimring, Nicolas Ledante, Dmitri Volfson and Carlos Par´es for interesting discussions on the application of the model. We also thank Anael Lemaitre by his interesting discussion about the definition of the stress tensor for the grain-fluid mixture layer. We thank Kolumban Hutter and the other anonymous reviewer for their delaited comments that has helped to improve the paper. This research has been partially supported by the Spanish Government Research projects MTM2006-08075 and MTM2006-01275. By the ACI Nouvelles Interface de Mathematiques (CNRS), ACI Jeunes Chercheurs (CNRS), ANR Blanche, BLAN-06-1-140039, ACI NIM contract no 03318 and the ANR contract no ANR-06-BLAN-0414, ACI jeunes chercheurs ”Analyses math´ematiques de param´etrisations en oc´eanographie”, Project r´egion Rhˆone-Alpes ”Mod´elisation d’avalanches”. References [1] C. Ancey, Plasticity and geophysical flows: A review J. Non-Newtonian Fluid Mech (2006). [2] A. Aradian, E. Raphael, P.G. de Gennes, Surface flow of granular materials: A short introduction to some recent models. C.R. Physique 3, 187-196 (2002). [3] I.S. Aranson, L. S. Tsimring, Continuum theory of partially fluidized granular flows. Phys. Rev. E. vol. 65, 061303 (2002). [4] G.S. Beaves, D.D. Joseph, Boundary conditions at a naturally permeable wall, J. Fluid Mech., 30, part I, 197–207 (1967). [5] F. Bouchut, E.D. Fern´andez-Nieto, A. Mangeney, P.Y. Lagree On new erosion models of Savage-Hutter type for avalanches. Preprint (2007). 37 [6] F. Bouchut, A. Mangeney-Castelnau, B. Perthame, J.P Vilotte, A new model of Saint Venant and Savage-Hutter type for gravity driven shallow flows. C.R. Acad. Sci. Paris, Ser I 336 531-536 (2003). [7] F. Bouchut, M. Westdickenberg, Gravity driven shallow water model for arbitrary topography. Comm. Math. Sci. 2(3), 359-389 (2004). [8] C. Cassar, M. Nicolas, O. Pouliquen, Submarine granular flows down inclined planes. Physics of Fluids, 17, 103301, DOI 10.1063/1.2069864, (2005). [9] M.J. Castro, T. Chac´on, E.D. Fern´andez-Nieto, C. Par´es On well-balanced finite volume methods for non-conservative non-homogeneous systems. SIAM J. Sci. Comput. 29(3): 10931126, (2007). [10] M.J. Castro, A.M. Ferreiro, J.A. Garc´ıa, J.M. Gonz´alez, J. Mac´ıas, C. Par´es, M.E. V´azquez, On the numerical treatment of wet/dry fronts in shallow flows: application to one-layer and two-layer systems. Match Comp. Model. 42 (3-4): 419-439 (2005). [11] M.J. Castro, J.A. Garc´ıa, J.M. Gonz´alez and C. Par´es, A parallel 2d finite volume scheme for solving systems of balance laws with nonconservative products: application to shallow flows. Comp. Meth. Appl. Mech. Eng., 196, 2788-2815 (2006). [12] T. Chac´on, A. Dom´ınguez, E.D. Fern´andez. A family of stable numerical solvers for Shallow Water equations with source terms. Comp. Meth. Appl. Mech. Eng. 192: 203–225, 2003. [13] T. Chac´on, A. Dom´ınguez, E.D. Fern´andez. Asymptotically balanced schemes for non-homoegeneous hyperbolic systems – application to the Shallow Water equations. C.R. Acad. Sci. Paris, Ser. I 338: 85–90, 2004. [14] G. Dal Maso, P.G. LeFloch, F. Murat. Definition and weak stability of nonconservative products. J. Math. Pures Appl. 74: 483–548, 1995. [15] R.A. Dalrymple, S.T. Grilli, J.T. Kirby. Tsunamis and challenges for accurate modeling. Oceanography, 19:142:151, (2006). [16] B. Di Martino, F. Flori, C. Giacomoni, P. Orenga. Mathematical and Numerical Analysis of a Tsunami Problem. Math. Models and Methods in Apllied Sciences. Vol. 13, p. 1489-1514, (2003). [17] D. Dutykh and F. Dias. Water waves generated by a moving bottom. In Anjan Kundu, editor, Tsunami and nonlinear waves, Springer-Verlag (Geo Sc.), 38, 64–68, (2007). [18] D. Dutykh, Mod´elisation math´ematique des Tsunamis. Th`ese de doctorat de l’Ecole Normale Sup´erieure de Cachan, (2007). [19] I.V. Fine, A.B. Rabinovich, R.E. Thomson, E.A. Kulikov, Numerical Modeling of Tsunami generation by submarine and subaerial landslides. Submarine Landslides and Tsunamis 6988, Kluwer Acadmic Publishers (2003). [20] D.L. Geroge, R.J. LeVeque, Finite Volume methods and adaptive refinement for global tsunami propagation and local inundation. Science of Tsunami Hazards, 24, 319–329 (2006). [21] J.M.N.T. Gray, Granular flow in partially filled slowly rotating drums. J. Fluid Mech, 441, 1-29 (2001). [22] Ph. Heinrich, A. Mangeney, S. Guibourg, R. Roche, G. Boudon, J.L. Cheminee, Simulation of water waves generated by a potential debris avalanche in Montserrat, Lesser Antilles, Geophys. Res. Lett., 25(19), 3697-3700 (1998). 38 [23] Ph. Heinrich, S. Guibourg, A. Mangeney, R. Roche, Numerical modeling of a landslidegenerated tsunami following a potential explosion of the Montserrat volcano, Phys. Chem. Earth, 24(2), 163-168 (1999). [24] Ph. Heinrich, A. Piatanesi, H. H´ebert, Numerical modelling of tsunami generation and propagation from submarine slumps: the 1998 Papua New Guinea event. Geophys. J. Int. 145, 97-11, (2001) [25] K. Hutter, K. J¨ohnk, B. Svendsen, On interfacial transition conditions in two phase gravity flow, Z angew Math Phys. 45 (1994). [26] R.M. Iverson, R.P. Denlinger, Flow of variably fluidized granular masses across threedimensional terrain. J. of Geoph. Res. 106, B1, 537-552, (2001). [27] D.V. Khakhar, A.V. Orpe, P. Andres´en, J.M. Ottino, Surface flow of granular materials: model and experiments in heap formation. J. Fluid Mech. 441, 225-264 (2001). [28] Locat, J., and Mienert, J (Eds.), Submarine Mass Movements and Their Consequences, Kluwer Academic Publishers, Dordrecht, 540 pp (2003). [29] Lucas, A., and Mangeney, A., Mobility and topography effects for large Valles Marineris landslides on Mars, Geophys. Res. Lett., 34, L10201, doi:10.1029/2007GL029835. [30] A. Mangeney, Ph. Heinrich, R. Roche, G. Boudon, J.L. Cheminee, Modeling of debris avalanche and generated water waves: Application to real and potential events in Montserrat, Phys. Chem. Earth, 25(9-11), 741-745, (2000). [31] A. Mangeney-Castelnau, J.P. Vilotte, M.O. Bristeau, B. Perthame, F. Bouchut, C. Siomeoni, S. Yerneni, Numerical modeling of avalanches based on Saint Venant equations using a kinetic scheme. J. Geoph. Res. vo. 108, n. B11, 2527 (2003). [32] A. Mangeney-Castelnau, F. Bouchut, T.P. Vilotte, E. Lajeneusse, A. Aubertin, M. Pirulli, On the use of Saint-Venant equations to simulate the spreading of a granular mass, J. Gephys. Res. 110 (B9), B09103 (2005). [33] A. Mangeney-Castelnau, F. Bouchut, N. Thomas, J.P. vilotte, M.O. Bristeau, Numerical modelling of self-channeling granular flows and of their leveel channel deposists, J. Geophs. Res.(SE), 112, F2017, 21pp. (2007). [34] T. Morales de Luna, A Saint Venant model for gravity driven shallow water flows with variable density and compressibility effects. Math. Comput. Modelling 47, 3-4, 436-444 (2008). [35] N. Nomanbhoy, K. Satake, Generation mechanism of tsunamis from the 1883 Krakatau eruption, Geopys. Res. Lett., 22, 509-512 (1995). [36] C. Par´es, M.J. Castro 2004, On the well-balance property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow-water systems. ESAIM: M2AN, 38(5):821– 852 (2004). [37] M. Pelanti, F. Bouchut, A. Mangeney, A Roe-Type Scheme for Two-Phase Shallow Granular Flows with Bottom Topography,in press in Math. Model. Numeric. Analy. (ESAIM:M2AN) (2008). [38] M. Pirulli, M.O. Bristeau, A. Mangeney, A and C. Scavia, The effect of the earth pressure coefficients on the runout of granular material, Environ. Modell. and Soft., 22(10), 14371454 (2007). 39 [39] M. Pirulli and A. Mangeney, Result of Back-Analysis of the Propagation of Rock Avalanches as a Function of the Assumed Rheology, Rock Mech. Rock Engng., 41(1), 59-84 (2008). [40] E. B. Pitman, l. Le, A two-fluid model for avalanche and debris flows, Phil. Trans. R. Soc. A, 363:1573-1601 (2005). [41] Pudasaini, S., and Hutter, K., Avalanche Dynamics, Springer-Verlag, 602p (2007). [42] Pudasaini, S., Wang Y., Hutter K., Modelling debris flows down general channels. Natural Hazards and Earth System Sciences, 5, 799-819 (2005). [43] M. Ricchiuto, R. Abgrall, H. Deconinck, Appliation of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes. J. Comp. Physics 222, 287-331 (2007) [44] J. Sander, K. Hutter, Experimental and computational study of channelized water waves generated by a porous body, Acta. Mech., 115, 133-149 (1996). [45] S.B. Savage, K. Hutter, The dynamics of avalanches of granular materials frominitiation to run-out, Acta Mech. 86, 201-223 (1991). [46] J.B. Schijf and J.C. Schonfeld, Theoretical considerations on the motion of salt and fresh water, in. Proc. Minn. Int. Hydraulics Conv., joint meeting IAHR Hydro. Div. ASCE., 321–333 (1953). [47] B. Svendsen, T. Wu, K. Johnk, K. Hutter, On the Role of Mechanical Interactions in the Steady-State Gravity Flow of a Two-Constituent Mixture down an Inclined Plane, Porc. R. Soc. Lond., Vol 452, 1189-1205 (1996). [48] Vanneste, M., Mienert, J., and Bunz, S., THe Hinlopen Slide: A giant, submarine slope failure on the northern Svalbard margin, Arctic Ocean, Earth Plan. Sci. Lett., 245, 373-388, (2006). [49] M. Wieland, J.M.N.T. Gray, K. Hutter, Channelized free-surface flow of cohesionless granular avalanches in a chute with shallow lateral curvature. J. Fluid Mech. 392, 73-100 (1999) [50] T. Wu and K. Hutter, On the Role of the Interface Mechanical Interaction in a GravityDriven Shear Flow of an Ice-Till Mixture. Transport in Porous Media, vol 34, 3-14 (1999) [51] T. Wu, K. Hutter, B. Svendesen, On shear flow of a satured ice-sediment mixture with thermodynamic equilibrium pressure and momentum exchange. Proc. R. Soc. Lond., A454, 71-88 (1998) 40