On new erosion models of Savage-Hutter type for avalanches

On new erosion models of Savage-Hutter type for avalanches F. Bouchut∗ , E.D. Fern´andez-Nieto† , A. Mangeney‡& P.-Y. Lagr´ee§ October 22, 2007 Abstract In this work we study the modeling of one-dimensional avalanche flows made of a moving layer over a static base, where the interface between the two can be time dependent. We propose a general model, obtained by looking for an approximate solution with constant velocity profile to the incompressible Euler equations. This model has an energy dissipation equation that is consistent with the depth integrated energy equation of the Euler system. It has physically relevant steady state solutions, and, for constant slope, it gives a particular exact solution to the incompressible hydrostatic Euler equations. Then, we propose a simplified model, for which the energy conservation holds only up to third-order terms. Its associated eigenvalues depend on the mass exchange velocity between the static and moving layers. We show that a simplification used in some previously proposed models gives a nonconsistent energy equation. Our models do not use, nor provide, any equation for the moving interface, thus other arguments have to be used in order to close the system. With special assumptions, and in particular small velocity, we can nevertheless obtain an equation for the evolution of the interface. Furthermore, the unknown parameters of the model proposed by Bouchaud, Cates, Ravi and Edwards (BCRE model [7]) can be derived. For the quasistationary case we compare and discuss the equation for the moving interface with Khakhar’s model [13]. ∗D´epartement de Math´ematiques et Applications, ´ Ecole Normale Sup´erieure & CNRS, 45 rue d’Ulm, 75230 Paris cedex 05, France ([email protected].fr) †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 Mod´elisation Physique et Num´erique, IPGP, 4 pl. Jussieu, 75232 Paris cedex 05, France ([email protected]ussieu.fr) Now at Institut for NonLinear Science, UCSD, 9500 Gilman Drive, La Jolla, CA 92093-0402, USA §Laboratoire de Mod´elisation en M´ecanique UMR 7607k, 4 pl. Jussieu, 75232 Paris cedex 05, France ([email protected]) 1 1 Introduction This work focuses on the mathematical modeling of erosion processes in granular flows. Granular avalanches generally involve flowing zones and zones at rest. Static particles can be put into motion by the flowing particles (i.e. erosion) whereas flowing particles may suddenly stop (i.e. deposition). This static/flowing transition is typical of granular matter. Many industrial or geophysical applications are related to this problem. For example, when granular materials are mixed for industrial purposes, the presence of static and flowing zones could significantly affect the efficiency of the mixing process. This is a major issue in particular when rotating cylinders are used to deal with mixing of cereals or pills. In geophysics, erosion processes are expected to play a key role in the dynamics of snow or rock avalanches in mountainous and volcanic context [Sovilla et al., 2006; Pudasaini and Hutter, 2007]. Mangeney et al., 2007b show that erosion processes may change drastically the dynamics of an avalanche flowing over an erodible bed. Actually, when entrainment of granular material lying on the erodible bed is effective, a decelerating avalanche could become a surge moving at constant velocity along the slope. The theoretical description and physical understanding of these processes are still challenging open problems. In some applications, e.g. avalanches in silos or surface flows over a pile of granular material, the no-flow and flowing zones can be fairly clearly separated into a granular layer flowing over a static layer of grains. The interface between the static and the flowing material changes with time due to the exchange of grains between the static and the flowing layers. In the numerical modeling of avalanche dynamics, the existence of mass exchange between the no-flow and flowing zones is generally not taken into account. The models used in geophysics are based on the pioneering work of Savage and Hutter [22] where the whole column of granular material is assumed to be flowing. The mass and momentum equations are averaged over the depth of the granular material and a scaling analysis is performed with respect to the aspect ratio of the flowing mass, considered to be small. Under this hypothesis the avalanche is modeled by a Saint Venant type system derived in a reference frame linked to an inclined plane. The unknowns are the thickness and the mean (depth-averaged) velocity of the flowing mass. A constitutive relation based on Mohr-Coulomb plasticity theory is imposed, making it possible to relate the normal stresses through a coefficient involving the so-called basal and internal friction angles and to deduce a friction term at the base of the flow. The resulting system, derived for granular flows over an inclined plane, will be denoted here as the S-H model (for Savage-Hutter). Wieland et al. [24] generalized the S-H model for granular flows over surfaces with small lateral curvature (Pudasaini and Hutter [21]). More recently Bouchut et al. proposed two more general one-dimensional Savage-Hutter models [5]. The first one is valid for small variations of the local slope angle and makes it possible to derive an exact energy equation. The second model is developed for general slopes. In [6], Bouchut and Westdickenberg generalize the previous models for small or for general slope 2 variation in two dimensions. The S-H model has been successfully applied to simulate a wide range of laboratory experiments [21], as complex as well as self-channeling flows and the leveechannel morphology of their deposit [19]. However when applied to the collapse of granular columns [18], several limitations of the model have been detected. Actually, experimental results and discrete element simulations show that a static/flowing interface appears and changes with time during the collapse. Taking into account this static/flowing transition in a depth-averaged model requires the determination of an equation for the evolution of the moving interface. In other words, one has to specify the entrainment/deposition rate, which is the velocity perpendicular to the interface of mass that changes from the static to the flowing side or vice versa of the interface. Its parametrization is insufficiently known owing to the currently unknown appropriate physical or mechanical processes that govern the static/flowing transition. In order to avoid prescribing the evolution of the interface by ad-hoc physical or mechanical considerations, several attempts have been made to close the equations by imposing a given profile of the horizontal velocity and/or a given velocity at the static/flowing interface. The purpose of this work is to derive mathematically accurate S-H type models for a one-dimensional flow of a layer made of a granular material moving over a static layer made of the same material. The closure relation making it possible to derive an equation for the static/flowing interface will be discussed and compared with different previously proposed models ([7], [8], [2], [13], [12]). In particular, the existence of an energy equation for the different systems will be investigated. The proposed models will be discussed using the following notation (see Figure 1). The horizontal coordinate is denoted by xand time by t. By b(t, x) we denote the vertical thickness of the static layer (i. e. the position of the static/flowing interface), by h(t, x) the thickness of the moving layer, measured perpendicularly to the unknown interface, and by θ(t, x) the angle between the horizontal and the tangent vector to the interface, according to Figure 1 (indeed, tan θ=∂xb). The variable Zis a coordinate normal to the interface, attached to a point ~ Xin the moving layer. One of the first models devoted to this erosion problem, proposed by Bouchaud, Cates, Ravi and Edwards (BCRE model, see [7], [8]), is simply based on the principle of mass conservation for density preserving materials. The variation of his driven by an advection equation with a volume exchange E(t, x) between the static and the flowing layer, while the time variation of bcompensates the volume exchange, and is equal to −E. The volume exchange E(t, x) and the horizontal velocity component of the moving grains, Vd, are both parameters of the model,        ∂th cos θ+∂x(hVd) = E(t, x), ∂tb=−E(t, x). (1.1) 3 X x θ z Z b( )t, x h db x dX d t, x ( ) X Figure 1: Interface between the static and the moving layers, and change of variable Physically, E(t, x) is a volume exchange per unit time and unit horizontal line element. Although the BCRE model does not ensure conservation of momentum, it has been applied to a lot of practical problems. A difficulty of this model is to adjust the parameters E(t, x) and Vd. Phenomenological relations have been proposed for E(t, x), defined as a function of the angle θminus a neutral angle θn(if θ > θnthen erosion occurs, deposition otherwise), see [7]. The exchange term reads E(t, x) = γh(θ−θn),(1.2) where γis an empirical constant (see also [1], [2]). The model proposed by Boutreux, Rapha¨el and DeGennes (BRDG model, see [8]) differs in the form of the exchange term, where his replaced by a thickness λsmaller than h, E(t, x) = γλ(θ−θn).(1.3) The idea is that only a very thin layer, of the order of the grain diameter, is involved in the exchange between the static and moving layers. Both expressions (1.2) and (1.3) can be written as E(t, x) = Vup(θ−θn),(1.4) where Vup is a parameter of the model that has the dimension of a velocity. In [2] the BCRE and BRDG models are derived from the Saint-Venant equations, making it possible to deduce the parameters involved in these models when the velocity profile is assumed to be linear with a constant slope. From this assumption, an equation for the change in time of band hcan be deduced. Concretely, Vd is simply obtained from the mass conservation equation when a velocity profile is assumed and Vup is shown to depend on the specific velocity profile when a quasistationary assumption is made (θ≈θn). Khakhar et al. in [13], [14] derive a model based on the Saint Venant equations by also assuming a linear velocity profile. The 4 equations are closed by assuming continuity of the shear stress at the static/flowing interface (this is only in conformity with the jump conditions of momentum, if E ≡ 0). Based on experimental data, the friction term in the flowing layer at the interface is written as the sum of a Coulomb term plus a Bagnold term whereas the friction is assumed to reduce to the Coulomb friction in the static layer at the same interface. From the equality of both expressions the value of the slope of the velocity profile can be deduced, U=U1Zwith U1=f(θ, θs, θm),(1.5) where θsis the static angle of repose and θmthe ”maximum angle of repose” (see [13]). Note that in [13] and [2], the pressure term ∂x(gh2/2) arising in the Saint Venant systems has been neglected. In [12], Gray deduces a S-H type model for the erosion problem in avalanches in the case of rotating drums. He uses similar arguments as in the original S-H model, but by depth-averaging the equation only within the moving layer. The system of coordinates is fixed over an inclined plane with a fixed angle. For the kinematic condition he supposes that the particles on the surface that separates the moving and the static layers move with a velocity which is the sum of the velocity of the particles in the moving layer and an exchange velocity (db) in the normal direction to the surface. From this kinematic condition he deduces the exchange velocity db for granular flows in rotating drums. For avalanches over an inclined plane one has db= cos θ ∂tb. After depth integration, a source term depending on dbappears in the mass balance equation. Moreover, a term proportional to U|b∂tbarises in the momentum equation, where U|bis the velocity of the grains in the moving layer at the static/flowing interface. This term corresponds to the impulse of the entrained mass that must instantly assume the velocity of the moving avalanche at the interface. In [13] the authors suppose that the velocity profile is linear with the velocity equal to zero at the static/moving interface, i.e. U|b= 0. Instead, by referring to experimental data, Gray in [12] sets U|b=u, the averaged velocity over the moving layer. The question is whether these closure arguments, in particular using the velocity profile and/or the velocity U|b, are justified. A possible way to address this problem is to look for an approximate solution to the Euler equations that satisfies the required conditions. When simple solutions cannot be calculated, another criterion is to check that the assumptions lead to a system of equations which exhibits good mathematical properties, i.e. to have the right steady states and to satisfy an energy dissipation inequality. The main model that we propose, (3.8)-(3.9) or its simplified version (3.20)-(3.21), is coherent with these two criteria and is based on the assumption of a uniform velocity profile. We have not been able to offer such a justification for a model that is based on the assumption of a linear profile of the velocity. In Subsection 3.6 we nevertheless discuss the energy equation obtained for the different models by introducing a coefficient δin front of the term u∂tbthat appears in the momentum equation in conservative form. The value of 5 this coefficient depends on the assumptions related to the velocity profile and/or the velocity at the static/flowing interface. The associated energy equation is shown to be fully energy conservative (i.e. has a neutral balance independently of b(t, x)) only if δ= 1/2. However, the only value that we have been able to justify in terms of mathematical properties is δ= 1. Moreover, we show that if the term u∂tbin the momentum equation is neglected (i.e. δ= 0), as is for most models proposed in the literature, the associated energy equation is not consistent with the depth-averaged energy equation of the Euler equations. In order to close depth-averaged avalanche models, the ”partial fluidization model” proposed in [3] could be a good candidate avoiding the difficulties encountered in the methods discussed above. Actually, Aranson and Tsimring present a theory based on the introduction of an order parameter ρthat describes the transition between the static and the moving phases without any particular assumption on the velocity profile and/or the velocity at the static/flowing interface. This parameter is introduced in the hydrodynamic equation and a coupled equation for ρ is proposed following the theory of phase transition of Ginzburg and Landau, as presented by Landau and Lifshitz in [16]. Based on this theory, Mangeney et al. in [20] succesfully simulate laboratory experiments of granular flows over erodible bed. The remainder of the paper is organized as follows. In Section 2 we perform a change of variables, starting from the incompressible Euler equations over a moving bottom. Actually, two variable transformations are needed in order to avoid the difficulty that the variation in time of the local coordinates is not known. Moreover, the mass and energy equations associated to the Euler system are written down. In Section 3 we show how an approximate solution to the Euler system can be built in new coordinates, by assuming a constant velocity profile. We deduce a new model that satisfies a conservative energy inequality without any error term, in the sense that for regular solutions the energy equation is exactly obtained by a combination of the equations of the model. Moreover, it gives the expected steadystates corresponding to zero velocity (material at rest). For constant slope, the model gives an exact solution to the hydrostatic Euler equations. Friction terms can also be introduced in the model. An interesting observation is that this new model contains first-order terms that do not appear in previously proposed models. Then, we introduce a simplified model which is energy conservative only up to thirdorder terms in the aspect ratio. Its associated characteristic velocities depend on the exchange velocity between the layers ∂tb. Finally, the developed models are compared to previously considered models in Section 4. The arguments proposed in the former studies are used in our model to derive an evolution equation for the interface b. In particular, we show that the parameters of the BCRE model found by [2] through the derivation of the Saint-Venant equations change by a factor 3/2 if the term U|b∂tb, usually neglected in the momentum equation, is taken into account. 6 2 Change of variables In this section we perform a change of variables from Cartesian coordinates to the coordinates related to the interface. This enables us to reformulate the incompressible Euler equations in a more suitable form. We also reformulate the energy equation associated to the Euler equations in order to be able to compare the integrated energy with the energy associated to any integrated model. We consider a one-dimensional flow comprising a moving layer over a static one, with a time-dependent interface between them. Time is denoted by t, the horizontal coordinate is denoted by x,b(t, x) is the vertical thickness of the static layer, h(t, x) is the thickness of the moving layer, measured perpendicularly to the unknown interface. The angle between the horizontal and the tangent vector to the interface is denoted by θ(t, x), and we notice that tan θ=∂xb. (2.1) Let ~ Xbe an arbitrary point in the avalanche domain bounded by the basal and free surfaces. We consider the coordinate Zthat measures the position of ~ Xinside the moving layer, in the direction normal to the interface, as displayed on Figure 1. Thus, 0< Z < h(t, x).(2.2) Then the relations between the Cartesian coordinates (t, ~ X) and the coordinates (t, x, Z) related to the interface is ~ X=x−Zsin θ(t, x), b(t, x) + Zcos θ(t, x).(2.3) Notice that (x, b(t, x)) is a point on the interface. We consider also a curvilinear variable X(t, x) measuring the arc length along the interface. Thus at fixed t, the mapping x7→ Xis a change of variable with ∂xX=1 cos θ(t, x),(2.4) which, together with (2.1) yields the relation (at fixed t) (db/dx)2+1 = (dX/dx)2, as illustrated in Figure 1. However, since there are several ways to choose the origin on the curve, that could make different dependencies with respect to time, we cannot say a priori what is the value of ∂tX. Nevertheless, from (2.1) we get ∂tθ= cos2θ ∂2 txb, ∂Xb|t= sin θ, (2.5) where the notation |tmeans that the derivative is taken at fixed t. 7 2.1 Transformation of the incompressible Euler equations We denote the velocity field of the moving layer by ~ U. Then the Euler equations for the moving layer can be written in Cartesian coordinates (t, ~ X) as ∇~ X·~ U= 0,(2.6) b ∂t~ U+~ U·∇~ X~ U+∇~ X·P=−∇~ X(~g ·~ X),(2.7) where b ∂t=∂t|~ X,~g = (0, g) and gis the gravity constant. By Pwe denote the pressure tensor, the negative of the stress tensor, divided by density, P=px x px z pz x pz z ,(2.8) with px z =pz x. Here we allow a full matrix, but later on we shall consider only the ”true Euler” case where Pis isotropic, e.g. Pis a scalar times the identity tensor. This system is completed with a kinematic law for the evolving free surface, the free surface is advected by the material velocity ~ U, (2.9) and with dynamic boundary conditions. At the free surface we just set P~ν = 0 at Z=h(t, x),(2.10) where ~ν is a unit vector normal to the free surface. For the bottom interface, the incompressibility of the fluid comes into play, and we have to give precise assumptions on the fluid density, that we have not introduced yet. We assume that the moving part of the fluid is incompressible with uniform density ρm, while the part of the fluid that is at rest is also incompressible, with uniform density ρr. Having different densities ρmand ρr, the model allows more generality, and it is physically more relevant, since experimental observations indicate that ρr> ρm. The physical interpretation of these assumptions is the following. On the one hand, when ρm=ρr, the whole fluid (comprising the moving part and the part at rest) is incompressible with uniform density. On the other hand, when ρr> ρm, the fluid that is deposited (e.g. stops moving and stays at rest) has to increase its density from ρmto ρr, thus is being compressed through the interface; while the fluid that is eroded (e.g. starts moving), has to diminish its density from ρrto ρm, thus is being decompressed through the interface. The dynamic condition at the interface corresponding to this situation with two densities ρm,ρris deduced from the mass conservation, that can be written b ∂tρ+∇~ X·(ρ~ U) = 0,(2.11) where ρtakes the value ρmin the moving part of the fluid, and ρrin the part at rest. The velocity ~ Uis set to 0 in the part at rest. Writing down the Rankine-Hugoniot condition through the interface gives the boundary condition [ρ],[ρ~ U]·~ N= 0,(2.12) 8 where [...] denotes the jump of a quantity through the interface, and ~ Ndenotes a time-space normal to the interface. Since the interface is defined by z−b(t, x) = 0, one can take ~ N= (−∂tb, −∂xb, 1) = (−cos θ ∂tb, −sin θ, cos θ)/cos θ. Denoting by ~n = (−sin θ, cos θ) the space unit vector normal to the interface, oriented upwards, the interface dynamic boundary condition (2.12) can be written ρm~ U·~n =−(ρr−ρm) cos θ ∂tbat Z= 0.(2.13) For completeness, let us mention that from (2.6)-(2.7) one classically deduces the energy equation, valid if Pis isotropic, e.g. Pis a scalar times the identity tensor, b ∂t |~ U|2 2+~g ·~ X!+∇~ X· |~ U|2 2~ U+P~ U+ (~g ·~ X)~ U!= 0.(2.14) In what follows we perform a change of variables from Cartesian coordinates (t, ~ X) to coordinates attached to the interface, first to the variables (t, X, Z), and then to the variables (t, x, Z). The velocity components in the frame associated to the interface will be denoted by (U, W ), U W=cos θsin θ −sin θcos θ~ U. (2.15) The new stress tensor is P=cos θsin θ −sin θcos θPcos θ−sin θ sin θcos θ=PXX PXZ PZX PZZ .(2.16) We observe that, as pxz =pzx, then PXZ =PZX . Let us consider first the variables (t, X, Z). In order to avoid confusion in differentiations, we shall denote the time in these variables by τ, so that ∂/∂τ means that we differentiate with respect to time at Xand Zfixed. We first compute the Jacobian matrix of the transformation ∇τ,X,Z(t, ~ X). According to (2.3), (2.4) and (2.5) we have ∇X,Z ~ X=Jcos θ−sin θ Jsin θcos θ, J = 1 −Z∂Xθ= det ∇X,Z ~ X.(2.17) Thus, ∇~ X(X, Z) = (∇X,Z ~ X)−1=1 Jcos θsin θ −Jsin θ J cos θ.(2.18) Then ∇τ,X,Z(t, ~ X) = 1 0 B∇X,Z ~ X, B =∂τ~ X, (2.19) and by inversion ∇t, ~ X(τ, X, Z) = 1 0 −(∇X,Z ~ X)−1B(∇X,Z ~ X)−1.(2.20) 9 A similar computation for W, again with the aid of Lemma 2.1, transforms (2.52) into ∂tJW cos θ+∂x(J(b ∂tX−∂tX) + U)W+PXZ  +∂ZJ cos θ(b ∂tZ+W)W+JPZZ cos θ+J cos θ∂Zg(b+Zcos θ) =−U∂x(cos θ ∂tb) + U∂xθ−PXX ∂xθ. (2.61) If we combine this again with (2.58) as before and use (2.53), we obtain the nonconservative equation J cos θ∂tW+ (Zcos2θ ∂2 txb−sin θ ∂tb+U)∂xW+JW cos θ−∂tb∂ZW +∂xPXZ +J cos θ∂Zg(b+Zcos θ) + PZZ =−U∂x(cos θ∂tb) + U∂xθ+ (PZZ −PXX )∂xθ. (2.62) In summary, the governing equations, comprising all the balance laws of mass and momentum and referred to the independent variables (t, x, Z) is defined: (i) in conservative form by (2.54), (2.59), (2.61), and (ii), in nonconservative form by (2.54), (2.60), (2.62). 2.6 Energy equation In what follows we consider the case where Pis isotropic, P=pzzI. Then we have also P=pzzI, and for simplicity we shall denote pzz by P. Starting from the energy equation (2.14) written in Cartesian coordinates ~ Y= (t, ~ X), we perform the change of variables to ~ ξ= (τ, X, Z). Applying Lemma 2.1 with (2.38) and appropriate choices for ~ Φ leads to ∂τJU2+W2 2+~g ·~ X+∂X(Jb ∂tX+U)U2+W2 2+~g ·~ X+PU +∂ZJ(b ∂tZ+W)U2+W2 2+~g ·~ X+PJW= 0. (2.63) Then, applying Lemma 2.1 again with ~ Y= (τ, X, Z), ~ ξ= (t, x, Z), and using (2.55), we obtain the balance equation for the energy referred to the (t, x, Z) coordinates as follows, ∂tJ cos θU2+W2 2+g(b+Zcos θ) +∂x(Zcos2θ ∂2 txb−sin θ ∂tb+U)U2+W2 2+g(b+Zcos θ)+PU +∂ZJ cos θ(−cos θ ∂tb+W)U2+W2 2+g(b+Zcos θ)+PJW cos θ= 0. (2.64) 16 2.7 Integration and kinematic boundary condition In this subsection we provide formulas for quantities which are integrated with respect to Z. Moreover, we derive the kinematic boundary condition associated to the free surface, that gives the evolution of h(t, x). According to (2.9), the free surface is advected by the material velocity ~ U. Consider a material point ~ X(s), where sis the time, that is moving with the flow, and that stays on the free surface. Then d~ X ds (s) = ~ U(s, ~ X(s)).(2.65) In the variables (t, x, Z), the free surface is characterized by Z=h(t, x). Thus Z(s, ~ X(s)) = h(s, x(s)),(2.66) where (s, x(s)) are the (t, x) coordinates corresponding to (s, ~ X(s)). Differentiating (2.66) with respect to sgives ∇t, ~ XZ(s, ~ X(s)) 1 ~ U(s, ~ X(s)) =∂th(s, x(s))+∂xh(s, x(s))dx ds(s).(2.67) We compute dx ds =∂τx+ (∂Xx)dX ds , dX ds =∇t, ~ XX(s, ~ X(s)) 1 ~ U(s, ~ X(s)) . (2.68) Using (2.18) and (2.23), equation (2.67) yields b ∂tZ+W=∂th+−cos θ ∂tX+ cos θ(b ∂tX+U/J)∂xh. (2.69) With (2.28), (2.29) we get −cos θ ∂tb+W=∂th+cos θ J(Zcos2θ ∂2 txb−sin θ ∂tb+U)∂xh. (2.70) This must hold at any point on the free surface, therefore the kinematic condition is that for any (t, x), J|h cos θ∂th+ (hcos2θ ∂2 txb−sin θ ∂tb+U|h)∂xh=J|h cos θ(−cos θ ∂tb+W|h),(2.71) where |hmeans that we take the value at Z=h(t, x). Now, taking into account the boundary condition (2.13), which can be written as WZ=0 =−ρr ρm−1cos θ ∂tb, (2.72) 17 integrating of the incompressibility equation (2.58) for Zbetween 0 and h(t, x), and using the kinematic condition (2.71), gives ∂tZh 0 J cos θdZ +∂xZh 0 (Zcos2θ ∂2 txb−sin θ ∂tb+U)dZ =−ρr ρm ∂tb. (2.73) More explicitly, the mass conservation reads ∂th cos θ−h2 2∂xθ+∂xh2 2cos2θ ∂2 txb−hsin θ ∂tb+hU=−ρr ρm ∂tb, (2.74) where U=1 hZh 0 UdZ. (2.75) Notice that with known (2.58) and (2.72), the mass conservation (2.74) is equivalent to the kinematic condition (2.71). We can proceed similarly for the energy equation. Assuming that Pis isotropic, we integrate the energy equation (2.64) for Zbetween 0 and h(t, x). Using the kinematic condition (2.71) and the boundary conditions (2.72) and (2.10), we obtain the depth-integrated energy equation ∂tZh 0 J cos θU2+W2 2+g(b+Zcos θ)dZ +∂xZh 0(Zcos2θ ∂2 txb−sin θ ∂tb+U)U2+W2 2+g(b+Zcos θ)+P UdZ =− 1 2U2 Z=0 +1 2ρr ρm−1cos θ ∂tb2 +g b!ρr ρm ∂tb−ρr ρm−1PZ=0 ∂tb. (2.76) It is not possible to obtain an explicit formula for this integrated equation, unless a particular profile of the velocity in Zis assumed. The right-hand sides in (2.74) and (2.76) deserve some comments. For the mass conservation law (2.74), the term ρr ρm∂tbcan be put on the left-hand side. Then, multiplying the result by ρm, what one obtains is the time variation of the mass of moving fluid plus the mass of fluid at rest, which is the total mass of the fluid. For the energy equation (2.76), one can similarly put the term ρr ρmgb∂tb=∂t(ρr ρmgb2/2) on the left-hand side, since after multiplication by ρm, the term ρrgb2/2 represents the potential energy of the static part of the fluid. However, even when ρm=ρr, there remains the term −1 2U2 Z=0∂tbon the right-hand side, that cannot be put in conservative form. This term represents the variation of the kinetic energy attached to the displacement of the interface b(t, x). When ρm6=ρr, the kinetic energy at the interface is completed with the normal term W2 Z=0, which is non zero, and an additional term involving PZ=0 remains on the right-hand side, related to the mass exchanges through the interface. 18 2.8 Hydrostatic model Classically, the hydrostatic assumption consists in removing the normal acceleration in the momentum equation. Since Wis the normal velocity to the bottom interface, doing this means to remove the first line in (2.62). Therefore, the hydrostatic model, formulated with respect to the coordinates (t, x, Z), comprises equations (2.54), (2.60), and ∂xPXZ +J cos θ∂Zg(b+Zcos θ) + PZZ  =−U∂x(cos θ∂tb) + U∂xθ+ (PZZ −PXX )∂xθ. (2.77) The kinematic and boundary conditions remain unchanged for this system. If we assume that Pis isotropic, one can directly combine (2.54), (2.60) and (2.77) to obtain a hydrostatic energy equation. More explicitly, we can proceed as follows. We first transform the equation b ∂t(~g ·~ X) + ∇~ X·0 = 0 to the (t, x, Z) variables; this yields the identity ∂tJ cos θg(b+Zcos θ)+∂x(Zcos2θ ∂2 txb−sin θ ∂tb)g(b+Zcos θ) +∂Z−J∂tb g(b+Zcos θ)= 0, (2.78) that can also be checked directly. Then, we add up (2.60) times Uplus (2.77) times Wplus (2.54) times U2/2 + g(b+Zcos θ) + Pplus (2.57) times U2/2 plus (2.78). This gives the hydrostatic energy equation ∂tJ cos θU2 2+g(b+Zcos θ) +∂x(Zcos2θ ∂2 txb−sin θ ∂tb+U)U2 2+g(b+Zcos θ)+PU +∂ZJ cos θ(−cos θ ∂tb+W)U2 2+g(b+Zcos θ)+PJW cos θ= 0, (2.79) the only difference with (2.64) being that the term W2/2 has disappeared. Indeed, the same computation for the non-hydrostatic model gives (2.64). Therefore, for the hydrostatic model, the mass conservation still reads (2.74), while the integrated energy equation takes the form ∂tZh 0 J cos θU2 2+g(b+Zcos θ)dZ +∂xZh 0(Zcos2θ ∂2 txb−sin θ ∂tb+U)U2 2+g(b+Zcos θ)+P UdZ =−1 2U2 Z=0 +g bρr ρm ∂tb−ρr ρm−1PZ=0 ∂tb. (2.80) 19 3 A new erosion model for avalanches In this section we propose a new erosion model for avalanches, derived from our formulation (2.54), (2.60), (2.62) for the incompressible equations subject to the boundary conditions (2.72), (2.10) and the mass conservation (2.74). We first deal with the case where the stress tensor Pis isotropic, then in Section 3.3 non-diagonal terms are considered in order to include friction effects. Our procedure consists, as usual, in assuming for the velocity U(t, x, Z) a particular profile in Z. An approximate solution to the incompressible system is searched for with this particular profile, retaining only lower order terms relatively to the aspect ratio characterizing the shallowness of the moving mass. Indeed, we retain all first and second-order terms, and some third-order terms that are necessary to obtain a model that verifies exactly an energy inequality. Interestingly, the resulting model contains some first-order terms that do not appear in previously proposed models. 3.1 Formal derivation In what follows we consider a moving layer over a given variable bottom interface b(t, x), as described at the beginning of Section 2. Thus, according to (2.1), θis also known. We denote by εthe aspect ratio between the characteristic lengths normal and parallel to the interface, which is supposed to be a small parameter. Thus we have the shallowness assumption h=O(ε).(3.1) We look for a solution with a velocity profile almost constant in Z, U(t, x, Z) = u(t, x) + O(ε2).(3.2) Moreover, we assume a small variation in space of the angle θand a small variation in time of b, ∂xθ=O(), ∂tb=O().(3.3) Under these hypotheses, we make the following approximations. We assume that there are no singularities in time or space; i.e. ∂tand ∂xare formally bounded operators. On the other hand, the shallowness assumption leads to ∂Z=O(1/). At first, since 0 < Z < h, (3.1) ensures that Z=O(ε). Notice that (2.54), together with (2.72), can be used to define W(t, x, Z), once U(t, x, Z) is known. From this we deduce that W=O(ε), ∂tW=O(ε), ∂xW=O(ε). With (3.3) we have J= 1 + O(ε2), and (2.62), (3.3) give ∂ZgZ cos θ+P=O(ε),(3.4) implying together with (2.10) the near hydrostatic pressure P=gcos θ(h−Z) + O(ε2).(3.5) 20 Substituting this expression into (2.60), we obtain, on using (3.2), 1−Zcos θ ∂xθ cos θ∂tu+ (Zcos2θ ∂2 txb−sin θ ∂tb)∂xu+∂x(u2 2+gh cos θ+gb) = O(ε2). (3.6) We observe that the terms depending on Z, i.e. Zcos θ ∂xθand Zcos2θ ∂2 txb, are second-order terms, according to (3.3). Therefore, they can be neglected; thus we get an equation in the variables (t, x) only. This justifies the compatibility of the choice (3.2) with the Euler system. However, although the two previously mentioned terms depending on Zcan be neglected, we keep them because they are necessary in order to obtain a system which verifies exactly an energy inequality as will be discussed in the following section. Indeed, in order to remove the dependence in Z, we can take Z=h/2 in (3.6) without changing the order of approximation. We close the system by rewriting (2.74), where U=u+O(ε2). This gives an error in O(ε3) in (2.74). Finally, by expanding (2.76) we obtain the following equation for the energy, ∂th cos θ−h2 2∂xθu2 2+gb+gh2+ρr ρmb2 2−gh3 3cos θ ∂xθ +∂xu2 2+gh cos θ+gbhu +h2 2cos2θ ∂2 txb−hsin θ ∂tbu2 2+gb +gh3 3cos3θ ∂2 txb−gh2 2sin θcos θ ∂tb =−u2 2 ρr ρm ∂tb−ρr ρm−1gh cos θ ∂tb+O(ε3). (3.7) 3.2 Properties of the model In this subsection we study the properties of the model obtained as explained above, and dropping the error terms. The model is defined by the two mass and momentum equations coming from (2.74) and (3.6) (with Z=h/2), ∂th cos θ−h2 2∂xθ+∂xh2 2cos2θ ∂2 txb−hsin θ ∂tb+hu=−ρr ρm ∂tb, (3.8) 1 cos θ−h 2∂xθ∂tu+h 2cos2θ ∂2 txb−sin θ ∂tb∂xu+∂xu2 2+gh cos θ+gb= 0.(3.9) Notice that when ∂tb= 0, the system differs slightly from the model proposed in [5]. As usual, (3.9) is valid for smooth solutions, but for possibly discontinuous solutions, one should rather write the momentum equation in conservative form ∂t h cos θ−h2 2∂xθu+∂xh2 2cos2θ ∂2 txb−hsin θ ∂tb+huu+gh2 2cos θ =−gsin θh cos θ−h2 2∂xθ−ρr ρm u∂tb, (3.10) 21 obtained by multiplying (3.8) by u, (3.9) by hand adding the results. The term gh2 2sin θ ∂xθon the right-hand side of (3.10) was introduced by Bouchut et al. in [5]. Here it plays the same role, to maintain a conservative energy balance and to have the lake at rest solution. Theorem 3.1 System (3.8)-(3.9) has the following properties. (i) It admits an energy dissipation inequality ∂th cos θ−h2 2∂xθu2 2+gb+gh2+ρr ρmb2 2−gh3 3cos θ ∂xθ +∂xu2 2+gh cos θ+gbhu +h2 2cos2θ ∂2 txb−hsin θ ∂tbu2 2+gb +gh3 3cos3θ ∂2 txb−gh2 2sin θcos θ ∂tb≤ −u2 2 ρr ρm ∂tb−ρr ρm−1gh cos θ ∂tb (3.11) which is an identity for smooth solutions, (ii) it has the lake-at-rest solution u= 0,hcos θ+b=Cst, for arbitrary given b(t, x)if ρm=ρr, or for bsatisfying ∂tb≡0if ρm6=ρr, (iii) under assumptions (3.1), (3.3), it provides in the limit ε→0an approximation up to error terms of order ε3in (3.8) and of order ε2in (3.9), with constant velocity profile in the normal variable, to the free surface incompressible Euler equations (2.54), (2.60), (2.62), (2.72), (2.10), (2.74), (iv) for constant slope (∂tθ= 0,∂xθ= 0), system (3.8)-(3.9) gives an exact solution to the free surface incompressible Euler system with hydrostatic assumption (2.54), (2.60), (2.77), (2.72), (2.10), (2.74). Proof. For (i), a long calculation shows that for smooth solutions, by multiplying (3.8) by u2/2+gh cos θ+gb, (3.9) by hu and adding the results, using (2.1), (2.5), we obtain exactly the equality in (3.11). Since the details of the computation are very tedious, we prefer to omit them, and rather give an abstract argument. Consider a solution (h, u) to (3.8)-(3.9). Set U(t, x, Z) = u(t, x), define Wby (2.54) with the boundary condition (2.72), and take P=gcos θ(h−Z), that satisfies the boundary condition (2.10). Then, the mass conservation (2.74) holds, and only the velocity equations for Uand Wfail to hold in order to have a solution to the Euler equations. We have the equation ∂Z(g(b+Zcos θ)+P) = 0, thus we can say that the hydrostatic equation (2.77) holds modulo the term −U∂x(cos θ∂tb) + U∂xθ. For U, we can say that (2.60) holds modulo the term W∂x(cos θ∂tb) + U∂xθ, with a right-hand side (h/2−Z)∂xθ ∂tu−cos2θ ∂2 txb ∂xu,(3.12) 22 just because (3.9) has been obtained by taking Z=h/2 in (2.60). Performing the same computation as for the hydrostatic model, we can derive an energy equation. We observe that the terms involving ∂x(cos θ∂tb) + U∂xθsimplify, thus we get the energy equation (2.79) with the right-hand side utimes (3.12). Integrating this equation for Zbetween 0 and hthen gives (2.80), since the integral in Zof (3.12) vanishes. Finally the evaluation of (2.80) gives (3.11). The property (ii) can be obtained as follows. We have to prove that if u= 0 and hcos θ+b=Cindependent of xand t, then we have a solution to (3.8)- (3.9), provided that ( ρr ρm−1)∂tb≡0. The second equation (3.9) holds trivially, thus we just have to prove the first, (3.8). We first observe that ~ U= 0, P= gC −~g ·~ Xgives a solution to (2.6)-(2.7). Since the free surface is defined by h(t, x) such that hcos θ+b=C, it is horizontal. Thus (2.9) is satisfied. The boundary condition (2.13) holds according to the assumption, and the pressure Palso vanishes at the free surface, ensuring that (2.10) holds. Therefore, we have a solution to the incompressible Euler equations with boundary conditions. Formulating this solution in the (t, x, Z) coordinates reads in particular as equation (2.74), proving the result. Another proof is to differentiate the relation hcos θ+b=Cwith respect to tand x, and to put the results directly in (3.8). Notice that when ρm=ρrwith bdepending on time, the lake-at-rest solution is not strictly speaking a steady state, since halso depends on time. This is due to the fact that in this case, the interface is artificial and does not have a physical meaning, since the parts of the fluid above and below it are both at rest. The property (iii) follows directly from the derivation above. For (iv), consider a solution (h, u) to (3.8)-(3.9). As in the proof of (i), let U(t, x, Z) = u(t, x), define Wby (2.54) with the boundary condition (2.72), and take P=gcos θ(h−Z). According to the assumptions that ∂tθ= 0 and ∂xθ= 0 (but without assuming that ∂tb= 0), the terms ∂x(cos θ∂tb) + U∂xθand (3.12) vanish. Thus the computations made in the proof of (i) show that we have an exact solution to the hydrostatic Euler system.  3.3 Friction terms In this subsection we explain how from the incompressible system with non-isotropic stress tensor we can derive friction terms in the proposed model (3.8)-(3.9). In the system (2.54), (2.60), (2.62), the stress tensor Pneeds to be defined. Making the same scale analysis as previously performed for (2.62) to obtain (3.4), we obtain ∂ZgZ cos θ+PZZ=O(ε), provided that PXZ =O(ε). Thus, PZZ =gcos θ(h−Z) + O(ε2).(3.13) In order to specify PXX it is necessary to give a constitutive law. With the purpose 23 of simplicity in this work we set PXX =PZZ, although the case PXX =KPZZ can be considered, K=Kact/pass following the Savage-Hutter model [22]. The friction effects are introduced in the shear stress PZX. Following the works of Savage-Hutter [22] and Gray [12], and in line with (2.10), we impose as boundary condition PZX to vanish at the free surface, and a Coulomb friction term at the flowing edge of the interface, (PZX )|Z=h= 0,(PZX)|Z=0 =−µsgn(U)(PZZ)|Z=0 ,(3.14) where µ= tan θsand θsis the basal friction angle. Now, if we assume that µ=O(ε), we obtain PZX =O(ε2). Arguing as in the beginning of Section 3, the only new term appearing in (2.60) is ∂Z(JPZX /cos θ). Thus, averaging (2.60) between Z= 0 and Z=h, we are led to 1 cos θ−h 2∂xθ∂tu+h 2cos2θ ∂2 txb−sin θ ∂tb∂xu+∂xu2 2+gh cos θ+gb =PZX |Z=0 hcos θ+O(ε2). (3.15) However, according to (3.14) and (3.13), we have PZX|Z=0 =−µgh cos θsgn(u) + O(ε3). Therefore, neglecting terms in ε2in (3.15), we obtain together with (2.74) the system ∂th cos θ−h2 2∂xθ+∂xh2 2cos2θ ∂2 txb−hsin θ ∂tb+hu=−ρr ρm ∂tb, (3.16) 1 cos θ−h 2∂xθ∂tu+h 2cos2θ ∂2 txb−sin θ ∂tb∂xu +∂xu2 2+gh cos θ+gb=−gµ sgn(u). (3.17) A more precise argument is indeed to say that taking into account the conditions on PZX , a possible ansatz is PZX =−µsgn(u)gcos θ(h−Z) + O(ε3).(3.18) Then we obtain ∂ZJPZX cos θ=gµ sgn(u) + O(ε2).(3.19) Substituting in (2.60) gives (3.17), up to terms of order ε2. 3.4 Simplified model The model (3.8)-(3.9) is a bit complicated, thus it is tempting to try to simplify it, by dropping small terms. However, doing this, we break the nice properties stated in Theorem 3.1. However, it is worthwhile to state the result. The simplified system is defined as ∂th cos θ+∂x−hsin θ ∂tb+hu=−ρr ρm ∂tb, (3.20) 24 1 cos θ∂tu+−sin θ ∂tb∂xu+∂xu2 2+gh cos θ+gb= 0.(3.21) The corresponding conservative momentum equation is ∂th cos θu+∂x−hsin θ ∂tb+huu+gh2 2cos θ =−gsin θh cos θ−h2 2∂xθ−ρr ρm u∂tb, (3.22) obtained by multiplying (3.20) by u, (3.21) by hand adding the results. Theorem 3.2 The system (3.20)-(3.21) has the following properties. (i) It admits an energy dissipation inequality ∂th cos θu2 2+gb+gh2+ρr ρmb2 2 +∂xu2 2+gh cos θ+gbhu +−hsin θ ∂tbu2 2+gb −gh2 2sin θcos θ ∂tb≤ −u2 2 ρr ρm ∂tb−ρr ρm−1gh cos θ ∂tb +1 2gh2∂xθ ∂tb−∂xb ∂tθ (3.23) which is an identity for smooth solutions, (ii) it has the lake-at-rest solution u= 0,hcos θ+b=Cst, whenever ∂tb≡0, (iii) under assumptions (3.1), (3.3), it provides in the limit ε→0an approximation up to error terms of order ε3in (3.20) and of order ε2in (3.21), with constant velocity profile in the normal variable, to the free surface incompressible Euler equations (2.54), (2.60), (2.62), (2.72), (2.10), (2.74), (iv) for constant slope (∂tθ= 0,∂xθ= 0), system (3.20)-(3.21) gives an exact solution to the free surface incompressible Euler system with hydrostatic assumption (2.54), (2.60), (2.77), (2.72), (2.10), (2.74). Proof. Property (iii) follows from the fact that we only dropped terms of order ε3in (3.20) and ε2in (3.21) as compared to (3.8)-(3.9). Property (iv) is also obvious since there is no difference between (3.20)-(3.21) and (3.8)-(3.9) when ∂tθ= 0, ∂xθ= 0. For (ii), we observe that ∂tb= 0 implies that ∂tθ= 0. Thus, differentiating the identity hcos θ+b= Cst with respect to time yields ∂th= 0. We deduce obviously that (3.20)-(3.21) hold. 25 the velocity of grains in the upper layer (see [2], [7], [8]). Owing to equation (1.4), the BCRE model can be written as    ∂tb=−Vup(θ−θn), ∂th cos θ+∂x(hVd) = Vup(θ−θn),(4.16) where Vup(θ−θn) is the exchange term between the static and the mobile layers. By θnone denotes a neutral angle, which separates the erosion and deposition profile. The factor Vup has the dimension of a velocity and is a parameter of the model. For more details see [2], [7], [8]. In [2], Aradian et al. relate this model to the Saint Venant model to derive the parameters Vup and Vdinvolved in the BCRE model. Nevertheless, as usual, they neglect the term U|b∂tbin the momentum equation. We show here that including this term can change the determination of the parameter Vup by a factor 3/2. Let us start as in the derivation of the friction proposed in Subsection 3.3, but instead of a constant velocity profile, we rather take, as in [2], U=−Γ0Z, (4.17) where Γ0=pg/d and dis the grain diameter. Thus U < 0 which is coherent with the situation, when 0 ≤θ < π/2. Observe that (4.17) gives ∂tU= 0, ∂xU= 0. Then, since u≡U=−Γ0h/2, averaging (2.60) between Z= 0 and Z=hgives ∂xgh cos θ+gb+ Γ0∂tb=PZX|Z=0 hcos θ+ (1 −δ)u h∂tb+O(ε2),(4.18) with PZX |Z=0 =−µsgn(U)gh cos θ. (4.19) In the above, δis a parameter, formally inserted to identify (4.18) with (3.35). This has to be completed by the mass conservation equation obtained from (2.74), ∂th cos θ−h2 2∂xθ+∂xh2 2cos2θ ∂2 txb−hsin θ ∂tb−Γ0 h2 2=−∂tb. (4.20) Neglecting second-order terms gives ∂th cos θ+∂x−Γ0 h2 2=−∂tb+O(ε2).(4.21) Combining (4.21) and (4.18) where we neglect the term ∂x(gh cos θ), we obtain      ∂th cos θ+∂x−Γ0 h2 2=g Γ0cos θ(3 −δ)/2(sin θ−µcos θ), ∂tb=−g Γ0cos θ(3 −δ)/2(sin θ−µcos θ). (4.22) 32 If µ= tan θsthen sin θ−µcos θ=1 cos θs sin(θ−θs).(4.23) Then, if θ≈θsthe following system is obtained,      ∂th cos θ+∂x−Γ0 h2 2=g Γ0cos2θs(3 −δ)/2(θ−θs), ∂tb=−g Γ0cos2θs(3 −δ)/2(θ−θs). (4.24) If we compare (4.16) and (4.24), following [2], we obtain the identifications Vd=−Γ0 h 2, θn=θs,Vup =g Γ0cos2θs(3 −δ)/2.(4.25) We observe that since u=−Γ0h/2, we have Vd=u. Finally, we observe that Vup, the uphill velocity of exchange from the static to the mobile layers, is modified by the choice of δ. In fact, one has (Vup)δ=1 =3 2(Vup)δ=0. Acknowledgments The authors are very grateful to Kolumban Hutter for his detailed comments that helped a lot in improving the paper. In particular, the generalization to take into account the case when ρm6=ρr has been performed according to his suggestion. We thank the members of GDR MIDI for fruitful discussions. This work has been partially supported by the ACI NIM contract no 03318, the ANR contract no ANR-06-BLAN-0414, the ACI Jeunes Chercheurs ”G´eomorphog´en`ese, Erosion et Transport Granulaire”, the Spanish Government Research project MTM2006-08075. IPGP contribution number 2292. References [1] A. Aradian, Quelques probl`emes d’interface mole. Ph.D. Thesis Paris VI (2001). [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. 65, 061303 (2002). [4] A. Bon, Avalanche d’un milieu granulaire. Master 2, University Paris 6. [5] F. Bouchut, A. Mangeney-Castelnau, B. Perthame, J.-P. Vilotte, A new model of Saint Venant and Savage-Hutter type for gravity driven shallow water flows. C.R. Acad. Sci. Paris, Ser I 336, 531-536 (2003). 33 [6] F. Bouchut, M. Westdickenberg, Gravity driven shallow water model for arbitrary topography. Comm. Math. Sci. 2(3), 359-389 (2004). [7] J.P. Bouchaud, M.E. Cates, J.R. Prakash, S.F. Edwards. A model for the dynamics of sandpile surface. J. Phys. Paris I 4, 1383-1410 (1994). [8] T. Boutreux, E. Raphael, P.G. DeGennes, Surface flows of granular materials: a modified picture for thick avalanches. Phys. Rev. E 58, 4692-4700 (1998). [9] S. Douady, B. Andreotti, A. Daerr, On granular surface flow equations. Eur. Phys. J.B. 11, 131-142 (1999). [10] M. Papa, S. Egashira, T. Itoh, Critical conditions of bed sediment entrainment due to drebis flow, Natural Hazards and Earth System Sciences 4, 469-474 (2004). [11] J.-F. Gerbeau, B. Perthame, Derivation of viscous Saint Venant system for laminar shallow water; numerical validation. DCDS(B) 1 (1) (2001). [12] J.M.N.T. Gray, Granular flow in partially filled slowly rotating drums. J. Fluid Mech 441, 1-29 (2001). [13] 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). [14] D.V. Khakhar, A.V. Orpe, Surface granular flows: two related examples. Advances in Complex Systems 4(4), 407-417 (2001). [15] K.K.J. Kouakou, Instabilit´e des interfaces ´ecoulement/fond ´erodable : approches asymptotiques. Preprint University of Paris VI (2005). [16] L.D. Landau, and Lifshitz, E.M., Statistical Physics, Pargamon Press, New- York (1980) [17] A. Mangeney-Castelnau, J.-P. Vilotte, M.-O. Bristeau, B. Perthame, F. Bouchut, C. Simeoni, S. Yerneni, Numerical modeling of avalanches based on Saint Venant equations using a kinetic scheme. J. Geophys. Res. 108 (B11), 2527 (2003). [18] A. Mangeney-Castelnau, F. Bouchut, J.-P. Vilotte, E. Lajeunesse, A. Aubertin, M. Pirulli, On the use of Saint Venant equations to simulate the spreading of a granular mass, J. Geophys. Res. 110 (B9), B09103 (2005). [19] A. Mangeney, F. Bouchut, N. Thomas, J.-P. Vilotte, M.-O. Bristeau, Numerical modeling of self-channeling granular flows and of their levee/channel deposits, to appear in J. Geophys. Res. (2007). 34 [20] A. Mangeney, L.S. Tsimring, D. Volfson, I.S. Aranson, F. Bouchut, Avalanche mobility induced by the presence of an erodible bed and associated entrainment, Geophys. Res. Lett., in press (2007). [21] Pudasaini, S., and Hutter, K., Avalanche Dynamics, Springer-Verlag, 602p (2007). [22] S.B. Savage, K. Hutter, The dynamics of avalanches of granular materials from initiation to run-out, Acta Mech. 86, 201-223 (1991). [23] B. Sovilla, P. Burlando, and P. Bartelt, Field experiments and numerical modeling of mass entrainment in snow avalanches J. Geophys. Res., 111, F03007 (2006). [24] 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). 35