A two-phase shallow debris flow model with energy balance

ESAIM: M2AN 49 (2015) 101–140 ESAIM: Mathematical Modelling and Numerical Analysis DOI: 10.1051/m2an/2014026 www.esaim-m2an.org A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE F. Bouchut1,E.D.Fern ´ andez-Nieto2, A. Mangeney3,4 and G. Narbona-Reina5 Abstract. This paper proposes a thin layer depth-averaged two-phase model provided by a dissipative energy balance to describe avalanches of solid-fluid mixtures. This model is derived from a 3D two-phase model based on the equations proposed by Jackson [The Dynamics of Fluidized Particles. Cambridges Monographs on Mechanics (2000)] which takes into account the force of buoyancy and the forces of interaction between the solid and fluid phases. Jackson’s model is based on mass and momentum conservation within the two phases, i.e. two vector and two scalar equations. This system has five unknowns: the solid volume fraction, the solid and fluid pressures and the solid and fluid velocities, i.e. three scalars and two vectors. As a result, an additional equation is necessary to close the system. Surprisingly, this issue is inadequately accounted for in the models that have been developed on the basis of Jackson’s work. In particular, Pitman and Le [Philos.Trans.R.Soc.A363 (2005) 799–819] replaced this closure simply by imposing an extra boundary condition. If the pressure is assumed to be hydrostatic, this condition can be considered as a closure condition. However, the corresponding model cannot account for a dissipative energy balance. We propose here a closure equation to complete Jackson’s model, imposing incompressibility of the solid phase. We prove that the resulting whole 3D model is compatible with a dissipative energy balance. From this model, we deduce a 2D depth-averaged model and we also prove that the energy balance associated with this model is dissipative. Finally, we propose a numerical scheme to approximate the depth-averaged model. We present several numerical tests for the 1D case that are compared to the results of the model proposed by Pitman and Le. Mathematics Subject Classification. 65C20, 81T80, 91B74, 97M10. Received April 9, 2013. Revised May 5, 2014. Published online January 14, 2015. Keywords and phrases. Granular flows, two-phase flows, thin layer approximation, energy balance, non-conservative systems, projection method, finite volume schemes. 1Universit´e Paris-Est, Laboratoire d’Analyse et de Math´ematiques Appliqu´ees (UMR 8050), CNRS, UPEMLV, UPEC, 77454 Marne-la-Vall´ee, France. [email protected] 2Departamento de Matem´atica Aplicada I, Universidad de Sevilla. E.T.S. Arquitectura. Avda, Reina Mercedes, s/n. 41012 Sevilla, Spain. [email protected] 3Universit´e Paris Diderot, Sorbone Paris Cit´e, Institut de Physique du Globe de Paris, Seismology group, 1 rue Jussieu, 75005 Paris, France. [email protected] 4ANGE group INRIA, Jacques Louis Lions, CETMEF 5Departamento de Matem´atica Aplicada I, Universidad de Sevilla. E.T.S. Arquitectura. Avda, Reina Mercedes, s/n. 41012 Sevilla, Spain. [email protected] Article published by EDP Sciences c EDP Sciences, SMAI 2015 102 F. BOUCHUT ET AL. (b) (a) Figure 1. (a) Deposits of several debris flows in Iceland. (b) Close-up of a cross-section of the deposit of a debris flow covering a road in Canada. 1. Introduction Landslides, debris avalanches or debris flows play a key role in erosion processes on the surface of the Earth and other telluric planets. On Earth, they represent one of the major natural hazards. Gravitational instabilities are also closely related to volcanic, seismic and climatic activity and thus represent potential precursors or proxies for the change of these activities with time. Research involving the dynamic analysis of gravitational mass flows is advancing rapidly. One of its ultimate goals is to produce tools for detection of natural instabilities and for prediction of velocity and runout extent of rapid landslides. The theoretical description and physical understanding of these processes in a natural environment are still open and extremely challenging problems for earth scientists, giving rise to equally challenging mechanical, mathematical and numerical issues. In recent years, significant progress in the mathematical, physical and numerical modelling of gravitational flows has made it possible to develop and use numerical models to investigate geomorphological processes and assess risks related to such natural hazards. However, key questions still remain unanswered, for instance concerning the reason for the high mobility of natural landslides (e.g. [21,22]). Severe limitations prevent a full understanding of physical processes involved in landslide dynamics and the development of tools for detection of instabilities and prediction of their velocity and extent. Indeed, numerical models do not take into account complex natural phenomena such as the static/flowing transition in granular flows or the co-existence and interaction of fluid (water, gas) (e.g. [7,8,14,18,23,26–28,42,43]). Water is almost always involved in natural landslides (e.g. [15,16,29]) (Fig. 1b). Interaction forces between the solid and fluid (water) phases may play an important role in flow mobility and deposit extent. Different approaches can be used to simulate fluid-solid mixtures, extending from discrete element models based for example on contact dynamics or molecular dynamics (e.g. [32,47]), and taking into account individual particles, to continuum models that deal with a fluid phase and a solid phase. The discrete element approach is hard to use in geophysical applications due to the high computational costs required to take into account the broad-size distribution of particles in real flows, which is critical in such simulations. Existing models used to describe the behaviour of fluid-solid mixtures are mainly based on Jackson’s model [19]. This model takes into account solid and fluid stresses, the interaction force between the fluid and solid phases and the buoyancy force, through mass and momentum conservation within the two phases. This model thus involves four equations (two scalar and two vector equations). However, the system has five unknowns: the solid volume fraction, the solid and fluid pressures and the solid and fluid velocities (three scalars A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 103 and two vectors). As a result, an additional equation is necessary to close the system. Surprisingly, this issue is inadequately accounted for in the models that have been developed on the basis of Jackson’s work. Solving the 3D two-phase equations leads to high computational costs. For this reason, mostly depth-averaged models have been proposed to deal with natural geophysical flows (e.g. [10,35,36,38]). Iverson [15]wasthefirstto address the need to include interstitial fluid effects in the constitutive behaviour of the mass flow and developed a thin layer model for a solid-fluid mixture moving on realistic terrain, under the simplifying assumptions of constant porosity and equality of the fluid and solid velocity. The flow is described by a single set of equations for the density and momentum of the mixture, which is formally represented by a single-phase model with a stress term accounting for contributions from the two constituents. Due to the lack of an explicit equation for the pore fluid pressure in this model, a pore pressure advection-diffusion equation was added based on experimental measurements. Various versions and applications of this grain-fluid mixture model have since been presented (e.g. Pudasaini et al. [39]; Georges and Iverson, [10]). Taking another step forward, Pitman and Le proposed in [38] a novel depth-averaged two-fluid model for debris flows, based on Jackson’s model, that contains mass and momentum equations for both the fluid and solid phases, thus providing equations for the velocities of the two phases and for solid volume fraction. In the model proposed by Pitman and Le and the modified version proposed by Pelanti et al. [36], the authors do not provide a closure equation for the two-phase model. On the other hand, they impose two boundary conditions involving vanishing surface tension conditions at the free surface, i.e. the pressure of both the solid and the fluid phases vanish at the free surface. Two kinematic boundary conditions are also imposed at the free surface, because the two phases are assumed to fill a common domain, this gives an overdetermined problem at the free surface. However, in the thin layer approximation, because of the hydrostatic pressure assumption, the extra boundary condition makes it possible to express a depth-averaged model, even though no closure relation for the whole system is provided. However, boundary conditions obviously do not replace a closure equation inside the domain. This artificial compensation of the missing closure equation by overdetermined boundary conditions leads to a physically irrelevant energy equation in the Pitman−Le model (see Sect. 4.1). A physically meaningful energy equation is essential to obtain realistic models. A key issue in two-phase flow models is thus to propose a suitable closure relation that is compatible with the energy balance. Some new and very useful ways to close the system of equations have been proposed by Roux and Radjai [44], Pailha and Pouliquen [35] and George and Iverson [10]. The general idea is to take into account the dilation/compression of the granular phase and its interaction with the pressure of the fluid filling the pores of the granular material. Indeed, these effects have been shown to be crucial at the initiation of mass destabilization and to have a strong impact on the generated flow dynamics [e.g. [17,40]]. Roux and Radjai [44] proposed an equation to describe the evolution of the volume fraction and of the shear stress in a granular material in terms of the shear-induced dilatancy (a property of the granular material related to its dilation when the material is submitted to a shear force). Pailha and Pouliquen used this equation to close their model, based on the two-phase approach proposed by Jackson [35]. They also imposed that both the solid and fluid pressures vanish at the free surface. Moreover, they introduced a closure equation, related to dilatancy effects (Eq. (3.18) in [35]). But, as the resulting system is overdetermined, a condition had to be relaxed. Indeed, they relaxed mass conservation for one of the two phases, that they justified by assuming that the thickness of the flowing mixture is nearly constant. Alternatively, George and Iverson [10] derived a model using the mass and momentum equations of the mixture. In their model, the unknowns are the total height and velocity of the mixture, the solid volume fraction and the pore fluid pressure. As a closure relation, they used a slightly different equation than that proposed in [44] to describe dilatancy effects, that includes the time derivative of the effective normal stress and the pore fluid pressure. This relation is derived from the mass conservation of the solid phase by assuming that the averaged mixture velocity is equal to the averaged solid velocity (Eqs. (6) and (7) of their paper) and a Darcy law. However, the final model does not impose explicitly the mass conservation of the solid phase. We propose here to solve the mass and momentum equations of both phases, together with the relevant number of boundary conditions and a closure equation that provides a possibly physically relevant energy 104 F. BOUCHUT ET AL. equation. In a first step toward this objective, we use the simplest closure equation (i.e. incompressibility of the solid phase). We impose a vanishing stress condition at the free surface for the mixture (not for each phase) and kinematic surface boundary conditions (the two phases are supposed to fill the same domain), forming a well-posed 3D system. The analysis of the hydrostatic approximation suggests that a variable related to the pressure field remains in the thin-layer asymptotics. On choosing a static constraint as a closure relation, this extra variable can be determined as the associated Lagrange multiplier. The resulting model has a built-in energy balance. 2.The3Dtwo-phasemodel In this section we present the three-dimensional model used to describe the mixture of solid and fluid materials. Note that we do not consider here the role of the air (i.e. a third phase) that can be critical in some cases due to capillary forces, especially at the laboratory scale [16]. As a result, these equations are only valid when the granular media is saturated with fluid so that there is no air within the pores of the granular material. In Section 2.1, the mass and momentum equations of Jackson’s model are presented and a closure equation is proposed. In Section 2.2, the boundary conditions are described. In Sections 2.3 and 2.4, we express the drag force and the assumptions concerning the stress tensor. Finally, in Section 2.5,weexpressthecompletemodel in local coordinates. 2.1. Mass and momentum equations We consider geophysical mass flows made of a mixture of solid and fluid materials. The two fluid model presented below is derived in the Jackson’s book [19]. It is based on the dynamics of an assembly of solid particles immersed in a Newtonian fluid. The two-fluid model is obtained by averaging in the whole region the fundamental equations for both components, the fluid and particles. Namely, the Navier-Stokes equation for the motion of the fluid and the equations of linear and angular momentum for each particle for the solid part. These equations are coupled by the no-slip boundary condition imposed on the surface of each particle. The mass and momentum conservation equations for fluid and particle phases are deduced by an averaging procedure. But, some terms linked to the microscopic level of the individual particles are neglected. Consequentely, after the averaging procedure there are more unknowns than equations in the derived system. And a closure for the system must be set. The two-phase model is defined by the following mass and momentum equations for the solid and fluid phases: ∂t(ρsϕ)+∇·(ρsϕv)=0,(2.1a) ∂t(ρf(1 −ϕ)) + ∇·(ρf(1 −ϕ)u)=0,(2.1b) ρsϕ(∂tv+(v·∇)v)=−∇ · Ts+f0+ρsϕg,(2.2a) ρf(1 −ϕ)(∂tu+(u·∇)u)=−∇ · Tf−f0+ρf(1 −ϕ)g,(2.2b) where the subscript “s” refers to the solid phase and the subscript “f” refers to the fluid phase. The velocities are vfor the solid phase and ufor the fluid phase. Tdenotes the stress tensor and ρthe density. Acceleration due to gravity is denoted by gand f0represents the average value of the resultant force exerted by the fluid on a solid particle. The solid volume fraction is ϕ. For monodisperse beads, the maximal volume fraction is ϕmax ≃0.6, while it can be higher than 0.9 for highly polydisperse materials because the small particles can fill the pore space between larger particles ([3,13,48]). The solid fraction is practically never equal to 1. The case of dry granular flows can be obtained by setting all the variables related to the fluid phase (fluid stress and f0) A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 105 to zero and volume fraction to one in equations (2.1a)and(2.2a). The minus sign on the stress tensor terms agrees with the sign convention used in soil mechanics, where stress is defined as positive in compression. Note that both the grain density ρsand the fluid density ρfare constant, so that each material is incompressible. However, the density of the solid phase ϕρs(i.e. density of the total amount of grains per unit volume) and the density of the fluid phase (1 −ϕ)ρf(density of the total amount of fluid filling the pores of the granular assembly per unit volume) can change because ϕvaries with space and time. In this sense, the solid and fluid phase could be compressible. Note that the combination of (2.1a)and(2.1b) defines mass conservation for the mixture: ∂t(ρm)+∇·(ρmvm)=0,(2.3) where ρm=ρsϕ+ρf(1 −ϕ)andvm=ρsϕv +ρf(1 −ϕ)u ρsϕ+ρf(1 −ϕ) are respectively the density and velocity of the mixture. Multiplying (2.1a)byρfand (2.1b)byρsgives: ∇·(ϕv +(1−ϕ)u)=0.(2.4) This relation is different from the one expressing incompressibility of the mixture because it does not imply that ∇·vmis equal to zero. The averaged value of the interaction force between fluid and particle is collected in f0. This force is decomposed into the sum of the buoyancy force fBand all remaining contribution f. The main components of this force fare a term depending on the particle concentration and the relative velocity (u−v)−the drag force −, a term depending on the concentration and the relative acceleration −the virtual mass force −and the third contribution due to the force normal to the direction (u−v)−the lift force −. Several expressions for the buoyancy force are discussed in [19]. In the simplest case this force is written as fB=−ϕ∇pfwith pfthe fluid pressure that resume the force exerted by the fluid at rest on an immersed body that is also at rest. The generalization to more general motions leads us to the expression fB=−ϕ∇Tf, however the approximation of −ϕ∇pfgives equivalent equations assuming that the action on the particles due to the gradient of the deviatoric part of Tfis collected by the other terms of the total force f.Sowewrite f0=fB+f=−ϕ∇pf+f. (2.5) In the case when the inertia associated to the relative motion of fluid and particles can be neglected, referred to as the short relaxation time approximation, the virtual mass force may be neglected compared to the drag force in the equations of motion. Regarding to the lift force, the algebraic expression for this contribution normal to the relative velocity is uncertain because it takes quite different forms for different flow regimes. This force will not be considered in this paper. Thus, we assume that fcan be expressed simply by the drag force. The drag force acts in the direction of the relative velocity (u−v) and also depends on the particle concentration. So in general it can by written as: β(ϕ, |u−v|)(u−v). For small values of |u−v|, this force is proportional to the relative velocity so we can write f=˜ β(u−v) (2.6) ˜ βbeing the drag coefficient (see Sect. 2.3). The notation ˜ βis used to distinguish this coefficient from the drag coefficients denoted by βin some other publications (see for example [35]). The effective stress tensors Tsand Tfare related to the interactions between fluid and particles and between particles themselves. In [19] is pointed out the difficulty of writing them in terms of averaged variables in order to close the system. Furthermore an extensive discussion is included to give some explicit and empirical closures for different regimes in terms of the Stokes number. In this work, we take a symmetric solid stress tensor Ts and denote psits first invariant, i.e., the pressure of the solid phase (see Sect. 2.4). 106 F. BOUCHUT ET AL. The viscosity of the fluid acts at the “macroscopic” scale through viscous terms of order μU/L2,whereU and Lare characteristic values of respectively the fluid velocity and flow length. On the other hand, the fluid viscosity acts at the “microscopic” scale during the relative motion between the fluid phase and the granular porous media commonly described by the Darcy law. This microscopic contribution is of the order of μΔU/κ, where κis the intrinsic hydraulic permeability of the granular media and ΔU is the typical relative velocity of the fluid phase with respect to the solid phase. Here we assume that the “macroscopic” viscous forces related to the fluid are negligible, so that the fluid stress tensor reduces to the pressure term, ∇·Tf=∇pf.(2.7) By substituting these expressions into (2.2a)and(2.2b), we obtain the system (2.1a), (2.1b), and ρsϕ(∂tv+(v·∇)v)=−∇ · Ts−ϕ∇pf+f+ρsϕg,(2.8a) ρf(1 −ϕ)(∂tu+(u·∇)u)=−(1 −ϕ)∇pf−f+ρf(1 −ϕ)g.(2.8b) This system of equations is the same as the system considered in [16,38]. Only the boundary conditions are different from those used here. As discussed above, this system of four equations (2.1a), (2.1b), (2.8a), (2.8b) has five unknowns ϕ,Ts,pf, uand v. To close the system, we propose to add a supplementary scalar equation, based on the physical processes involved. Starting from the simplest closure relation, we propose to impose the incompressibility of the solid phase: ∇·v=0.(2.9) In real granular materials the dilatancy effect may induce changes of the volume of the solid phase, even if the mass of the granular material remains constant. This means that the divergence of the velocity of the solid phase vmay not be zero (see [11]). The compression/dilation of the granular phase changes the interstitial fluid pressure that in turn couples with the solid momentum equations. This coupling appears in the non-hydrostatic pressure terms (see [30]), not included in the approximations made in this work. The consistency of the whole model can be evaluated by the local energy balance equation. To obtain it, we multiply (2.8a), (2.8b)byvand urespectively, combine with (2.1a)and(2.1b), and add the results. This yields ∂tρsϕ|v|2 2+ρf(1 −ϕ)|u|2 2+∇·ρsϕ|v|2 2v+ρf(1 −ϕ)|u|2 2u =−v·(∇·Ts)−ϕv +(1−ϕ)u·∇pf+f·(v−u)+ρsϕv +ρf(1 −ϕ)u·g. (2.10) Denoting Xthe space position and once again using (2.1a)and(2.1b)alongwith(2.4), we obtain ∂tρsϕ|v|2 2+ρf(1 −ϕ)|u|2 2−(g·X)ρsϕ+ρf(1 −ϕ) +∇·ρsϕ|v|2 2v+ρf(1 −ϕ)|u|2 2u−(g·X)ρsϕv +ρf(1 −ϕ)u +pfϕv +(1−ϕ)u+Tsv =(Ts−psId) : ∇v+ps∇·v+f·(v−u),(2.11) where psdenotes the solid particles pressure. From equation (2.6), the drag contribution f·(v−u) is non-positive. With the assumption that the solid phase is incompressible, the second term on the right-hand side ps∇·vis equal to zero and it is natural to assume that the friction dissipation (Ts−psId) : ∇vis non-positive. As a result, the sum of the three terms in the right-hand side of (2.11) is non-positive. A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 107 The model defined by (2.1a), (2.1b), (2.8a), (2.8b) with closure (2.9) has a locally dissipative energy balance (2.11). Note that in the initial system considered by Pitman and Le, the term ps∇·vdoes not vanish and we cannot ensure the non-positiveness of the right-hand side term in (2.11). We will show in Section 4.1 that the term resulting from the closure equation also makes it possible to obtain a dissipative energy balance in the model. 2.2. Boundary conditions 2.2.1. At the free surface We consider the usual geometric setting, which is that the mixture lies in a spatial domain limited by a fixed topography at the bottom and by a free surface at the top. We assume that the fluid and the solid fill the same domain that is moving with the velocity of both. This gives the simultaneous kinematic conditions (1,u)·N=0,(1,v)·N=0 atthefreesurface,(2.12) where N=(Nt,N X) is the time-space normal. It can be rewritten u·NX=v·NX=−Ntat the free surface.(2.13) Note that this is a strong assumption that plays a key role in the derivation of the equations and in the resulting model presented below. In [16,35,38], both the fluid and the solid pressures are set to zero at the free surface. However, as discussed in the introduction, only one dynamic boundary condition can be imposed at thefreesurfaceofthemixture: (Ts+pfId)NX=0 atthefreesurface.(2.14) Remark 2.1 (about the total stress tensor).To obtain the total stress for the mixture we can combine equations (2.8a)and(2.8b). From here the total stress for the mixture becomes more complicated than the sum of the two stress tensors for each phase. Namely, it can be written as T=Ts+Tf+T=Ts+pfId +T,withT a contribution coming from the non-linear convective terms written through the relative velocities of the solid and the fluid with respect to the velocity of the mixture vm: T=−ρsϕ(v−vm)(v−vm)−ρf(1 −ϕ)(u−vm)(u−vm). Nevertheless, for many geophysical flows one can assume that this term is negligible, by assuming in particular that the relative velocity of the fluid with respect to the solid is small compared to the solid velocity (see p. 540 of [16] for details.). Thus, we can see condition (2.14) as a simplification where the total stress of the mixture at the free surface for our system is defined as the sum of the fluid and solid phase stress tensors. 2.2.2. At the bottom The conditions at the bottom are classically the non-penetration conditions u·n=0,v·n= 0 at the bottom,(2.15) where nis the upward space unit normal (i.e. the normal to the topography). This must be completed by further conditions for the solid, in particular we consider a Coulomb friction law, following [45] Tsn−(Tsn)·nn=−tan δsign(v)(Tsn)·nat the bottom,(2.16) where δis the intergranular Coulomb friction angle. 108 F. BOUCHUT ET AL. Remark 2.2. The system (2.1a), (2.1b), (2.8a), (2.8b), (2.9) with the boundary conditions (2.13), (2.14), (2.15), (2.16), is formally well-posed. Moreover, we can check that the previous boundary conditions ensure that all the boundary contributions vanish in the energy balance of the model, except the one coming from the Coulomb condition (2.16), which dissipates at the bottom. The main difference between this system and those considered by Pitman and Le (see [38]) and Pailha and Pouliquen (see [35]) is the definition of the boundary conditions. Instead of considering that the total pressure vanishes at the free surface (Eq. (2.14)), they consider that both the pressure of the solid phase and the pressure of the fluid phase vanish at the free surface. Pitman and Le do not consider any closure equation, consequently we cannot check the well-posedness of this system. Pailha and Pouliquen consider a closure equation in terms of the divergence of the solid phase velocity. Nevertheless, given that the system is overdetermined in this case and they relax the mass conservation of one of the two phases. 2.3. Assumptions concerning the drag force Different empirical relations are proposed in the literature for the drag force. As already mentioned, the drag force expression is assumed to be f=˜ β(u−v).(2.17) The drag coefficient ˜ βcan be defined in different ways: •Pitman and Le [38] used the drag force proposed by Richardson and Zaki (see [41]): ˜ β=(ρs−ρf)ϕg vT(1 −ϕ)m−1,(2.18) where vTis the terminal velocity of an isolated representative solid particle falling in the fluid under gravity. This force has been calculated by Richardson and Zaki, based on laboratory experiments measuring vTand vS,wherevSis the sedimentation velocity of the dispersion of particles in a fluid. Experiments give the empirical law: vS=(1−ϕ)nvT. The value of the empirical exponent nlies in the range [2.4,4.65]. Pitman and Le [38] (Appendix A) show that m=n−2, so that m∈[0.4,2.65]. Depending on the respective roles of viscous and inertial forces, vS/vTdepends or does not depend on the Reynolds number (see [41] for more details). For example, with the typical values of the experiments done by Iverson [17], the typical Reynolds number is Re=vTdρ f/μ≈50. From Table VI of Richardson and Zacki, this gives n≈3andthenm≈1. •Pailha and Pouliquen [35] use the following definition of the drag coefficient: ˜ β=(1−ϕ)2μ αd2,(2.19) μbeing the dynamic viscosity, dthe mean grain diameter and α=(1 −ϕ)3 150 ϕ2· This is derived from the Carman–Kozeny relation for the permeability of the porous media formed by the particles (see [12,33]). Another way to estimate ˜ βis to assume that the friction between the two phases is similar to the Darcy law. In debris flows, part of the vertical displacements and of the fluctuations of the horizontal displacement are induced by the dilation or the compaction of the granular media. These effects impact on the fluid pressure field that in turn affects the momentum conservation of the solid and fluid phases. The coupling between the fluid and solid phases comes from the drag force (see [35]). This can be understood by considering the deviation A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 109 from the hydrostatic fluid pressure. Let us denote pf=ph f+pe f,whereph fcorresponds to the hydrostatic fluid pressure, satisfying ∂zph f=−ρfgcos θ,andpe fis the excess pore-fluid pressure. If the right hand side of (2.8b) is considered predominant (small inertia), the horizontal variation of ph fis negligible and the gradient of the excess pore-fluid pressure pe fsatisfies ∇pe f=−˜ β (1 −ϕ)(u−v).(2.20) This formula has the same structure as the linear Darcian drag formula describing fluid flow within porous media. This law, considered in [10], relates u−vto the gradient of the excess pore-fluid pressure, ∇pe f=−μ κ(1 −ϕ)(u−v),(2.21) where μis the pore-fluid viscosity and κis the intrinsic hydraulic permeability of the granular debris. George and Iverson [10] point out that even if this linear drag formula may oversimplify the effects of complex phaseinteraction forces in debris flows, several research papers, such as [20,46], indicate that it probably provides a suitable first approximation. Comparing (2.20) with the Darcian law (2.21)leadsto: ˜ β=(1−ϕ)2μ κ,(2.22) κbeing the permeability of the granular media. The value of the effective permeability derived from (2.18) and (2.19), when compared to (2.22) gives respectively: •For Pitman and Le (2.18): κ=μvT(1 −ϕ)m+1 (ρs−ρf)gϕ ·(2.23) •For Pailha and Pouliquen (2.19): κ=d2(1 −ϕ)3 150ϕ2·(2.24) These two different values of permeability derived from (2.23)and(2.24) are compared in Figure 2.Forthis comparison we set μ=0.001 Pa s, d=10 −3m, (ρs−ρf) = 1500 kg m−3,g=9.81 m s−2,m=1,vT= 0.143 m s−1. We observe that both models give relatively close approximations of the permeability for values of ϕgreater than 0.4. George and Iverson [10] have simulated the experiments performed in [17]. In these experiments, the value of κwas approximately 10−12 m2, whereas the numerical simulations where performed with a constant value of κ≈10−8m2. Note that, with the definition of κgiven by (2.23)or(2.24), κ≈10−8m2when ϕ≈0.5. A value of κ≈10−12 m2corresponds to ϕ≈0.9. 2.4. Assumptions concerning the stress tensor To obtain the final system, a constitutive relation should be stated for the fluid and granular phases. •Fluid stress tensor Tf. As mentioned before, we assume that the fluid stress tensor can be expressed by the fluid pressure: Txy f=Txz f=Tyz f=0,T xx f=Tyy f=Tzz f=pf.(2.25) •Solid stress tensor Ts. We assume that all its components are proportional to the normal stress perpendicular to the topography, i.e. the stress component Tzz s, Tjk s=αjkTzz s,j,k=x, y, z. (2.26) 116 F. BOUCHUT ET AL. In this section, we will first establish a local energy equation for this model and then describe some of its properties. 4.1. Local energy In the following lines we prove that the model (4.1) is compatible with a dissipative energy balance. First, from the mass equations (4.1a), (4.1b)wehave ∂th+div(hϕv +h(1 −ϕ)u)=0,(4.2) ∂t(h(ρsϕ+ρf(1 −ϕ))) + div(h(ρsϕv +ρf(1 −ϕ)u)) = 0.(4.3) We can write sin θ(1,0)t=cosθ∇˜ bwith ˜ b=xtan θ, so that the sin θterms in (4.1c)and(4.1d)canbe grouped with the ∇bterms to give ∇(b+˜ b). Then we multiply equation (4.1c)by(hv)and(4.1d)by(hu) and sum up the results. Using the mass equations to simplify the left-hand side, we obtain ∂tρsϕh|v|2 2+ρf(1 −ϕ)h|u|2 2+divρsϕh|v|2 2v+ρf(1 −ϕ)h|u|2 2u =−(1 −ϕ)h(u−v)·∇(pfbed −ρfghcosθ)   (a) −ghcosθρsϕv +ρf(1 −ϕ)u·∇(b+˜ b+h)   (b) −h2 2gcos θ(ρs−ρf)v·∇ϕ   (c) −˜ βh|u−v|2−|v|tan δϕ(ρs−ρf)ghcosθ. Our objective is to compute each term on the right-hand side of the previous equation and try to write it as a time derivative or a divergence of something explicit. •Term (a). Using (4.1e), (a)=−div (1 −ϕ)h(u−v)(pfbed −ρfghcosθ).(4.4) •Term (b). Taking into account (4.3), (b)=−ghcosθ∇(b+˜ b)·(ρsϕv +ρf(1 −ϕ)u)−gcos θ∇h2 2·(ρsϕv +ρf(1 −ϕ)u) =−div gh(b+˜ b)cosθ(ρsϕv +ρf(1 −ϕ)u) +g(b+˜ b)cosθdiv h(ρsϕv +ρf(1 −ϕ)u) −div 1 2gh2cos θ(ρsϕv +ρf(1 −ϕ)u)+1 2gh2cos θdiv(ρsϕv +ρf(1 −ϕ)u) =−div ghcosθ(b+˜ b+h 2)(ρsϕv +ρf(1 −ϕ)u) −∂tgh(b+˜ b)cosθ(ρsϕ+ρf(1 −ϕ))+1 2gh2cos θdiv(ρsϕv +ρf(1 −ϕ)u). A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 117 •Term (c). (c)=−div 1 2gh2cos θ(ρs−ρf)ϕv+ϕ(ρs−ρf)gcos θdiv h2 2v. Gathering all the terms we get ∂tρsϕh|v|2 2+ρf(1 −ϕ)h|u|2 2+gh(b+˜ b)cosθ(ρsϕ+ρf(1 −ϕ)) +divρsϕh|v|2 2v+ρf(1 −ϕ)h|u|2 2u+(1−ϕ)h(u−v)(pfbed −ρfghcosθ) +ghcosθ(b+˜ b+h 2)(ρsϕv +ρf(1 −ϕ)u)+1 2gh2cos θ(ρs−ρf)ϕv =T1,(4.5) where T1canbeexpressedas T1=1 2gh2cos θdiv(ρsϕv +ρf(1 −ϕ)u)+ϕ(ρs−ρf)gcos θdiv h2 2v −˜ βh|u−v|2−|v|tan δϕ(ρs−ρf)ghcos θ. (4.6) The first term can be expressed as 1 2gh2cos θdiv(ρsϕv +ρf(1 −ϕ)u)= 1 2ghcosθdiv (h(ρsϕv +ρf(1 −ϕ)u)) −1 2ghcosθ(ρsϕv +ρf(1 −ϕ)u)·∇h. (4.7) However, according to (4.2)and(4.3)wehave ∂t1 2gh2cos θ(ρsϕ+ρf(1 −ϕ))=−1 2ghcos θdiv (h(ρsϕv +ρf(1 −ϕ)u)) −1 2ghcos θ(ρsϕ+ρf(1 −ϕ)) div (h(ϕv +(1−ϕ)u)) , Thus using this in (4.7)weobtain 1 2gh2cos θdiv(ρsϕv +ρf(1 −ϕ)u) =−∂t1 2gh2cos θ(ρsϕ+ρf(1 −ϕ))−1 2ghcosθ(ρsϕ+ρf(1 −ϕ)) div(hv) −1 2ghcos θ(ρsϕv +ρf(1 −ϕ)u)·∇h. 118 F. BOUCHUT ET AL. Adding the second term in (4.6) yields 1 2gh2cos θdiv(ρsϕv +ρf(1 −ϕ)u)+ϕ(ρs−ρf)gcos θdiv h2 2v =−∂t1 2gh2cos θ(ρsϕ+ρf(1 −ϕ))−1 2ghcosθ(ρsϕ+ρf(1 −ϕ)) div(hv) −1 2ghcosθ(ρsϕv +ρf(1 −ϕ)u)·∇h+ϕ(ρs−ρf)gcos θdiv h2 2v =−∂t1 2gh2cos θ(ρsϕ+ρf(1 −ϕ))−1 2ρfghcosθ(div(hv)+(ϕv +(1−ϕ)u)·∇h). (4.8) Then we can compute −1 2ρfghcosθ(div(hv)+(ϕv +(1−ϕ)u)·∇h) =−1 2ρfghcosθdiv(hv)−div 1 2ρfgh2cos θ(ϕv +(1−ϕ)u) +1 2ρfghcosθdiv (h(ϕv +(1−ϕ)u)) =−div 1 2ρfgh2cos θ(ϕv +(1−ϕ)u)+1 2ρfghcosθdiv (h(1 −ϕ)(u−v)) . Plugging this into (4.8)and(4.6), we obtain T1=−∂t1 2gh2cos θ(ρsϕ+ρf(1 −ϕ))−div 1 2ρfgh2cos θ(ϕv +(1−ϕ)u) −˜ βh|u−v|2−|v|tan δϕ(ρs−ρf)ghcos θ. Using this result in (4.5) finally yields the energy identity ∂tρsϕh|v|2 2+ρf(1 −ϕ)h|u|2 2+ghcos θ(b+˜ b+h 2)(ρsϕ+ρf(1 −ϕ)) +divρsϕh|v|2 2v+ρf(1 −ϕ)h|u|2 2u+(1−ϕ)h(u−v)pfbed −ρf(1 −ϕ)gh2cos θ(u−v) +ghcosθ(b+˜ b+h 2)(ρsϕv +ρf(1 −ϕ)u)+1 2ϕgh2cos θ(ρs−ρf)v +1 2ρfgh2cos θ(ϕv +(1−ϕ)u)=Re, (4.9) with Re=−˜ βh|u−v|2−|v|tan δϕ(ρs−ρf)ghcos θ. (4.10) Therefore the model (4.1) has a locally dissipative energy balance, since the residual Reis non-positive. A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 119 Remark 4.1. Identity (4.9) can be obtained (up to O(3)) by integration of (2.11) with respect to z.This shows that the left-hand side contains the physically relevant energy and energy flux. Remark 4.2. Let us recall that the Pitman−Le model [38] does not use any closure equation (4.1e). Instead, the Pitman−Le model, in the form proposed by Pelanti et al. [36], can be seen as (4.1)wherewesetpfbed =ρfghcosθ (or equivalently ps|b+h= 0 according to (3.10)). Consequently, the energy equation satisfied by the Pitman−Le model is (4.9) with a right-hand side Rethat is not always non-positive. The residual term for the Pitman−Le model is Re=−1 2ϕ(ρs−ρf)ghcos θdiv (h(1 −ϕ)(u−v)) −˜ βh|u−v|2−|v|tan δϕ(ρs−ρf)ghcos θ and has no fixed sign (we will study this term in Test 1 presented in Sect. 5.2.1). The intrinsic reason why the Pitman−Le model has a physically irrelevant energy equation is that it is derived from a 3D model that does not have an energy dissipation principle (see Eq. (2.11)). 4.2. Other properties Model (4.1) is a balance law type system with non-local terms related to pfbed .Notefirstthatitispossible to eliminate pfbed from the system. Indeed, pfbed appears only in (4.1c)and(4.1d). We can thus retain the sum of (4.1c)and(4.1d) and if we express ∇pfbed from (4.1c) for example and write that the curl of the result vanishes, we obtain the missing relation. Proposition 4.3. System (4.1)has the following properties. (i) The two mass equations are conservative. The momentum equations take the quasi-conservative form ρs(∂t(hϕv)+div(hϕv ⊗v)) = h(1 −ϕ)∇pfbed −(1 −ϕ)ρfgcos θ∇h −ϕρsgcos θ∇(b+h)−1 2(ρs−ρf)ghcos θ∇ϕ −ϕρsgsin θ(1,0)t+˜ β(u−v) −sign(v) tan δϕ(ρs−ρf)g cos θ, (4.11a) ρf(∂t(h(1 −ϕ)u)+div(h(1 −ϕ)u⊗u)) = h−(1 −ϕ)∇pfbed −(1 −ϕ)ρfgcos θ∇b −(1 −ϕ)ρfgsin θ(1,0)t−˜ β(u−v).(4.11b) The total momentum takes the conservative form ρs(∂t(hϕv)+div(hϕv ⊗v)) + ρf(∂t(h(1 −ϕ)u)+div(h(1 −ϕ)u⊗u)) =−∇ (ρsϕ+ρf(1 −ϕ))gh2 2cos θ−gcos θ(ρsϕ+ρf(1 −ϕ)) h∇(b+˜ b) −sign(v) tan δϕ(ρs−ρf)gh cos θ. (4.12) (ii) The thickness hremains non-negative, and 0≤ϕ≤1. 120 F. BOUCHUT ET AL. (iii) Special solutions are the steady states at rest, characterized by u=v=0,b+˜ b+h=Cst, ϕ =Cst, (4.13) where we recall that ˜ b=xtan θ. Indeed this is a solution to our system with pfbed =ρfghcosθ. (iv) The classical single fluid shallow water system is obtained when u=v,ρf=ρs=ρand pfbed =ρgh cosθ, and either ϕ=1in equation (4.11a)or ϕ=0in equation (4.11b). 5. Numerical approximation In this section we describe a numerical method to approximate the proposed two-phase model (4.1)inone dimension. Then we perform different tests, including a comparison with the solution provided by the Pitman−Le model. We focus on the one-dimensional situation. As pointed out previously, the model can be rewritten in terms ofthesolidpressureatthefreesurfaceps|b+h, the fluid pressure at the free surface pf|b+h=−ps|b+hor the solid pressure at the bed psbed, instead of the fluid pressure at the bed pfbed,via relations (3.9)and(3.10). In this section we consider the formulation in terms of the solid pressure at the free surface ps|b+h,whichcanbe written ∂t(hϕ)+∂x(hϕv)=0,(5.1a) ∂t(h(1 −ϕ)) + ∂x(h(1 −ϕ)u)=0,(5.1b) ∂t(hϕv)+∂x(hϕv2)=−h(1 −ϕ)∂xψ−ϕgh cos θ∂ x(b+h) −1 2(1 −r)gh2cos θ∂ xϕ −ϕgh sin θ+ˆ βh(u−v), −sign(v) tan δgcosθ(1 −r)hϕ, (5.1c) ∂t(h(1 −ϕ)u)+∂x(h(1 −ϕ)u2)=h r(1 −ϕ)∂xψ−(1 −ϕ)ghcos θ∂ x(b+h) −(1 −ϕ)ghsin θ−1 r˜ βh(u−v),(5.1d) ∂x(h(1 −ϕ)(u−v)) = 0,(5.1e) where r=ρf/ρs,ˆ β=˜ β ρs ,ψ=ps|b+h ρs =ρfghcos θ−pfbed ρs ·(5.2) If we consider that the drag coefficient ˜ βis defined by (3.12), then ˆ β=(1 −r)ϕg vT(1 −ϕ)m−1·(5.3) 5.1. Numerical method We apply a splitting algorithm, similar to the Teman-Chorin method for incompressible Euler equations, in order to impose the constraint (5.1e). We observe that at the first step, when we neglect the extra unknown ψ(which can be seen as a Lagrange multiplier associated with the constraint) in (5.1a)−(5.1d), we obtain the Pitman−Le model in the form proposed in [36], which is hyperbolic whenever u−vis not too large. We consider the space domain [0,L] divided in cells Ii=(xi−1/2,x i+1/2). For simplicity, we assume that these cells have a constant size Δx. We define xi+1 2=iΔx and xi=(i−1/2)Δx, the center of the cell Ii.Let Δt be the time step and define tn+1 =tn+Δt.Wis the vector of the following unknowns of the problem, W=[hϕ, h(1 −ϕ),hϕv,h(1 −ϕ)u].(5.4) A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 121 Therefore Wn idenotes the approximation provided by the numerical scheme of the cell averages of the solution, Wn i∼ =1 Δx xi+1/2 xi−1/2 W(tn,x)dx, (5.5) and by ψn i+1/2, an approximation of ψ(tn,x i+1/2). Assuming that the values of Wn iare known, the system can be discretized in two steps. •First step. We compute the state W∗=[h∗ϕ∗,h ∗(1−ϕ∗),h ∗ϕ∗v∗,h ∗(1−ϕ∗)u∗] by a semi-implicit discretization for the drag W∗ i=Wn i−Δt ΔxLWn i−1,Wn i,Wn i+1,Δ(b+˜ b)i−1/2,Δ(b+˜ b)i+1/2 +0,0,Δtˆ β∗ ih∗ i(u∗ i−v∗ i),−Δtˆ β∗ i h∗ i r(u∗ i−v∗ i),(5.6) where L(Wn i−1,Wn i,Wn i+1,Δ(b+˜ b)i−1/2,Δ(b+˜ b)i+1/2) defines the space discretization operator applied to model (5.1a)−(5.1d)withψ=0andˆ β= 0. In this work, we have considered the generalized Roe method proposed in [34]. Another possibility would be to use the relaxation solver proposed by Pelanti et al. [37]. •Second step. In order to enforce the constraint, we set hn+1 i=h∗ iand ϕn+1 i=ϕ∗ iand vn+1 i,un+1 iand ψn+1 i+1/2 are solutions to the following coupled system, ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (hϕv)n+1 i=(hϕv)∗ i−Δt Δx(1 −ϕ∗ i)h∗ i(ψn+1 i+1/2−ψn+1 i−1/2), (h(1 −ϕ)u)n+1 i=(h(1 −ϕ)u)∗ i+Δt Δx(1 −ϕ∗ i)h∗ i r(ψn+1 i+1/2−ψn+1 i−1/2), (h(1 −ϕ)(u−v))n+1 i+1 −(h(1 −ϕ)(u−v))n+1 i=0. (5.7) By extracting vn+1 iand un+1 ifrom the two first equations of (5.7) and by substitution in the third equation, we obtain the following system with unknowns {ψn+1 i+1/2}i, −a∗ i+1ψn+1 i+3/2+(a∗ i+a∗ i+1)ψn+1 i+1/2−a∗ iψn+1 i−1/2=h(1 −ϕ)(u−v)∗ i+1 −h(1 −ϕ)(u−v)∗ i,(5.8) with a∗ i=Δt Δxh∗ i(1 −ϕ∗ i)1 r+1−ϕ∗ i ϕ∗ i·(5.9) Thus, in this second step, we must solve system (5.8) (with Dirichlet boundary conditions, ψ=0)toobtain {ψn+1 i+1/2}iand use these values to update vn+1 iand un+1 iby the two first equations of (5.7). The obtained scheme is obviously well-balanced with respect to the steady states at rest (4.13) if the hyperbolic solver Lis well-balanced, and preserves the natural bounds h≥0and0≤ϕ≤1. 5.2. Numerical tests We will now present some numerical tests in order to compare the solution of the proposed model (5.1a)−(5.1e) with the solution of the modified Pitman and Le problem proposed in [36] (with the same drag coefficient ˜ β). We simulate the collapse of a column made of a mixture of grains and water first over a horizontal plane and then over an inclined plane, a situation widely investigated for dry granular flows (see for example [8,25,27]). First, we simulate the flow of the mixture over a horizontal bed for the two different drag forces given by (2.18) and (2.19). In the second test, we simulate the flow of the mixture over an inclined bed of constant slope for a fixed choice of these parameters. 122 F. BOUCHUT ET AL. As general considerations, we fix the CFL number as 0.8, acceleration due to gravity g=9.81 m s−2and the material densities ρf= 1000 kg m−3and ρs= 2500 kg m−3, respectively. Therefore the ratio of densities is r=0.4. 5.2.1. Test 1: Flat bottom In this experiment, the space domain is Ω=[0,10]m and we consider 200 points. At time t= 5 s, the solid phase is stopped and some small velocities appear in the fluid phase. The initial conditions are defined as follows h(t=0s)=0.5m4m≤x≤6m 0.1motherwise ;u(t=0s)=v(t=0s)=0ms −1;ϕ(t=0s)=ϕ0. For the initial solid volume fraction, we consider two values of ϕ0(0.3 and 0.6) to see the effect the initial state of mixture fluidization. We consider the intergranular Coulomb angle to be δ=18 o. The objective of this test is to check the influence of the drag force and initial solid volume fraction on the flow and deposit. Remember that the drag friction laws used here (see Sect. 2.3)aregivenby: f=˜ β(u−v), where the drag coefficient ˜ βcan be set according to: •Richardson and Zaki [41]: ˜ β=(ρs−ρf)ϕg vT(1 −ϕ)m−1,with m∈[0.4,2.65]. •Pailha and Pouliquen [35]: ˜ β=150μϕ2 d2(1 −ϕ). We also set vT=0.143 m s−1,μ=10 −3Pa s,d=10 −3m and we vary the coefficient m, using the values m=0.4,1and2.65. In the following, we will refer to “RZ” and “PP” for the Richardson and Zaki and for the Pailha and Poliquen drag forces respectively. Influenceofthedragforce. Figures 4and 5show the thickness of the mass (i.e. at time t= 5 s) and the associated volume fraction simulated with different drag forces both with the Pitman−Le (PL) model and with the new model proposed here. At that time, the solid phase has completely stopped. In Figure 4we also compare these two-phase flow models with the results obtained with the Savage-Hutter model, where the fluid phase is not considered (i.e. dry granular flows). Even for ϕ0=0.6, Figure 4show the strong influence of the fluid phase on the avalanche thickness profile. The first observation is that the PL model and our model have the same qualitative behaviour. However, for ϕ0=0.3, the PL model is more sensitive to the different drag forces introduced in the model (Figs. 4aand5a). In particular, the final volume fraction is higher for the PP drag force than for the RZ drag force at a centered mass interval, x∈(2m, 8m). They coincides just in the center of mass (x=5m) for the case m=0.4, reaching up to ϕ=0.65. For the RZ drag forces, the variation of the volume fraction in space is smoother than for the PP drag force. For our model, the drag force only slightly affects the results for ϕ0=0.3. For ϕ0=0.6, the sensitivity of the two models to the different drag forces is qualitatively similar even though the volume fraction calculated with the PL model is still more sensitive than that calculated with our model (Figs. 5b,d). For ϕ0=0.6, the volume fraction at the center of the column reaches very high values with the PL model A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 123 Figure 4. Test 1: thickness profile of the mass h(m) as a function of the distance x(m) at time t= 5 s, when the granular phase has already stopped, for the collapse of a rectangular granular mixture over a horizontal layer made of the same mixture, simulated with different friction laws (“RZ” refers to the Richardson and Zaki and “PP” to the Pailha and Pouliquen drag forces). The initial volume fractions are: (left) ϕ0=0.3; (right) ϕ0=0.6. The thickness profile obtained using the dry granular flow model of Savage and Hutter (obtained by setting all the terms related to the fluid phase equal to zero) is also represented for comparison. (about 0.9), while ϕ<0.75 with our model. The overall lower sensitivity of our model to the different drag forces suggests that the difference between the velocities of the two phases is lower with our model than with the PL model, as shown in Figure 8. Comparison of the two models. Let us now compare the two models for a given friction law, i.e., the Richardson and Zaki law with m=1 (corresponding to the data used in [10]). We also consider two different values for the Coulomb friction angle: δ=18 ◦and δ=28 ◦. For both the PL model and our model, the higher the initial fluidization, the larger the spreading of the material and the smaller the aspect ratio of the deposit (Fig. 6). Furthermore, the volume fraction is highest at the center of mass and decreases toward the front, leading to ϕ<ϕ 0at the front. 124 F. BOUCHUT ET AL. Figure 5. Test 1: solid volume fraction of the mass ϕas a function of the distance x(m) at time t= 5 s, when the granular phase has already stopped, for the collapse of a rectangular granular mixture over a horizontal layer made of the same mixture, simulated with different friction laws (“RZ” refers to the Richardson and Zaki and “PP” to the Pailha and Pouliquen drag forces). The initial volume fractions are: (left) ϕ0=0.3; (right) ϕ0=0.6. At time t= 5 s, when the solid phase is already at rest, the thickness of the mass is very similar in both models, even though the maximum thickness is slightly smaller with our model (Fig. 6). However, the distribution of the phases (i.e. volume fraction) is different. As observed previously, the volume fraction is more uniformly distributed in the simulations with our model. The peak of high volume fraction at the center of mass is higher for the PL model and the decrease in volume fraction toward the front of the mass is larger than in our model (Fig. 6c). The velocity of the fluid phase at an intermediate time t=1.5 s, when the granular phase is still flowing, is slightly higher with our model for δ=18 ◦for both ϕ0=0.3andϕ0=0.6 (see Figs. 7aand7b). For δ=28 ◦,the fluid velocities are almost the same in the two models towards the front and the fluid velocity with our model is lower around the center than that calculated with the PL model. The difference between the two models is greater for the solid velocity for both ϕ0=0.3andϕ0=0.6, (see Figs. 7cand7d). For the proposed model, the solid phase moves faster than for the PL model. We also observe that for δ=18 ◦, the velocities of both phases are greater than for δ=28 ◦. Moreover, for larger values of δ, the difference between the velocities of the A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 125 Figure 6. Test 1: comparison between the solutions obtained with the Pitman−Le and proposed models for the thickness of the mass h(m) and the volume fraction ϕas functions of the distance x(m)attimet= 5 s when the solid phase has already stopped, for the collapse of a rectangular granular mixture over a horizontal layer made of the same mixture. The initial volume fractions are: (left) ϕ0=0.3; (right) ϕ0=0.6. Here the friction law is the Richardson and Zaki drag force (RZ) with m=1. twophasesisgreater(seeFig.8). This difference is much larger in the PL model than in our model, leading to higher drag forces in the PL model. This explain why our model is less sensitive to the definition of the drag force (see Fig. 4). Evolution in time. Let us now look at the changes of the different quantities with time (Figs. 9−13 in which times t=1,2,3and 5 s are represented with different colours). For these simulations, we use the Richardson and Zaki friction law with m= 1 and the Coulomb friction angle δ=18 ◦. Note that even though the mass profiles change with time in a similar way for the two models (Fig. 9), there is strong difference between the two models for the changes of the volume fraction with time (Fig. 10), 132 F. BOUCHUT ET AL. Figure 13. Test 1: the variable ψ(x, t)inm 2s−2in the proposed model at different times for the collapse of a rectangular granular mixture over a horizontal layer made of the same mixture. The initial volume fractions are: (left) ϕ0=0.3; (right) ϕ0=0.6. Here the friction law is the Richardson and Zaki drag force (RZ) with m=1. Figure 14. Test 1: values of the terms involved in the residual energy (see Eq. (5.10)) as a function of the distance x(m), for the collapse of a rectangular granular mixture over a horizontal layer made of the same mixture at the different times t=0.5kwithk=1s,...10 s. The initial volume fraction is ϕ0=0.3 and the friction law is the Richardson and Zaki drag force (RZ) with m=1. A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 133 Figure 15. Test 1: values of the total residual term as a function of the distance x(m) for the collapse of a rectangular granular mixture over a horizontal layer made of the same mixture at the different times t=0.5kwithk=1s,...10 s. The initial volume fraction is ϕ0=0.3and the friction law is the Richardson and Zaki drag force (RZ) with m=1. Figure 16. Test 2: comparison of the thickness profiles of the mass h(x, t) in meters as a function of the distance x(m) for the Pitman−Le and proposed models, at time t=1,3,5,10 s (at t= 10 s the solid phase has already stopped), for the collapse of a rectangular granular mixture over an inclined layer (θ=10 ◦) made of the same mixture. The initial volume fraction is ϕ0=0.6. Here the friction law is the Richardson and Zaki drag force (RZ) with m=1. 134 F. BOUCHUT ET AL. Figure 17. Test 2: comparison of the solid volume fraction of the mixture ϕ(x, t) as a function of the distance x(m) for the Pitman−Leandproposedmodels,attimet=1,3,5,10 s (at t= 10 s the solid phase has already stopped), for the collapse of a rectangular granular mixture over an inclined layer (θ=10 ◦) made of the same mixture. The initial volume fraction is ϕ0=0.6. Here the friction law is the Richardson and Zaki drag force (RZ) with m=1. added to Jackson’s model to obtain a well-posed system. This may not be apparent with a depth-averaged model with hydrostatic pressure, such as the one proposed by Pitman and Le [38]. Indeed, if we assume hydrostatic pressure for both phases, they are related by their boundary condition. In this case, imposing zero atmospheric pressure for both phases can be seen as the corresponding closure equation. Nevertheless, the model that is deduced does not have a dissipative energy balance. The main difference between the model that we propose in this paper and the Pitman−Le model comes from the boundary condition on the free surface. In the proposed model, we only impose that the sum of the pressures of the two phases is zero and not each of them. This introduces new unknown in the simplified model. As a closure equation for the 3D system, we consider incompressibility of the solid phase. This closure relation is consistent with the hydrostatic pressure assumption. The numerical tests presented here show that, overall, the changes of the profiles of the flowing mass with time are similar for the Pitman−Le model and the model proposed here. The qualitative behaviour of the solid volume fraction and A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 135 Figure 18. Test 2: comparison of the velocity of the fluid phase u(x, t)inms −1as a function of the distance x(m) for the Pitman−Le and proposed models, at time t=1,3,5,10 s (at t=10s the solid phase has already stopped), for the collapse of a rectangular granular mixture over an inclined layer (θ=10 ◦) made of the same mixture. The initial volume fraction is ϕ0=0.6. Here the friction law is the Richardson and Zaki drag force (RZ) with m=1. the solid and fluid velocities is the same for both models. However, with the model presented here, the solid volume fraction varies less, the solid phase velocity is generally higher and the difference between the velocities of the two phases is smaller, leading to smaller drag forces between the two phases. This induces significant differences in the profile of the spreading mass and of the deposit with runout distances more than 10% larger and velocities that could be more than 20% heigher in the simple test of granular collapse over inclined plane performed here. While it is quite difficult to measure experimentally or in the field the fluid and solid velocities, the use of seismic waves generated by debris flows or avalanches may be a new way to discriminate between these two models ([9,31]). An advantage of our model is that the closure equation (i.e. incompressibility of the solid phase) is explicitly imposed, making it possible to derive physical interpretation of our results while in the PL model, the behavior is dictated by the imposed zero-pressure at the surface of each phase, without any description of the mechanical 136 F. BOUCHUT ET AL. Figure 19. Test 2: comparison of the velocity of the solid phase v(x, t)inms −1as a function of the distance x(m) for the Pitman−Leandproposedmodels,attimet=1,3,5,7s(att=10s the solid phase has already stopped), for the collapse of a rectangular granular mixture over an inclined layer (θ=10 ◦) made of the same mixture. The initial volume fraction is ϕ0=0.6. Here the friction law is the Richardson and Zaki drag force (RZ) with m=1. properties of the solid phase. It is interesting to see the significant differences between the two models even though the surface pressure of each phase is very small in our model. This analysis is largely driven by the kinematic boundary conditions that imposes the two phases to fill the same domain. However in debris flows, the fluid phase surface can be higher or lower than the solid phase surface, due to the relative motion between these two phases. This is expected to be significant in particular when compression/dilation of the solid phase occur. Because the models seem to be very sensitive to what happens at the surface, even small variations of the fluid and solid surfaces could have a strong impact on the results. Further analysis of the equations should be performed with relaxation of these boundary conditions and including a more realistic closure relation related to the compression/dilation of the granular phase. A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 137 Figure 20. Test 2: comparison of the velocity difference of the fluid u(x, t) and solid phase v(x, t)inms −1as functions of the distance x(m) for the Pitman−Le and proposed models, at time t= 1 and 5 s, for the collapse of a rectangular granular mixture over an inclined layer (θ=10 ◦) made of the same mixture. The initial volume fraction is ϕ0=0.6. Here the friction law is the Richardson and Zaki drag force (RZ) with m=1. Figure 21. Test 2: the mass thickness h(x, t) in meters as a function of the distance x(m) (a-b) and of the solid volume fraction ϕ(x, t) as a function of the distance x(m) (c-d) at different times for the collapse of a rectangular granular mixture over an inclined layer (θ=10 ◦)made of the same mixture. The initial volume fraction is ϕ0=0.6. Here the friction law is the Richardson and Zaki drag force (RZ) with m=1. 138 F. BOUCHUT ET AL. Figure 22. Test 2: the velocity of the fluid phase u(x, t)inms −1as a function of the distance x(m) (a-b) and of the solid phase v(x, t)inms −1as a function of the distance x(m) (c-d) at different times for the collapse of a rectangular granular mixture over an inclined layer (θ=10 ◦) made of the same mixture. The initial volume fraction is ϕ0=0.6. Here the friction law is the Richardson and Zaki drag force (RZ) with m=1. Figure 23. Test 2: left: the surface pressure of the solid phase ψ(x, t)inm 2s−2as a function of the distance x(m) at different times for the collapse of a rectangular granular mixture over an inclined layer (θ=10 ◦) made of the same mixture. Right: comparison between ψand g(b+h1+h2)att=10s. A TWO-PHASE SHALLOW DEBRIS FLOW MODEL WITH ENERGY BALANCE 139 Acknowledgements. We thank Farhang Radjai, Renaud Toussaint and Olivier Pouliquen for interesting discussions during this work. We thank Nicolas Mangold and Oldrich Hungr for leading the field trips in which we took the pictures in Figure 1. The work of E.D. Fern´andez-Nieto and G. Narbona-Reina was partially supported by the Spanish Government and FEDER through the Research projects MTM2009-07719 and MTM2012-38383-C02-02 and by the Andalusian Government through the project P11RNM7069. Anne Mangeney and Francois Bouchut were funded by the ANR PLANETEROS and LANDQUAKES programmes. References [1] T.B. Anderson and R. Jackson, A fluid mechanical description of fluidized beds. Ind. Eng. Chem. Fundam. 6(1967) 527–539. [2] K. Anderson, S. Sundaresan and R. Jackson, Instabilities and the formation of bubbles in fluidized beds. J. Fluid Mech. 303 (1995) 327–366. [3] T. Aste, Circle, sphere, and drop packings. Phys. Rev. E 53 (1996) 2571. [4] F. Boyer, O. Pouliquen and E. Guazzelli, Dense suspensions in rotating-rod flows: normal stresses and particle migration. J. Fluid Mech. 686 (2011) 5–25. [5] F. Bouchut and M. Westdickenberg, Gravity driven shallow water models for arbitrary topography. Commun. Math. Sci. 2 (2004) 359–389. [6] F. Bouchut, A. Mangeney-Castelnau, B. Perthame and J.-P. Vilotte, A new model of Saint Venant and Savage-Hutter type for gravity driven shallow water flows. C.R. Acad. Sci. Paris, S´er. I 336 (2003) 531–536. [7] F. Bouchut, E. Fernandez-Nieto, A. Mangeney and P.-Y. Lagr´ee, On new erosion models of Savage-Hutter type for avalanches. Acta Mech. 199 (2008) 181–208. [8] M. Farin, A. Mangeney and O. Roche, Dynamics, deposit and erosion processes in granular collapse over sloping beds. J. Geophys. Res. 119 (2013) 504–532. [9] P. Favreau, A. Mangeney, A. Lucas, G. Crosta and F. Bouchut, Numerical modeling of landquakes. Geophys. Res. Lett. 37 (2010) L15305. [10] D.L. George and R.M. Iverson, A two-phase debris-flow model that includes coupled evolution of volume fractions, granular dilatancy, and pore-fluid pressure. ItalianJ.Engrg.Geol.Environ. (2011) DOI:10.4408/IJEGE.2011-03.B-047. [11] M. Goodman and S. Cowin, Two problems in the gravity flow of granular materials. J. Fluid Mech. 45 (1971) 321–339. [12] J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics. Kluwer Academic Publishers (1983). [13] H.J. Herrmann, R. Mahmoodi Baram and M. Wackenhut, Searching for the perfect packing. Phys. A 330 (2003) 77–82. [14] O. Hungr and S.G. Evans, Entrainment of debris in rock avalanches: An analysis of a long run-out mechanism. Bull. Geol. Soc. Am. 116 (2004) 1240–1252. [15] R.M. Iverson, The physics of debris flows, Rev. Geophys. 35 (1997) 245–296. [16] R.M. Iverson, R.P. Denlinger, Flow of variably fluidized granular masses across three-dimensional terrain 1: Coulomb mixture theory. J. Geophys. Res. 106 (2001) 537–552. [17] R.M. Iverson, M. Logan, R.G. LaHusen and M. Berti, The perfect debris flow? aggregated results from 28 large-scale experiments. J. Geophys. Res. 115 (2010) F03005. DOI:10.1029/2009JF001514. [18] R.M. Iverson, M.E. Reid, M. Logan, R.G. LaHusen, J.W. Godt and J.P. Griswold, Positive feedback and momentum growth during debris-flow entrainment of wet bed sediment. Nature Geoscience 4(2011) 116–121. [19] R. Jackson, The Dynamics of Fluidized Particles. Cambridges Monographs on Mechanics (2000). [20] G. Jhonson, M. Massoudi and K.R. Rajagopal, A review of interaction mechanisms in fluidsolid flows. Technical Report. DOE/PETC/TR-90/9, U.S. Dep. of Energy, Pittsburgh Energy Tech. Ctr., Pittsburgh, USA (1990). [21] F. Legros, The mobility of long-runout landslides. Eng. Geol. 63 (2002) 301–331. [22] A. Lucas and A. Mangeney, Mobility and topographic effects for large Valles Marineris landslides on Mars. Geophys. Res. Lett. 34 (2007) L10201. [23] A. Lucas, A. Mangeney and J.P. Ampuero, Frictional weakening in landslides on Earth and on other planetary bodies. Nature Commun. 5(2014) 3417. [24] A. Mangeney-Castelnau, J.P. Vilotte, M.O. Bristeau, B. Perthame, F. Bouchut, C. Simeoni and S. Yernini, Numerical modeling of avalanches based on Saint-Venant equations using a kinetic scheme. J. Geophys. Res. 108 (2003) 2527. [25] A. Mangeney-Castelnau, F. Bouchut, J.P. Vilotte, E. Lajeunesse, A. Aubertin and M. Pirulli, On the use of Saint-Venant equations for simulating the spreading of a granular mass. J. Geophys. Res. 110 (2005) B09103. [26] A. Mangeney, L. S. Tsimring, D. Volfson, I.S. Aranson and F. Bouchut, Avalanche mobility induced by the presence of an erodible bed and associated entrainment. Geophys. Res. Lett. 34 (2007) L22401. [27] A. Mangeney, O. Roche, O. Hungr, N. Mangold, G. Faccanoni and A. Lucas, Erosion and mobility in granular collapse over sloping beds. J. Geophys. Res.-Earth Surf. 115 (2010) F03040. [28] A. Mangeney, Landslide boost from entrainment. Nature Geosci. 4(2011) 77–78. [29] N. Mangold, A. Mangeney, V. Migeon, V. Ansan, A. Lucas, D. Baratoux and F. Bouchut, Sinuous gullies on Mars: Frequency, distribution, and implications for flow properties. J. Geophys. Res. Planets 115 (2010) E11001. [30] C. Meruane, A. Tamburrino, O. Roche, On the role of the ambient fluid on gravitational granular flow dynamics. J. Fluid. Mech. 648 (2010) 381–404. 140 F. BOUCHUT ET AL. [31] L. Moretti, A. Mangeney, Y. Capdeville, E. Stutzmann, C. Christian Huggel, D. Schneider and F. Francois Bouchut, Numerical modeling of the Mount Steller landslide flow history and of the generated long period seismic waves. Geophys. Res. Lett. 39 (2012) L16402. [32] M.J. Niebling, E.G. Flekkoy, K.J. Mˆaloy and R. Toussaint, Mixing of a granular layer falling through a fluid.Phys. Rev. E 82 (2010) 011301. [33] M. Ouriemi, P. Aussillous and E. Guazzelli, Sediment dynamics. Part I: Bed-load transport by shearing flows. J. Fluid Mech. 636 (2009) 295–319. [34] C. Par´es and M.J. Castro, On the well-balance property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow-water systems. ESAIM: M2AN 38 (2004) 821–852. [35] M. Pailha and O. Pouliquen, A two-phase flow description of the initiation of underwater granular avalanches. J. Fluid Mech. 633 (2009) 115–135. [36] M. Pelanti, F. Bouchut and A. Mangeney, A Roe-type scheme for two-phase shallow granular flows over variable topography. ESAIM: M2AN 42 (2008) 851–885. [37] M. Pelanti, F. Bouchut and A. Mangeney, A Riemann solver for single-phase and two-phase shallow flow models based on relaxation. Relations with Roe and VFRoe solvers. J. Comput. Phys. 230 (2011) 515–550. [38] E.B. Pitman and L. Le, A two-fluid model for avalanche and debris flows. Philos.Trans.R.Soc.A363 (2005) 1573–1601. [39] S.P. Pudasaini, Y. Wang and K. Hutter, Modelling debris flows down general channels. Natural Hazards Earth System Sci. 5 (2005) 799–819. [40] L. Rondon, O. Pouliquen and P. Aussillous, Granular collapse in a fluid: role of the initial volume fraction. Phys. Fluids 23 (2011) 073301. [41] J.F. Richardson and W.N. Zaki, Sedimentation and Fluuidization: part I. Trans. Inst. Chem. Eng. 32 (1954) 35–53. [42] O. Roche, M. Attali, A. Mangeney and A. Lucas, On the run out distance of geophysical gravitational flows: insight from fluidized granular collapse experiments. Earth Planet. Sci. Lett. 311 (2011) 375–385. [43] O. Roche, Y. Ni˜no, A. Mangeney, B. Brand, N. Pollock and G.A. Valentine, Dynamic pore pressure variations induce substrate erosion by pyroclastic flows, Geology 41 (2013) 1107–1110. [44] S. Roux and F. Radjai, Texture-dependent rigid plastic behavior. In Proc.ofPhysicsofDryGranularMedia1997. Edited by H.J. Herrmann. Kluwer. Carg`ese, France (1998) 305–311. [45] S.B. Savage and K. Hutter, The motion of a finite mass of granular material down a rough incline. J. Fluid Mach. 199 (1989) 177–215. [46] U.E. Shamy and M. Zhegal, Coupled continuum discrete model for saturated granular materials. J. Engrg. Mech. 131 (2005) 413–426. [47] V. Topin, F. Dubois, Y. Monerie, F. Perales, A. Wachs: Micro-rheology of dense particulate flows: application to immersed avalanches. J. Non-Newtonian Fluid 166 (2011) 63–72. [48] C. Voivret, F. Radjai and J.Y. Delenne, M.S. El Youssoufi, Space-filling properties of polydisperse granular media. Phys. Rev. E76 (2007) 021301.