A two-phase two-layer model for fluidized granular flows with dilatancy effects

A two-phase two-layer model for fluidized granular flows with dilatancy effects Fran¸cois Bouchut∗ , Enrique D. Fern´andez-Nieto† , Anne Mangeney‡§ , Gladys Narbona-Reina† Abstract We propose a two-phase two-thin-layer model for fluidized debris flows that takes into account dilatancy effects, based on the closure relation proposed by Roux and Radjai (1998). This relation implies that the occurrence of dilation or contraction of the granular material depends on whether the solid volume fraction is respectively higher or lower than a critical value. When dilation occurs, the fluid is sucked into the granular material, the pore pressure decreases and the friction force on the granular phase increases. On the contrary, in the case of contraction, the fluid is expelled from the mixture, the pore pressure increases and the friction force diminishes. To account for this transfer of fluid into and out of the mixture, a two-layer model is proposed with a fluid layer on top of the two-phase mixture layer. Mass and momentum conservation are satisfied for the two phases, and mass and momentum are transferred between the two layers. A thin-layer approximation is used to derive average equations, with accurate asymptotic expansions. Special attention is paid to the drag friction terms that are responsible for the transfer of momentum between the two phases and for the appearance of an excess pore pressure with respect to the hydrostatic pressure. For an appropriate form of dilatancy law we obtain a depth-averaged model with a dissipative energy balance in accordance with the corresponding 3D initial system. Keywords: Fluidized granular flows, two-phase, dilatancy, two-layer, depth-averaged model, critical volume fraction, excess pore pressure 1 Introduction Gravity driven flows such as debris flows, sub-aerial and submarine landslides play a key role in erosion processes on the Earth’s surface. They represent one of the major natural hazards threatening life and property in mountainous, volcanic, seismic and coastal areas, as shown recently by the debris flows that occurred in Uganda and Brazil in 2010, causing 400 and 350 deaths respectively and displacing several hundred thousand inhabitants. One of the ultimate goals of landslide studies is to produce tools for the prediction of velocity and runout extent of rapid landslides. Developing a theoretical description and physical ∗Universit´e Paris-Est, Laboratoire d’Analyse et de Math´ematiques Appliqu´ees (UMR 8050), CNRS, UPEM, UPEC, F-77454, Marne-la-Vall´ee, France ([email protected]) †Departamento de Matem´atica Aplicada I, Universidad de Sevilla. E.T.S. Arquitectura. Avda, Reina Mercedes, s/n. 41012 Sevilla, Spain ([email protected], [email protected]) ‡Universit´e Paris Diderot, Sorbone Paris Cit´e, Institut de Physique du Globe de Paris, Equipe de Sismologie, 1 rue Jussieu, 75005 Paris, France ([email protected]u.fr) §ANGE team, INRIA, CETMEF, Lab. J. Louis Lions, Paris, France 1 understanding of the associated processes in a natural environment remains an unsolved and extremely challenging problem in Earth science, mechanics and mathematics. Recent progress in the mathematical, physical and numerical modelling of gravity driven flows has led to the development and use of numerical models for investigating geomorphological processes and assessing risks related to such natural hazards. However, severe limitations prevent us from fully understanding the physical processes acting in natural flows and from predicting landslide dynamics and deposition. One of the important issues is that existing models do not accurately account for the co-existence and interaction of fluid (water and gas) and solid granular phases within the flowing mass, which play a key role in natural gravity related instabilities. Water is almost always present in natural landslides and the frequently resulting debris flows (mixture of water and grains) are often highly destructive. The interaction between the fluid and granular phases within a saturated mixture essentially depends on the fluid pressure, also called pore pressure, that determines the effective friction force acting on the granular medium (e.g. Jackson 2000; Iverson 2000, 2005). Since the pioneered work of (Reynolds 1885), a large amount of studies have been dedicated to dilatancy effects in granular materials and to their interaction with pore fluid pressure, solid pressure and strain rates (e.g. Schofield and Wroth 1968; Jackson 1983; Vardoulakis 1986; Bolton 1986; Mitchell 1993; Wood 1990). A change in the fluid pressure may result from a dilation of the granular phase, that induces a sucking of the fluid within the mixture and a diminution of the fluid pressure, thereby increasing the effective friction on the granular phase. On the other hand, a contraction of the granular phase induces an expulsion of the fluid from the mixture and an increase of the fluid pressure, thereby decreasing the effective friction. This process is sometimes called “pore pressure feedback” (Iverson 2005). Contraction of a grain-fluid mixture may lead to liquefaction of the mixture. Dilation and contraction occur in response to a deformation of the granular medium, and in particular to shearing. Indeed, a densely packed granular assembly (high solid volume fraction) must dilate to be sheared, in order for the grains to have room enough to move one with respect to the other. On the other hand, a loosely packed assembly contracts in response to shearing. These processes play a dramatic role in the dynamics of fluidized granular flows, from their initial destabilization to their final deposition (Andreini et al. 2013; Iverson et al. 2010; George and Iverson 2011; Montserrat et al. 2012; Rondon et al. 2011; Iverson and George 2016). Taking into account dilatancy effects in numerical models of granular flows is a crucial issue. However, solving the complete 3D equations of granular mass motion, with sufficient resolution to describe the real topography, requires prohibitive computational costs. For this reason, it is necessary to write simplified models. A class of efficient techniques, developed and successfully employed to reproduce a large range of experimental and geological observations, makes use of a depth-averaged continuum description, based on the thin-layer approximation (i.e. the thickness of the flowing mass is assumed to be small compared to its downslope extension) [Savage and Hutter, 1989]. This leads to the assumption that the velocity normal to the topography is small compared to the downslope velocity. Taking into account two-phase grain-fluid mixtures and dilatancy in the thin-layer approximation raises significant mathematical difficulties because of the need for a consistent description of these effects within this approximation. In particular, contraction-dilation induces a relative motion of the fluid and solid phases in the direction normal to the topography, that is formally small in the thin-layer asymptotic expansion. The drag friction force between the fluid and solid phases is however strong enough to make it important to take this relative motion into account in the asymptotic model as detailed in this paper. 2 The solid-fluid mixture models described in the literature are generally based on Jackson’s model (Jackson 2000) that describes the main interactions between the two phases, such as buoyancy and drag frictional forces. Setting apart rheological laws, the main equations in Jackson’s model are mass and momentum conservation for the two phases, thus eight scalar equations. It has nine principal unknowns: the solid volume fraction, the solid and fluid pressures and the components of solid and fluid velocities. As a result, a scalar closure equation is necessary to complete the model. Several depth-averaged thin-layer models have been deduced from Jackson’s model (e.g. Pitman and Le 2005; Pelanti et al. 2008; Pailha and Pouliquen 2009; Kowalski and McElwaine 2013; Iverson and George 2014). Pitman and Le (2005) followed by Pelanti et al. (2008) replaced the closure relation by an extra boundary condition at the free surface. This leads to an overdetermined problem at the free surface (two kinematic conditions), and to an underdetermined problem inside the domain. However, given the hydrostatic pressure assumption, a depth-averaged model can be obtained since the disappearence of the normal variable gives a kind of equivalence between a boundary condition and a closure relation inside the domain. The lack of a relevant closure equation leads to a non-dissipative energy balance in the Pitman and Le model, as well as in its variants. Moreover, these models do not take into account dilatancy effects. See (Bouchut et al. 2015) for more details on the different methods used to tackle this problem and on the validity of the proposed closure relations. A crucial point in order to obtain a realistic model is that the energy balance associated with the model must be physically relevant. A main objective here is to propose a closure equation that gives such an energy balance, at least in the case when a simplified rheology is taken. Along this line, in our previous work (Bouchut et al. 2015) we proposed a depth-averaged two-phase debris flow model that gives a dissipative energy balance. In that model, the closure equation is simply the incompressibility of the solid phase – in the sense of cancellation of the dilation rate (divergence of the solid velocity) – so that dilatancy is not accounted for. Moreover, in order to avoid overdetermined boundary conditions, only the sum of the solid and fluid normal stresses is set to zero at the free surface, instead of both separately. We propose here to close Jackson’s model by including dilatancy effects, based on the model proposed by Roux and Radjai (1998) for dry granular flows derived from critical state mechanics (e.g. Schofield and Wroth 1968; Wood 1990). In this model, the dilation rate is directly related to the volume fraction and is taken to be equal to ˙γtan ψ, where ˙γis the shear rate and ψis the “dilation angle” that depends on the volume fraction. This description of dilatancy has been used in (Pailha and Pouliquen 2009) to develop a thin-layer depth-averaged two-phase model for immersed granular flows. In this configuration there is one moving surface for the mixture, and one fixed (say horizontal) surface for the fluid, thus there is no excess boundary condition. However the authors assumed uniformity in the slope aligned direction. In their final model, the dilatancy effect appears through an excess pore pressure term, in addition to the hydrostatic pressure. Other kinds of debris-flow models are based on the idea of a single-phase mixture model. One of the first such models was presented by Iverson (1997), followed by other versions proposed in (Iverson 2009; George and Iverson 2011; Iverson and George 2014; George and Iverson 2014), still based on a single-phase mixture model. As a result, the relative motion between the solid and fluid phases does not appear explicitly. The mass and momentum equations for the mixture are coupled to an advection-diffusion equation to describe the changes in pore pressure. The model in (Kowalski and McElwaine 2013) is also of this type, it uses a closure by the (Richardson and Zaki 1954) sedimentation law, which is an alternative way to formulate the relative motion of solid and fluid phases by an advection-diffusion equation on the volume fraction indeed. To close the system, Iverson & George assume that the mixture obeys a Darcy 3 law and they use a closure relation that takes into account the dilatancy effects. More precisely, they consider a modification of the Roux and Radjai dilatancy law in order to introduce the variations of the effective stress, already proposed in (Iverson 2009). In this case the dilation rate is given by ˙γtan ψ−αd dt (σ−pf), where αis the compressibility of the mixture, σthe total normal stress and pfthe fluid pressure. The definition of αis discussed in (Andreini et al. 2013). The aim of this paper is to establish a depth-averaged two-phase thin-layer model including dilatancy effects from Jackson’s model with the Roux and Radjai closure. It is a kind of extension with slope aligned variable dependency of the model of (Pailha and Pouliquen 2009), in the two thin-layer configuration. As opposed to previously cited works, and in order to be consistent with the physical processes described above, we consider an extra upper fluid layer, that allows the fluid to be expelled or sucked in from the mixture at its upper boundary. This also allows us to resolve the overdetermination at the boundary, because now there are two moving surfaces, and one kinematic condition for each of them. This is a key point in our approach. An accurate asymptotic analysis is performed to derive the depth-averaged system. We show that the effect of dilatancy on the fluid pressure appears through an extra contribution to the hydrostatic pressure, the so-called excess pore pressure. It is strongly related to the normal relative motion between the granular and fluid phases. We prove additionally that the proposed model satisfies a dissipative energy balance equation as well as the initial 3D starting system, under the assumption of a pressure dependent critical volume fraction. This is obtained via a compressible interpretation of our model. The paper is organized as follows. Section 2 describes the 3D starting mixture system together with closure equation and boundary conditions. The thin-layer model is derived in Section 3 where the scaling assumptions are specified. In Section 4 we discuss the properties of our thin-layer model and the differences with other models in the literature. Section 5 shows some preliminary numerical simulations in the uniform setting, and Section 6 presents our conclusions. Technical calculations are provided in several appendices. 2 Two-phase mixture model 2.1 Jackson’s model The starting point of our derivation is the same as in (Bouchut et al. 2015), i.e. the 3D model proposed by Jackson (2000) for flows of solid granular materials filled (saturated) with fluid. A related theory of mixtures is also developed in (Brenner 2009). The two mass conservation equations for the solid and fluid phases are, respectively, ∂t(ρsϕ) + ∇·(ρsϕv) = 0,(2.1a) ∂t(ρf(1 −ϕ)) + ∇·(ρf(1 −ϕ)u) = 0,(2.1b) and equations of momentum conservation for each phase are ρsϕ(∂tv+ (v·∇)v) = −∇·Ts+f0+ρsϕg,(2.2a) ρf(1 −ϕ)(∂tu+ (u·∇)u) = −∇·Tfm−f0+ρf(1 −ϕ)g.(2.2b) The velocities are vfor the solid phase and ufor the fluid phase, while Tsand Tfmdenote the (symmetric) stress tensors for the solid and the fluid, respectively. Moreover, the constant densities are denoted by ρsand ρf. Acceleration due to gravity is denoted by g,and f0represents 4 the average value of the resultant force exerted by the fluid on a solid particle. The solid volume fraction is ϕ. The combination of (2.1a) and (2.1b) yields the mass conservation for the mixture ∂t(ρm) + ∇·(ρmVm) = 0,(2.3) where ρm=ρsϕ+ρf(1 −ϕ), Vm=ρsϕv +ρf(1 −ϕ)u ρsϕ+ρf(1 −ϕ),(2.4) are the density and velocity of the mixture, respectively. Dividing (2.1a) by ρs, (2.1b) by ρfand adding the results gives ∇·(ϕv + (1 −ϕ)u) = 0,(2.5) that can be written also ∇·v=∇·((1 −ϕ)(v−u)). Note that this relation does not imply that ∇·Vmis equal to zero. According to Anderson and Jackson (1967) and as in (Bouchut et al. 2015), the force f0is decomposed into the sum of the buoyancy force fBand all remaining contributions f, f0=fB+f=−ϕ∇pfm+f, (2.6) where pfmis the fluid pressure in the mixture (pore pressure). The term fcombines the drag force, the lift force and the virtual mass force. Note that separation of the buoyancy force from the rest of inter-phase forces is not trivial as explained in (Jackson 2000). Here we assume that fcan be expressed simply by the drag force, thus f=˜ β(u−v),(2.7) ˜ βbeing the drag coefficient given as in (Pailha and Pouliquen 2009; Iverson and George 2014) by ˜ β= (1 −ϕ)2ηf κ,(2.8) where ηfis the dynamic viscosity of the fluid and κis the hydraulic permeability of the granular aggregate, that depends on ϕ. By substituting (2.6) into (2.2a) and (2.2b), we obtain ρsϕ(∂tv+ (v·∇)v) = −∇·Ts−ϕ∇pfm+f+ρsϕg,(2.9a) ρf(1 −ϕ)(∂tu+ (u·∇)u) = −∇·Tfm+ϕ∇pfm−f+ρf(1 −ϕ)g.(2.9b) Note that adding (2.9a) and (2.9b) and taking into account (2.1a), (2.1b) yields the conservation of total momentum ∂tρsϕv +ρf(1 −ϕ)u+∇·ρsϕv ⊗v+ρf(1 −ϕ)u⊗u+Ts+Tfm =ρsϕ+ρf(1 −ϕ)g.(2.10) We shall assume rheologies of the form Ts=psId + e Ts, Tfm=pfmId + e Tfm,(2.11) where psand pfmare the total pressures for the solid and fluid within the mixture, respectively, and e Ts,e Tfmneed to be defined, according to rheological assumptions. With (2.11), the fluid momentum equation (2.9b) involves the pressure contribution −(1 −ϕ)∇pfm. It is important to see that, since the factor 1 −ϕappears also on the left-hand side, the velocity uhowever only feels the term −∇pfm, as expected since pfmis the pore pressure. The interpretation of the 5 solid momentum equation (2.9a) is that the solid feels the buoyancy term −ϕ∇pfmand the solid pressure term −∇ps. The latter pressure ps(also called effective normal stress) represents only the effects of grains interactions, and its gradient can be evaluated in practice by measuring the force exerted on a grid immersed into the mixture. The system of eight scalar equations (2.1a), (2.1b), (2.9a), (2.9b) has nine scalar unknowns ϕ,ps,pfm, and the components of uand v. Thus, as exposed in (Bouchut et al. 2015), it is not closed, and this is due to the averaging process used for its deduction (see (Jackson 2000) for details). Therefore, a closure relation is needed, under the form of an additional scalar equation that should be imposed, based on the physical processes involved. A possible closure is to impose the incompressibility of the solid phase, ∇ · v= 0, considered in the previous work (Bouchut et al. 2015). But in real granular materials the dilatancy effects, due to geometrical congestion, may induce changes of the solid dilation rate ∇ · v, even if the mass of the granular material remains constant. This effect has to be included in the model instead of incompressibility. 2.2 Closure and energy balance The energy balance associated to Jackson’s system can be written, as in (Bouchut et al. 2015), ∂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 +pfmϕv + (1 −ϕ)u+e Tfmu+Tsv =Ts:∇v+e Tfm:∇u+f·(v−u), (2.12) where Xdenotes the space position. The friction effects give naturally a dissipative term f· (v−u)≤0, and it is also natural to assume that e Tfm:∇u≤0. The sign of Ts:∇vremains however undetermined. Since by (2.11) Ts:∇v=ps∇·v+e Ts:∇v, (2.13) and it is also natural to have e Ts:∇v≤0, it remains the term ps∇·v. As mentioned above, the closure relation that states the incompressibility of the solid phase ∇·v= 0 gives a consistent energy balance and the model of (Bouchut et al. 2015), but does not take into account dilatancy. Thus we consider the following closure equation to Jackson’s model, involving the solid dilation rate ∇·v, ∇·v= Φ,(2.14) with Φ a function to be determined, that may depend on the unknowns of the system, as discussed in the next subsection. This kind of “weakly compressible” closure is considered in low Mach number flows, see for example (Penel et al. 2015). This equation (2.14) together with (2.1a), (2.1b), (2.9a), (2.9b), (and (2.11) with suitable definitions of e Ts,e Tfm), gives a closed system. Then in the right-hand side of (2.12) with the decomposition (2.13), only the first term psΦ is not always nonpositive. This term is further analyzed in Subsection 2.4. 2.3 Dilatancy in dense granular flows In the work of Roux and Radjai (1998), a model for introducing dilatancy effects into the behaviour of dry granular media is proposed. This effect is directly related to the changes 6 experimented by the solid volume fraction. In particular, the rate of volume change is given by ˙γtan ψ, where ˙γ=|Dv|is the norm of the strain rate Dv = (∇v+∇vt)/2, and ψis the so called “dilation angle”. This means more explicitly that ∂tϕ+v·∇ϕ=−ϕ˙γtan ψ. (2.15) From the mass equation (2.1a) we have ∂tϕ+v·∇ϕ=−ϕ∇·v, thus we can reformulate (2.15) as a relation between the solid dilation rate ∇·vand the dilation angle ψ, as ∇·v= ˙γtan ψ. (2.16) The dilation angle ψis in turn related to the solid volume fraction ϕ, and a linear approximation can be written ψ=a(ϕ−ϕeq c), with a > 0, and ϕeq cthe critical-state equilibrium compacity, that corresponds to the volume fraction obtained when a steady-state regime is reached (Schofield and Wroth 1968; Wood 1990). This critical-state compacity ϕeq cis generally a function of the solid pressure ps, of the shear rate ˙γ, and of the granular temperature, increasing with respect to ps. For the case of pores filled by fluid considered here, the granular temperature can be neglected. For steady granular flows it was shown in (GDR MiDi group 2004; Da Cruz et al. 2005; Cassar et al. 2005; Forterre and Pouliquen 2008) that ϕeq cis indeed a decreasing function of ˙γ/√psin the dry case (respectively of ˙γ/psin the wet case). This approach with critical-state compacity ϕeq callows to recover the different behaviours of loose and dense granular media, according to the sign of ϕ−ϕeq c. Namely, for a dense packing ϕ > ϕeq c, one has a positive dilation angle, ψ > 0, that induces dilation of the granular medium, ∇ · v > 0, while for a loose packing ϕ < ϕeq c, one has a negative dilation angle, ψ < 0, that induces contraction of the granular medium, ∇·v < 0. This is valid as soon as ˙γ > 0, i.e. when a deformation occurs. Pailha and Pouliquen (2009) deal with the immersed granular flows system. They consider the precedent model where a linearization of tan ψis proposed, tan ψ=K(ϕ−ϕeq c),(2.17) K > 0 being a calibration constant (dilation constant). We adopt this dilation model to write ∇·v=K˙γ(ϕ−ϕeq c).(2.18) Thus the closure considered in this work for (2.14) is Φ = K˙γ(ϕ−ϕeq c).(2.19) As exposed by Iverson (2005), Schaeffer and Iverson (2008), there is a coupling between the dilatancy and the pore pressure, called “pore pressure feedback”. This effect plays an important role in the way a landslide starts, and then dramatically affects the flow dynamics. The formula (2.18) well reproduces the contraction-dilation effects (see Andreotti et al. 2011; Pailha and Pouliquen 2009), which are •If ϕ > ϕeq cthen the granular medium dilates (∇·v > 0) as soon as there is a deformation (˙γ > 0). Consequently, –the fluid must be sucked into the mixture, –the pore pressure decreases. 7 •If ϕ < ϕeq cthen the granular medium contracts (∇·v < 0) as soon as there is a deformation (˙γ > 0). Consequently, –the fluid must be expelled from the mixture, –the pore pressure increases. The type of closure (2.18) entails a modification of the coefficient of the Coulomb friction law that becomes tan(δ+ψ) instead of tan δ. By linearization, we can write an effective friction coefficient as tan δeff = tan δ+ tan ψ. (2.20) In the thin-layer expansion performed below, we neglect the deviatoric solid stress e Tsinside the mixture, and only consider the bottom solid friction with the friction coefficient tan δeff. Closure laws slightly different from (2.19) are considered in Subsection 4.6. 2.4 Interpretation as a compressible model We would like here to propose an interpretation of the Roux and Radjai (1998) dilatancy relation under the form (2.18) as a compressible model, that enables to write down a fully dissipative energy equation in the case when the critical-state compacity ϕeq cdepends only on the pressure ps, and not on ˙γ. We consider the critical volume fraction ϕeq cto be an increasing function of the solid pressure only, ϕeq c=ϕeq c(ps), bounded by some maximal value ϕmax (ϕmax ∼0.6 for monodisperse spherical grains, but ϕmax can be higher for real polydisperse materials). This function ϕ= ϕeq c(ps) can be defined by its inverse p=peq c(ϕ) (peq c(ϕ) being called the critical pressure), as for example peq c(ϕ) = Kϕγ/(ϕmax −ϕ)ι, for some coefficient K, and some exponents γ,ι. Particular dependencies of peq c(ϕ) in ϕappear for example in (Lee et al. 2015). Since the granular temperature is negligible in the present context of pores filled by fluid, the critical pressure peq c(ϕ) is only related to the deformation of the grains that are in contact. A formula valid in the context of granular mixtures is given in (Iverson and George 2014) as peq c(ϕ)/p0= exp((ϕ−ϕmin)/a)−1, with 0.01 ≤a≤0.05 and 10 Pa ≤p0≤1000 Pa, depending on the materials. Classically in thermodynamics, the mechanical internal energy Uis related to the pressure pand volume Vby the relation dU=−pdV. Here the specific volume (i.e. volume per mass unit) is 1/(ρsϕ), thus to the critical pressure peq c(ϕ) one can associate by this relation a specific internal energy (i.e. internal energy per mass unit) eeq c(ϕ). Since d(1/ϕ) = −dϕ/ϕ2we obtain the differential relation deeq c dϕ=peq c ρsϕ2.(2.21) Then writing the mass equation (2.1a) as ∂tϕ+v· ∇ϕ+ϕ∇ · v= 0, and multiplying it by deeq c/dϕ, we get ∂teeq c+v·∇eeq c+peq c ρsϕ∇·v= 0.(2.22) Multiplying this by ϕand using again (2.1a) yields ∂t(ϕeeq c) + ∇·(ϕeeq cv) + peq c ρs∇·v= 0.(2.23) 8 Adding this times ρsto the energy equation (2.12) gives ∂tρsϕ|v|2 2+ρf(1 −ϕ)|u|2 2−(g·X)ρsϕ+ρf(1 −ϕ)+ρsϕeeq c +∇·ρsϕ|v|2 2v+ρf(1 −ϕ)|u|2 2u−(g·X)ρsϕv +ρf(1 −ϕ)u +pfmϕv + (1 −ϕ)u+e Tfmu+Tsv+ρsϕeeq cv = (ps−peq c)∇·v+e Ts:∇v+e Tfm:∇u+f·(v−u). (2.24) Now, according to (2.18) and since ps−peq c(ϕ) and ϕ−ϕeq c(ps) have opposite signs because ϕeq c is an increasing function of ps, one has (ps−peq c)∇ · v≤0, and the energy balance equation (2.24) has a nonpositive right-hand side. This means that, as required by the laws of physics, the total mechanical energy of the system is dissipated. Another way to understand the dilatancy law (2.18) is to perform a further linear approximation K(ϕ−ϕeq c(ps)) ≃Kp(peq c(ϕ)−ps), which is valid for ϕand psfar from extreme values 0, ϕmax and 0,∞respectively, with Kpthe order of magnitude of K(dpeq c/dϕ)−1. Then the dilatancy law (2.18) is tranformed into ∇·v=Kp˙γ(peq c(ϕ)−ps),(2.25) which can be written also as ps=peq c(ϕ)−∇·v Kp˙γ.(2.26) When (2.26) is introduced into the stress in (2.9a) it gives a diffusion equation on the solid velocity v(or on its divergence), with diffusion coefficient that is induced by the Roux-Radjai dilatancy law, DRR = dpeq c dϕ K˙γρsϕ.(2.27) Indeed (2.26) appears clearly as a compressible rheological law with bulk viscoplastic term (because of ˙γin the denominator), that can be compared with (6) in (Lee et al. 2015). We can propose also a general closure law under the form of a critical pressure peq c(ϕ, ˙γ) that generalises (2.25), by defining Φ in (2.14) as Φ = Kp˙γpeq c(ϕ, ˙γ)−ps, peq c(ϕ, ˙γ) = maxpcompr(ϕ),ηf˙γ Ieq(ϕ),(2.28) where pcompr(ϕ) is a static compressible law, and Ieq(ϕ) is an equilibrium relating the volume fraction ϕto the inertial number I. It can be for example Ieq(ϕ) = ( ¯ϕstat c−ϕ)/K2in the context of (5.8). The compressible pressure pcompr can be taken as above pcompr(ϕ)/p0= exp((ϕ− ϕmin)/a)−1. 2.5 Domain and boundary conditions We assume that the mixture (0 < ϕ < 1) lies between a fixed bottom and an interface, and that between the interface and an upper free surface, there is only fluid (ϕ≡0), see figure 1. The thickness of the mixture layer is denoted by hm, the thickness of the fluid-only layer by hf, and the fixed bottom is defined by a function b. The fluid velocity in the top layer is denoted by uf, and in the mixture layer by u, while vdenotes the velocity of the solid phase. For other terms, we will use as general notation the 9 About the stress tensors Tk(k=s, fm, f), they are decomposed as Tk=pkId + e Tk,(3.34) and suitable rheological assumptions should be made to define e Tk. A general approach has been proposed in (Bouchut and Boyaval 2016) to deal with velocity profiles in the thin-layer asymptotics and in the case of Newtonian or non-Newtonian rheologies. Here, as in (Bouchut and Westdickenberg 2004), since we aim to represent only depth-average effects, we prefer to simplify the rheologies and replace the effect of the stress tensors inside the domain by boundary layers due to the friction conditions, namely (3.10), (3.12), (3.22), and also due to the momentum conservation (3.18), while we neglect viscous effects. Thus we shall assume that the stresses e Tk are O(ǫ2) far from the boundaries z=b, b+hmand can just be nonzero close to these boundaries. Indeed, because of the particular form of (3.10), (3.12), (3.22), (3.18), we assume that e Txz s,e Txz fm,e Txz fcan be O(ǫ) close to the boundaries z=b, b +hm, but are O(ǫ2) far from these boundaries,(3.35) while the other components satisfy e Txx k=e Tzz k=O(ǫ2) everywhere.(3.36) Regarding the drag term defined in (2.7), we have according to (2.8) ˜ β=¯ β(t, x)1 + O(ǫ2),(3.37) with ¯ β= (1 −¯ϕ)2ηf ¯κ.(3.38) We shall consider two possible sets of assumptions. (i) The drag term is quite strong, that is ¯ β∼ǫ−1.(3.39) Then since the drag force ˜ β(u−v) has to balance gravity terms, it necessarily remains bounded. This implies that after an eventual initial layer (i.e. a short time interval during which the initial value of ux−vxis damped), one has ux−vx=O(ǫ).(3.40) (ii) The drag term is moderate, that is ¯ β=O(1).(3.41) In this case one has just ux−vx=O(1), according to (3.26). Note that in both cases one has ¯ β(ux−vx) = O(1). The relevance of the assumptions (3.39) or (3.41) can be evaluated as follows. According to (2.9a), the effective drag friction coefficient for the solid phase is ˜ β/ρsϕ. The assumption (3.39) or (3.41) has to be evaluated in the corresponding unit, which means that we must evaluate the dimensionless number ¯ βτ/ρs¯ϕ, with τ=pL/g the reference time unit (see above). We compute using (3.38) ¯ βτ ρs¯ϕ=(1 −¯ϕ)2 ¯ϕ ηfτ ρs¯κ.(3.42) 16 We consider the values g= 9.81m/s2,ρs= 2500kg/m3. In the typical laboratory experimental context we can take as in (Pailha and Pouliquen 2009) ¯κ=d2(1 −¯ϕ)3/(150 ¯ϕ2) with dthe diameter of the (spherical) grains. The typical values d= 5 ×10−4m, L≈0.2m, ¯ϕ/(1 −¯ϕ)≈1 with ηf= 10−3Pa s for water give a slightly strong dimensionless drag coefficient ¯ βτ/ρs¯ϕ≈34. For natural landslides or large scale USGS debris flows (Iverson et al. 2010), one can take as in Iverson and George (2014) a grain-size variability empirical formula ¯κ=κ0exp((0.6−¯ϕ)/0.04) with κ0≈10−11m2. We choose L≈20m, ηf= 10−2Pa s for muddy water, ¯ϕ≈0.5, which gives a very strong dimensionless drag coefficient ¯ βτ/ρs¯ϕ≈2×104. We conclude that the assumption (3.39) is valid in the natural context, while (3.41) is more valid in the experimental context. However, (3.41) could be valid also in the natural context if the permeability is higher κ0≈10−7m2for highly mobile flows (Iverson and George 2014). 3.4 Averaged momentum equations In order to get the averaged momentum equations, we have first to get expressions for the pressures. Computations shown in Appendix B give the fluid pressure in the fluid-only layer pf=ρfgcos θ(b+hm+hf−z) + O(ǫ2) for b+hm< z < b +hm+hf,(3.43) and in the mixture layer pfm=ρfgcos θ(b+hm+hf−z) + pe fm+O(ǫ2) for b < z < b +hm,(3.44) where pe fm≡¯ β 1−¯ϕZb+hm z (uz−vz)(z′)dz′(3.45) is the excess pore pressure. In the expression (3.44) of the fluid pressure we can see that there is an extra contribution pe fmto the commonly found hydrostatic pressure (3.43). A similar contribution to the hydrostatic pressure of the fluid phase is found in (Pailha and Pouliquen 2009). This excess pore pressure term is induced by the normal displacement produced by the dilation-compaction of the granular material immersed into the fluid. As seen on (3.45), the excess pore pressure is negative if the granular material goes up with respect to the fluid (vz> uz), and positive in the converse case. It vanishes at z=b+hm. The solid pressure is given (see Appendix B) by ps= ¯ϕ(ρs−ρf)gcos θ(b+hm−z)−pe fm+O(ǫ2) for b < z < b +hm.(3.46) Its nonhydrostatic component is the opposite of that of pfmin (3.44). About the averaged tangential components of momentum equations, we have the momentum equation for the fluid-only layer ρf(∂tux f+ux f·∇xux f) = −ρfgcos θ∇x(b+hm+hf)−1 hf1 2ρfVf+ki(ux f−ux) −ρfgsin θ(1,0)t+O(ǫ2), (3.47) the momentum equation for fluid phase in the mixture ρf(1 −¯ϕ)∂tux+ux·∇xux=−(1 −¯ϕ)ρfgcos θ∇x(b+hm+hf)−(1 −¯ϕ)∇xpe fm −1 hm1 2ρfVf−ki(ux f−ux) + kbux −¯ β(ux−vx)−(1 −¯ϕ)ρfgsin θ(1,0)t+O(ǫ2), (3.48) 17 where ∇xpe fmis given by (B.16), and the momentum equation for the solid phase ρs¯ϕ(∂tvx+vx·∇xvx) = −¯ϕg cos θρs∇xb+hm) + ρf∇xhf−(ρs−ρf)gcos θhm 2∇x¯ϕ + (1 −¯ϕ)∇xpe fm−sgn(vx) tan δeff ps|b hm +¯ β(ux−vx)−¯ϕρsgsin θ(1,0)t+O(ǫ2), (3.49) where according to (3.46) the bottom value of the solid pressure is given by ps|b= ¯ϕ(ρs−ρf)gcos θhm−(pe fm)|b+O(ǫ2),(3.50) and according to (3.45) (pe fm)|b=¯ β 1−¯ϕZb+hm b (uz−vz)(z′)dz′.(3.51) 3.5 Evaluation of the excess pore pressure The excess pore pressure pe fmis involved in (3.48), (3.49) and represents physically important effects. Thus it is necessary to derive an expansion of pe fmup to O(ǫ2) error terms. Recalling the definition (3.45) of pe fm, we have thus to evaluate uz−vzup to O(ǫ2) errors. We use equations (3.5) and (3.8) to get the solid normal velocity, vz=vx·∇xb+ (z−b)(¯ Φ−∇x·vx) + O(ǫ3).(3.52) Next, adding the mass equations in the mixture (3.2a), (3.2b), we find ∇x·(ϕvx+ (1 −ϕ)ux) + ∂z(ϕvz+ (1 −ϕ)uz) = 0,(3.53) and using (3.8) and (3.9), we get ϕvz+ (1 −ϕ)uz= ( ¯ϕvx+ (1 −¯ϕ)ux)·∇xb−(z−b)∇x·( ¯ϕvx+ (1 −¯ϕ)ux) + O(ǫ3).(3.54) Then, subtracting (3.52) to (3.54) yields uz−vz= (ux−vx)·∇xb−z−b 1−¯ϕ¯ Φ + ∇x·(1 −¯ϕ)(ux−vx)+O(ǫ3).(3.55) The definition (3.45) of pe fmthen gives for b < z < b +hm pe fm=¯ β 1−¯ϕ(b+hm−z)(ux−vx)·∇xb −1 2 h2 m−(z−b)2 1−¯ϕ¯ Φ + ∇x·(1 −¯ϕ)(ux−vx)+O(ǫ4). (3.56) Noticing that with either assumptions (i) or (ii) we have ¯ β=O(ǫ−1) (because a bounded term gives also something bounded when multiplied by ǫ), we deduce the bottom value (pe fm)|band the average pe fmas (3.58) and (3.59) below. We can then consider two possible sets of expansions for the values of (pe fm)|b,pe fm: (I) The values of (pe fm)|b,pe fmare given simply by (pe fm)|b=−¯ β (1 −¯ϕ)2 h2 m 2¯ Φ + O(ǫ2),pe fm=−¯ β (1 −¯ϕ)2 h2 m 3¯ Φ + O(ǫ2).(3.57) 18 expansion (I) expansion (II) assumption (i) i.e. ¯ β∼ǫ−1relevant relevant assumption (ii) i.e. ¯ β=O(1) not relevant relevant Table 1: Relevance of the formulas (I) or (II) for the values of (pe fm)|band pe fm (II) The values of (pe fm)|b,pe fmare given by (pe fm)|b=¯ β 1−¯ϕhm(ux−vx)·∇xb−h2 m 2(1 −¯ϕ)¯ Φ+∇x·(1−¯ϕ)(ux−vx)+O(ǫ3),(3.58) pe fm=¯ β 1−¯ϕhm 2(ux−vx)·∇xb−h2 m 3(1 −¯ϕ)¯ Φ+∇x·(1−¯ϕ)(ux−vx)+O(ǫ3).(3.59) Indeed, (3.57) follows from (3.58), (3.59) by dropping O(ǫ2) terms (because with either assumptions (i) or (ii) we have ¯ β(ux−vx) = O(1)). Thus the relations (I) are just simplified lower order approximations of the relations (II). However under assumption (ii) i.e. (3.41), it is not appropriate to consider (I) because the leading term is also O(ǫ2). Thus in this case only (II) is relevant, and the errors in (3.58), (3.59) are indeed O(ǫ4) as shown by the above computations. The relevance of the expansions (I) or (II) is summarized on table 1. We observe on (3.55) and (3.57) that at leading order, as explained in the introduction, the relative velocity uz−vzand the excess pore pressure pe fmhave sign opposite to ¯ Φ. 4 The two-phase two-layer model In the previous section we have established a complete set of equations for our two-phase twolayer model. In this section we give the main properties of this system. 4.1 System and first properties The system of equations derived in Section 3 has three scalar unknowns ¯ϕ,hm,hf, and three vector unknowns vx,ux,ux f. Dropping the error terms, it can be written as follows. The mass conservation equations follow from (3.23)-(3.25) by dropping O(ǫ3) terms, ∂t( ¯ϕhm) + ∇x·( ¯ϕhmvx) = 0,(4.1) ∂t(1 −¯ϕ)hm+∇x·(1 −¯ϕ)hmux=−Vf,(4.2) ∂thf+∇x·(hfux f) = Vf.(4.3) We can eliminate the fluid mass exchange term Vfby writing the fluid total mass conservation. Adding the two last equations yields ∂t(1 −¯ϕ)hm+hf+∇x·(1 −¯ϕ)hmux+hfux f= 0.(4.4) Adding (4.1) we deduce also whole system volume conservation as ∂t(hm+hf) + ∇x·(1 −¯ϕ)hmux+ ¯ϕhmvx+hfux f= 0.(4.5) 19 The evolution equation (3.33) for ¯ϕis ∂t¯ϕ+vx·∇x¯ϕ=−¯ϕ¯ Φ.(4.6) Multiplying it by hmand subtracting the result to (4.1), it yields ∂thm+∇x·(hmvx) = hm¯ Φ.(4.7) Finally, combining it with (4.5) gives ∂thf+∇x·(1 −¯ϕ)hm(ux−vx) + hfux f=−hm¯ Φ.(4.8) Thus, regarding scalar equations we have to keep a set of three independent equations for the three independent unknowns ¯ϕ,hm,hf. This can be either (4.1), (4.4), (4.6), or (4.1), (4.4), (4.8), or (4.6), (4.7), (4.8), or (4.1), (4.7), (4.8). This has to be completed by (4.2) or (4.3) to define Vf, that can in fact be expressed without time derivative, since subtracting (4.8) to (4.3) yields Vf=−hm¯ Φ−∇x·(1 −¯ϕ)hm(ux−vx).(4.9) The momentum equations are given by (3.47), (3.48), and (3.49). Thus the model is reduced to the following set of equations: ∂t( ¯ϕhm) + ∇x·( ¯ϕhmvx) = 0,(4.10a) ρs¯ϕ(∂tvx+vx·∇xvx) = −¯ϕg cos θρs∇x(b+hm) + ρf∇xhf −(ρs−ρf)gcos θhm 2∇x¯ϕ+ (1 −¯ϕ)∇xpe fm −sgn(vx) tan δeff ¯ϕ(ρs−ρf)gcos θhm−(pe fm)|b+ hm +¯ β(ux−vx)−¯ϕρsgsin θ(1,0)t,(4.10b) ∂t(1 −¯ϕ)hm+∇x·(1 −¯ϕ)hmux=−Vf,(4.11a) ρf(1 −¯ϕ)∂tux+ux·∇xux=−(1 −¯ϕ)ρfgcos θ∇x(b+hm+hf) −(1 −¯ϕ)∇xpe fm −1 hm1 2ρfVf−ki(ux f−ux) + kbux −¯ β(ux−vx)−(1 −¯ϕ)ρfgsin θ(1,0)t,(4.11b) ∂thf+∇x·(hfux f) = Vf,(4.12a) ρf(∂tux f+ux f·∇xux f) = −ρfgcos θ∇x(b+hm+hf) −1 hf1 2ρfVf+ki(ux f−ux)−ρfgsin θ(1,0)t,(4.12b) ∂t¯ϕ+vx·∇x¯ϕ=−¯ϕ¯ Φ,(4.13) where we used the formula ps|b= ¯ϕ(ρs−ρf)gcos θhm−(pe fm)|bfrom (3.50), the average ∇xpe fm is computed by (B.16) i.e. ∇xpe fm=1 hm∇x(hmpe fm) + (pe fm)|b∇xb,(4.14) 20 and according to (3.57) and (3.58), (3.59), (pe fm)|b=−¯ β (1 −¯ϕ)2 h2 m 2¯ Φ,pe fm=−¯ β (1 −¯ϕ)2 h2 m 3¯ Φ for case (I),(4.15) (pe fm)|b=−¯ β 1−¯ϕh2 m 2 ¯ Φ + ∇x·(1 −¯ϕ)(ux−vx) 1−¯ϕ−hm(ux−vx)·∇xb, pe fm=−¯ β 1−¯ϕh2 m 3 ¯ Φ + ∇x·(1 −¯ϕ)(ux−vx) 1−¯ϕ−hm 2(ux−vx)·∇xb       for case (II). (4.16) We put a positive part (we denote the positive part of a number xby x+≡max(0, x)) in the bottom solid friction term in (4.10b) because otherwise we could have a negative value for ps|b. The coefficient ¯ βis defined in (3.38) i.e. ¯ β= (1 −¯ϕ)2ηf/¯κ, and the closure function ¯ Φ is defined in (3.32) i.e. ¯ Φ = K¯ ˙γ( ¯ϕ−¯ϕeq c). We observe that writing the linear combination ρsvx×(4.10a)+hm×(4.10b) +ρfux×(4.11a)+hm×(4.11b) +ρfux f×(4.12a)+hf×(4.12b) we obtain the total momentum conservation ∂tρs¯ϕhmvx+ρf(1 −¯ϕ)hmux+ρfhfux f+∇x·ρs¯ϕhmvx⊗vx +ρf(1 −¯ϕ)hmux⊗ux+ρfhfux f⊗ux f+gcos θ∇x(ρs−ρf) ¯ϕh2 m 2+ρf (hm+hf)2 2 =−sgn(vx) tan δeff ¯ϕ(ρs−ρf)gcos θhm−(pe fm)|b+−kbux −ρs¯ϕhm+ρf((1 −¯ϕ)hm+hf)gcos θ∇xb+gsin θ(1,0)t. (4.17) The system (4.10)-(4.14) has the following other properties. It is a quasilinear system in case (I), while in case (II) it has an extra second-order term involving ∇x·(1 −¯ϕ)(ux−vx)due to the term ∇x(hmpe fm) in (4.14), and also a nonlinearity in terms of ∇x·(1 −¯ϕ)(ux−vx)in the bottom solid friction term. Next, solid and fluid masses are conserved, according to (4.10a) and (4.11a)+(4.12a). The width of the mixture hmremains nonnegative because of (4.10a). The solid volume fraction ¯ϕremains between 0 and 1 because of (4.6) and (3.32), indeed the value ¯ϕeq cis an attractive value for ¯ϕ. However, there is no reason for the width of the fluid-only layer hfto remain nonnegative, and this is due to the fact that the fluid could be fully sucked into the granular material. Therefore, our model is valid as long as hfremains nonnegative. Otherwise, one should write down equations that include the case of a mixture layer topped by a dry granular layer, what we have not done here. The system has the solution at rest characterized by vx=ux=ux f= 0, ¯ Φ = 0, b+˜ b+hm=cst,hf=cst, ¯ϕ=cst, with ˜ b≡xtan θ. 4.2 Comparison with other debris flows models In this subsection we would like to explain the main differences between our model and other debris flow models in the literature that include excess pore pressure effects, namely those of Pailha and Pouliquen (2009) and Iverson and George (2014). 4.2.1 The Pailha and Pouliquen model In Pailha and Pouliquen (2009) a two-phase debris flows model is proposed. As in our model, it is based on the dilatancy law proposed by Roux and Radjai (2.16). In their case the granular assembly is immersed, meaning that there is a thin mixture layer and a fluid layer above it, but as opposed to us the fluid layer is not thin but is approximately at rest. The hydrostatic pore 21 hm hf Figure 2: Immersed configuration. The pure fluid layer is not thin, it has a horizontal free surface. pressure satisfies ∇phydro f=ρfg, thus it is as if one would have a horizontal free surface, see figure 2. For us according to (3.43), (3.44) phydro f=ρfgcos θ(b+hm+hf−z) corresponds to the inclined free surface. Moreover they make the assumption of uniformity in the slope aligned direction x. This leads to a simplification regarding the normal velocities ϕvz+ (1 −ϕ)uz= 0, which is coherent with (3.54) only when there is no xdependency. Then they use the following relations for the solid velocity, vz=K4|vx|tan ψ, ˙γ=3|vx| h,(4.18) for some constant K4of the order of unity. This gives vz=1 3K4h˙γtan ψ, and this relation is indeed related to (3.52) where only the term in ¯ Φ is considered (according to the xindependency), with ¯ Φ = 2 3K4˙γtan ψ. This formula is indeed identical to (2.19), (2.17) with the choice K4= 3/2. Then, the solid pressure at the bottom is given by ps bed =ϕ(ρs−ρf)gh cos θ+K4 3 ηf κh2˙γtan ψ, (4.19) where the coefficients κand ηfare still the hydraulic permeability of the granular aggregate and the dynamic viscosity of the fluid, that are related to the drag friction by (2.8). Taking hm≡h in the value of ps|bin (3.50) with (4.15), i.e. case (I) (note that without xdependency, case (II) reduces to case (I)), we obtain the same equation (4.19). At the end we observe that in the x-independent case and choosing in our model hf(t, x) corresponding to the immersed situation with an horizontal free surface b+hm(t) + hf(t, x) + xtan θ=cst as on figure 2, our hydrostatic pressure satisfies ∇phydro f=ρfg(use (3.1)), and our model reduces to the one of Pailha and Pouliquen (with K4= 3/2), with a few differences. The first difference is that they consider dependency in ˙γarising from constitutive relations for shear stresses (that we have neglected in (3.35), (3.36)) in their dilatancy law and in their bottom friction. Depth-averaged models with internal shear stresses are also considered in (Gray and Edwards 2014). In our case we can nevertheless put the dependency in ˙γdirectly in the definition of ¯ϕeq cthat arises in the dilatancy closure law (as done in (4.52)), and in a bottom viscous friction term as in (5.4). The second difference is that they miss an interface fluid momentum exchange term. Since the pure fluid part remains at rest, we have to neglect the friction between this pure fluid part 22 and the mixture, ki= 0. Nevertheless we see that there remains in the two fluid momentum equations (4.11b) and (4.12b) a term proportional to the mass exchange Vfthat is computed by (4.9). The presence of this term is related to the fact that because of dilatancy there is a mass exchange between the two fluid parts, that induces also a momentum exchange. Since the total momentum is conserved including the static pure fluid layer, the fluid velocity in the mixture has to increase if the mixture fluid mass decreases, and vice-versa. The final effect of this exchange term is therefore the following. If tan ψ > 0 then the fluid is sucked into the mixture Vf<0 and the downslope fluid velocity in the mixture layer has to diminish. On the contrary if tan ψ < 0 then the fluid is expelled from the mixture Vf>0 and the downslope fluid velocity in the mixture layer has to increase. The third difference is that Pailha and Pouliquen slightly simplify their model by replacing the conservation of solid mass ¯ϕhmby the conservation of hm, and in several places ¯ϕby a critical static value ¯ϕstat c. These changes induce some slight differences in dynamical behaviour, see Section 5. We conclude that our model is an extension of that of (Pailha and Pouliquen 2009) to the case with dependency in the slope aligned variable x. It includes a coherent description of mass and momentum exchanges accross the interface, and importantly the average ∇xpe fmof the downslope gradient of the excess pore pressure, that is not active in (Pailha and Pouliquen 2009). 4.2.2 The Iverson and George model The dynamics of debris flows is described in (Iverson and George 2014) by a single phase model including dilatancy effects. An evolution equation for the fluid pore pressure is established using a dilatancy law and a Darcy law. Namely, these two relations are used: Dilatancy empirical law: ∇·v= ˙γtan ψ−αd dt(σ−pf),(4.20) Darcy law: (1 −ϕ)(u−v) = −κ ηf∇pe f,(4.21) with d/dt =∂t+v· ∇,αthe mixture compressibility, σ=ps+pfthe total normal stress and pf=ρfgcos θ(h−z) + pe fthe pore fluid pressure. The coefficients κand ηfare again the hydraulic permeability of the granular aggregate and the dynamic viscosity of the fluid. The Darcy law enables to express Φ = ∇ · v=∇ · ((1 −ϕ)(v−u)) = ∇ · (κ ηf∇pe f), as in (Morris and Boulay 1999; Lhuillier 2009; Nott et al. 2011). Combining both equations, the following evolution equation for the excess pore pressure is deduced, d dtpe f−1 α∇·κ ηf∇pe f=−˙γtan ψ α+d dtσ−ρfgcos θ(h−z).(4.22) Thus, the excess pore pressure pe fobeys a diffusion like equation with diffusion coefficient DIG = κ/ηfα, or with (3.38) DIG =(1 −¯ϕ)2 α¯ β.(4.23) This diffusion is analogous to the one acting on the solid velocity v(or on its divergence) in (2.27), associated to our compressible interpretation when (2.26) is assumed, and which is related to the unmodified Roux-Radjai dilatancy law. Note that according to the Darcy law one has ∇ · v=∇ · (κ ηf∇pe f), relating both diffusion equations. The diffusion coefficient (4.23) is proportional to 1/α instead of dpeq c/dϕ in (2.27). Both diffusions act in the whole 23 3D space variable. Valid ranges for the mixture compressibility αare discussed in (Andreini et al. 2013). This coefficient strongly affects the pore pressure diffusion in the sense that large values of αtranslate into delayed pore pressure diffusion. The value considered by Iverson and George (2014) is α= 5 ×10−5Pa−1, which corresponds to a typical value for a sediment-water mixture compressibility, while a typical value for pure water compressibility is of the order of 5×10−10 Pa−1(see Andreini et al. 2013; Montserrat et al. 2012). Note that for this value, the dilatancy law (4.20) is very close to that of Roux and Radjai. An evolution equation for the pore fluid pressure at the bottom pf bed can be deduced from the average of (4.22) knowing boundary conditions on σ, see (Iverson and George 2014) for details. In that procedure, the diffusion in xis neglected, and only the diffusion in zis taken into account. Then at the averaged level the diffusive aspect disappears. It can be seen that in the case α= 0 the equation on pf bed simplifies to pf bed =ρfgcos θh −1 2 ηf κh2˙γtan ψ. As for the Pailha and Pouliquen model above this corresponds again to the fluid pressure obtained with our formula (4.15) for (pe fm)|b, i.e. case (I), with hm≡h,hf= 0. The first main difference between our model and that of Iverson and George is that we have two free surfaces of respective heights hmand hm+hf, while they consider a single virtual free surface of height h. Therefore we have two equations (one for each free surface), while Iverson and George have only one. The height his defined by total mass conservation as ρmh=hmρsϕ+ρf(1 −ϕ)+ρfhf,with ρm=ρsϕ+ρf(1 −ϕ).(4.24) Many differences between the models follow from this fact. For example in (Iverson and George 2014) the hydrostatic fluid pressure is ρfgcos θ(h−z), while in our model it is ρfgcos θ(hm+ hf−z) (bis taken 0 here). This is due to the fact that pfvanishes at the fluid free surface, but not at the virtual free surface. Another consequence is that from the knowledge of ϕand total mass ρmh=ρmhm+ρfhf, it is not possible to compute hmnor hf. In particular the Iverson and George model does not describe separately the solid and fluid masses. Our detailed description with two free surfaces enables to describe accurately the mass and momentum exchanges between the two layers and the mass conservation for each phase. The second main difference between our model and that of Iverson and George is that they introduce the compressibility αin the dilatancy law as (4.20). Instead, we keep the original Roux- Radjai law (2.16), and introduce the compressibility in the critical state relation ϕ=ϕeq c(ps,˙γ) or equivalently ps=peq c(ϕ, ˙γ). Indeed we believe that a functional relation between ϕ,psand ˙γ can hold only in the critical state (i.e. when ∇·v= 0), and not for general states; preventing to write the fundamental relation (3.12) of (Iverson and George 2014). A possible general form for peq c(ϕ, ˙γ) is (2.28). The situation when ˙γ= 0 corresponds to quasistatic dynamics, that is related to soil mechanics. Compressibility values measured in such situation thus correspond to the law peq c(ϕ, ˙γ= 0), i.e. pcompr(ϕ) in (2.28). The last main difference is in the fact that Iverson and George write the Darcy law (4.21) relating the relative velocity v−uto the gradient ∇pe fof the excess pore pressure. According to the value (3.38) of ¯ β, the zcomponent of this Darcy law identifies indeed with our definition (3.45) of pe f. Then it is important to notice that Iverson and George neglect the xcomponent of the Darcy law (4.21), leading to the approximation ux≃vx. It follows that Iverson and George have as unknowns a single mixture velocity and pe f, while we have two velocities, or equivalently a mean mixture velocity and a relative velocity vx−ux. Thus in their case they replace the relative velocity by the excess pore pressure. In our model we have in case (II) all the terms involving the difference ux−vxin the evaluation of the excess pore pressure (4.16). In particular the term with ∇x·(1−¯ϕ)(ux−vx)induces a diffusion term in our system, in the slope aligned variable x. Its strength can be evaluated by writing the equation on (1 −¯ϕ)(ux−vx) that can 24 be deduced from (4.10b), (4.11b) (see also (4.37)). The result is a diffusion coefficient given by Dtwo−layer =¯ βh2 m 31 ρf(1 −¯ϕ)+1 ρs¯ϕ.(4.25) This diffusion arises at the level of the averaged system, and has no analogue in the Iverson and George model. Nevertheless it can be compared to the diffusion coefficient (4.23) of the non integrated Iverson and George model. In our model the diffusion equation is on (1−¯ϕ)(ux−vx), while in the Iverson and George model it is stated on pe f, both being somehow related as we said. However the two approaches differ by the physical interpretation of the diffusion. In our model the diffusion comes from the relative momentum equation, and from the expansion of the excess pore pressure (4.16) that appears at the averaged level. In the Iverson and George model, the Darcy law (4.21) means that the time derivative in the equation on the relative velocity ux−vx is neglected. It is reintroduced via the mixture compressibility αin (4.20). This difference results in different diffusion coefficients (4.25), (4.23). The proportionality to ¯ βin (4.25) is quite natural, it means that the largest the drag ¯ βis, the strongest is the diffusion coefficient and the damping of the relative velocity. A strong drag thus leads to a fast convergence to the hydrostatic equilibrium. The Darcy law (4.21) can be recovered from our model by an asymptotic expansion, see Subsubsection 4.4.2. It shows that an important correction to it, due to solid friction, has to be taken into account, and moreover that the slope aligned component cannot be neglected. If we look more precisely at conservation equations, we can describe the differences between the model proposed by Iverson and George (2014) and the one proposed in this work. Let us compare the continuity equations. The model that we propose contains two mass conservation equations, (4.10a), (4.11a)+(4.12a), and a closure equation (4.13), or equivalently (4.7) or (4.8). The Iverson-George model is defined by only two continuity equations ∂th+∇x·(hvx m) = ρm−ρf ρm D, (4.26) ∂t( ¯ϕh) + ∇x·( ¯ϕhvx m) = −ρf ρm ¯ϕD, (4.27) with ρm,hdefined in (4.24), vx ma mixture velocity, and D=Zh 0 (∇·v)dz. Indeed D=h¯ Φ with our notations. These continuity equations can be recast as ∂t(ρmh) + ∇x·(ρmhvx m) = 0,(4.28) ∂t¯ϕ+vx m·∇x¯ϕ=−¯ϕD h.(4.29) The second equation (4.29) is very similar to our closure equation (4.13), and the mass conservation (4.28) is similar to our total mass conservation ρs×(4.10a)+ρf×(4.11a)+ρf×(4.12a). The difference is indeed in the fact that in our case the total mass flux ρs¯ϕhmvx+ρf(1 −¯ϕ)hmux+ ρfhfux finvolves the three velocities vx,ux,ux f, and not only the velocity vxthat appears in the closure equation (4.13). Knowing the Darcy law (4.21) that gives the difference between vxand ux, we see that the Iverson-George model differs from ours. Regarding momentum equations, (George and Iverson 2014) have the equation ∂t(hvx m) + ∇x·(hvx m⊗vx m) + Kanis∇x(gcos θh2 2) + h(1 −Kanis) ρm∇xpf bed =−hg sin θ(1,0)t+ρm−ρf ρm Dvx m−1 ρmps bed tan δeff sgn(vx m) + 2ηf 1−¯ϕ hvx m, (4.30) 25 In the simplified case when the critical-state compacity depends only on the solid pressure ϕeq c=ϕeq c(ps) as described in Subsection 2.4, the most simple closure is to take ¯ϕeq c=ϕeq c(phydro s), phydro s=1 2(ρs−ρf) ¯ϕg cos θhm.(4.54) We shall call this closure closure (B). Still the energy dissipation in (4.47) is not always negative, because of the excess term appearing in the value of psin (4.48). Still in the case when ϕeq c=ϕeq c(ps), we can propose here another definition of ¯ Φ that is more consistent than (4.51) in terms of energy. For this, we put the same structure as that of the energy dissipation (2.24) of the 3D model in the energy dissipation Rin (4.47). This means defining ¯ϕeq c=ϕeq c(ps), where as in Subsection 2.4 the function ϕ=ϕeq c(p) is the inverse of the function p=peq c(ϕ). This definition leads automatically to a nonpositive first term in the formula (4.47) for R. However, (4.51) then means defining ¯ Φ as a function of ps, which by (4.48) and (4.15) or (4.16) itself depends on ¯ Φ. To avoid dealing with a nonlinear equation, we can use the closure (2.25) instead of (2.18), which means that we replace (4.51) by ¯ Φ = Kp¯ ˙γ(peq c( ¯ϕ)−ps).(4.55) Then with the relations (4.48) and (4.15) or (4.16), we obtain the value of ¯ Φ, 1 + ¯ βKp¯ ˙γ (1 −¯ϕ)2 h2 m 3¯ Φ = Kp¯ ˙γ peq c( ¯ϕ)−1 2(ρs−ρf) ¯ϕg cos θhm −¯ β 1−¯ϕh2 m 3(1 −¯ϕ)∇x·(1 −¯ϕ)(ux−vx)−hm 2(ux−vx)·∇xb!. (4.56) This formula is for case (II), i.e. (4.16). Otherwise for case (I) and (4.15) the second line must be removed in the right-hand side of (4.56). We shall call this closure closure (C). With this formula for case (II), we get R≤0, i.e. full dissipativity of the model. Note however that the coefficient ¯ βKp¯ ˙γh2 mhas to be not too large for the value of ¯ Φ from (4.56) to remain consistent, which is the case since ¯ β=O(ǫ−1). A discussion on the well-posedness of our models related to energy considerations is proposed in Appendix C. Concerning the value of ˙γ, if we consider that ˙γ=|Dv|, with Dv = (∇v+ (∇v)t)/2, the shear component ∇xvzsatisfies ∇xvz=O(ǫ), while Dxvx=O(1), ∂zvz=O(1) according to the asymptotic assumptions (3.26). Thus the leading term is ∂zvx, that can be of order ǫ−1(note that the technical assumption (3.27) is not relevant for real flows). Therefore |Dv| ≃ |∂zvx|and a typical value of ˙γis as (4.18), ¯ ˙γ=|vx|/hm. Then ¯ ˙γ∼1/ǫ, but this is compatible with the assumptions (3.26) as long as tan ψ=O(ǫ), which is satisfied for physically relevant data. Concerning the choice between Model (I), i.e. (4.15), or Model (II), i.e. (4.16), we recall that they have been derived under the assumptions (i) i.e. (3.39)-(3.40) or (ii) i.e. (3.41), according to table 1. Whatever are the scaling assumptions (i) or (ii), the closure (4.16) (thus Model (II)) is always more accurate, since in this case we retain all the terms in the expansions of Subsection 3.5. The advantage of Model (I) is that it is simpler since it involves only first-order derivatives. It can be used with closure (B) as (4.54). On the contrary, Model (II) involves second-order viscoplastic-like terms, and it is more natural for it to use closure (C) as (4.55), (4.56) since it enforces the full dissipation of energy. The variable ux−vxobeys an equation that includes a diffusion very similar to the one in the Iverson and George model, but which results from different 32 physical assumptions (see Subsection 4.2). In particular, the time derivative simply results here from the mass and momentum equations while it is related to mixture compressibility in Iverson and George. The sign of the excess pore pressure pe fmis always the same as that of the normal relative velocity uz−vzbecause of (3.45). According to (4.15) for Model (I), this sign is opposite to ¯ Φ, which corresponds exactly to the pore pressure feedback described in Subsection 2.3. However, other terms involving the tangential relative velocity ux−vxalso come into play in (4.16) for Model (II) to determine if the fluid is transferred into or out of the mixture. The dynamical behaviour of our model is naturally induced by the dilatancy closure of Roux and Radjai. Namely, in the absence of external inflow, there is convergence to the hydrostatic equilibrium over a sufficient long time. This can be seen in terms of the volume fraction by its evolution equation (4.6) that holds for any variant of the model considered. In this equation, the right-hand side ¯ Φ is defined according to one of the closures (A), (B) or (C). Thus we conclude that either ¯ ˙γtends to zero, or ¯ϕtends to ¯ϕeq c, which is an attractive value for (4.6). In any case ¯ Φ tends to zero. In the case of Model (II) with closure (C), the fully dissipative nature of the system, as seen on the right-hand side (4.47) of the energy equation, leads to normal dissipation (C.8) tending to zero and thus A,Bdefined in (4.45) also tending to zero. In any case we conclude that (pe fm)|b,pe fmin (4.15) or (4.16) tend to zero. This means that pe fmtends to zero, and the pressures becomes hydrostatic. Several limit systems can be obtained from our model for particular values of the parameters. The first is obtained by taking ¯ϕ≡0, leading simply to the standard shallow water system for height hm+hfand velocity ux. A second system is obtained for ¯ϕ≡1, leading to ¯ Φ≡0 and the usual two-layer shallow water system. A third system is obtained by taking ¯ Φ = 0, ¯ϕ=cst,¯ β→ ∞ which leads to ux−vx= 0, the excess pore pressure vanishes as well as the mass exchange, and there remains a two-layer system with fluid above with density ρf, width hf, velocity ux fand mixture below with density ρs¯ϕ+ρf(1 −¯ϕ), width hm, velocity vx. This system is indeed identical to the one of (Fern´andez-Nieto et al. 2008) with topography linearized around a constant slope and earth pressure coefficient K= 1. Finally, a fourth system is obtained by taking ρf= 0. It leads to ¯ β= 0, pe fm≡0, thus (I) and (II) are identical. The term ¯ Φ can be taken as in closures (A), (B) or (C). This yields an apparently new thin-layer model for a dry granular material with dilatancy effects, described with the unknowns ¯ϕ,hmand vx. A limitation of our Models (I) and (II) is that they have been derived under the assumption that there always remains some pure fluid above the mixture, hf>0. As we have said, this condition can naturally cease to be valid after a finite time. Then it would be necessary to consider negative hf, meaning that there is a layer of dry granular material above the mixture layer. Our approach can be applied to this situation indeed, then both models arising from both configurations could in principle be patched together to obtain a model valid in any situation where hfcan change sign according to space and time evolution. Such generalized model is highly relevant for practical studies, but requires additional analytical computations, it is postponed to future work. 5 Numerical tests for uniform flows In this section we perform numerical tests to compare our model to the models of (Pailha and Pouliquen 2009; George and Iverson 2014), in the very simple case of spatially uniform and 33 immersed configuration studied in (Pailha and Pouliquen 2009) and illustrated on figure 2. 5.1 Our model in the uniform immersed configuration To simulate underwater granular flows, we take the upper pure fluid layer at rest ux f= 0 in our three-velocity model (4.10)-(4.14), all quantities except hfare independent of x, and hm(t) + hf(t, x) + xtan θ=cst,ki= 0, according to Subsubsection 4.2.1. The ycomponents are taken zero. The equations are then ∂t( ¯ϕhm) = 0, ∂t¯ϕ=−¯ϕ¯ Φ,(5.1) ρs¯ϕ∂tvx=−sgn(vx)τb hm +¯ β(ux−vx)−¯ϕ(ρs−ρf)gsin θ, (5.2) ρf(1 −¯ϕ)∂tux=1 2ρfVf−kbux hm−¯ β(ux−vx),(5.3) with Vf=−hm¯ Φ, τb= tan δeff ps|b+K1ηf¯ ˙γ, (5.4) ps|b= ¯ϕ(ρs−ρf)gcos θhm−(pe fm)|b,(pe fm)|b=−¯ β (1 −¯ϕ)2 h2 m 2¯ Φ,(5.5) where the second term in the definition of τbin (5.4) is a viscous friction term. Alternatively, in order to really get total momentum conservation, we could put a coefficient 1 instead of 1/2 for the mass transfer term in (5.3) to compensate for neglecting the similar term Vfwhen setting ux f= 0 in the pure fluid momentum equation (4.12b). In order to compare with the experiments from (Pailha and Pouliquen 2009), we consider the same drag coefficient ¯ β= (1 −¯ϕ)2ηf ¯κ,¯κ=(1 −¯ϕ)3d2 150 ¯ϕ2,(5.6) with dthe diameter of the grains, and the same closure ¯ Φ ¯ Φ = ¯ ˙γtan ψ, (5.7) tan ψ=K( ¯ϕ−¯ϕeq c),¯ϕeq c= ¯ϕstat c−K2 ηf¯ ˙γ ps|b ,(5.8) and with ¯ ˙γ= 3|vx| hm ,tan δeff = tan δ+ tan ψ. (5.9) Note that (5.5), (5.7) (5.8) amounts to solve a quadratic equation in order to find the value of ps|b. We take here kb= 0. 5.2 The Pailha-Pouliquen model It is defined (in the case K4= 3/2, K3=K) as our model above, except that Vfin (5.3) is taken 0 instead of −hm¯ Φ, the mass equation ∂t( ¯ϕhm) = 0 is simplified to ∂thm= 0, and ¯ϕis simplified to ¯ϕstat cin several places, namely on the left-hand side of the momentum equations (5.2), (5.3), in the hydrostatic gravity term in the definition of the pressure ps|bin (5.5), and in the last gravity term of (5.2). Note for this gravity term that our sign convention for θis opposite to that of (Pailha and Pouliquen 2009). We take here θ < 0. 34 5.3 Iverson-George type model in the uniform immersed configuration We consider an Iverson-George type model obtained from our uniform immersed system (5.1)- (5.3) by adding the two momentum equations and neglecting the difference ux−vx, and also by considering a compressible coefficient αthat induces a relaxation time on the excess pore pressure. We obtain the system ∂t( ¯ϕhm) = 0, ∂t¯ϕ=−¯ϕ¯ Φ,(5.10) ρs¯ϕ+ρf(1 −¯ϕ)∂tvx=−sgn(vx)τb hm +1 2ρfVf−kbvx hm−¯ϕ(ρs−ρf)gsin θ, (5.11) with still Vf=−hm¯ Φ, τb= tan δeff ps|b+K1ηf¯ ˙γ, (5.12) ps|b= ¯ϕ(ρs−ρf)gcos θhm−(pe fm)|b,(pe fm)|b=−¯ β (1 −¯ϕ)2 h2 m 2¯ Φ,(5.13) but now instead of (5.7), ps|bsolves α∂tps|b=¯ ˙γtan ψ−¯ Φ.(5.14) The relations (5.6), (5.8), (5.9) remain unchanged. The equations (5.13), (5.14) can also be recast as ¯ Φ = −2(1 −¯ϕ)2 ¯ βh2 m¯ϕ(ρs−ρf)gcos θhm−ps|b, α∂tps|b=¯ ˙γtan ψ+ 2(1 −¯ϕ)2 ¯ βh2 m¯ϕ(ρs−ρf)gcos θhm−ps|b. (5.15) 5.4 Time scales involved in the different models Two timescales appear in all the three above models, tcrit = 1/¯ ˙γ, tvisc =ρsh2 m/3K1ηf,(5.16) that are related to the volume fraction equation and the convergence to a critical state, and to the viscous solid friction, respectively. Our model and the Pailha-Pouliquen model also contain the timescale trel involved in (5.3), trel =ρs/¯ β, (5.17) that governs the relaxation of the relative velocity ux−vx→0 as far as kb= 0 and Vfis sufficiently small. These timescales are indeed related to the dimensionless numbers Fr2and Sintroduced in (Pailha and Pouliquen 2009, p. 126) as F r2∼tvisc/tcrit,S∼tvisc/trel. The timescale trel is not present in the Iverson-George model since it involves a single velocity. In the Iverson-George type model, equation (5.15) involves another characteristic timescale tIG =α¯ βh2 m 2(1 −¯ϕ)2.(5.18) The form of (5.15) implies that after a few tIG the relation (5.7) will be satisfied, thus reducing to our model (5.1)-(5.9) except the fact that there is only one velocity. Differences between the Iverson-George solution and ours in the present uniform setting therefore only exist for small times less than a few tIG or a few trel. 35 The previous analysis of timescales may however be modified by the coupling terms that imply a feedback of one component to the other (to be rigorous, one should compute the eigenvalues of the jacobian matrix of the system, giving the inverse of the characteristic timescales; here we only have estimated the diagonal terms). In particular there is a coupling between the momentum equations (5.2), (5.3). There may be also a time delay depending on initial data, that can be necessary in order to approach a neighbourhood of the steady solution. This delay can be quite long since it involves nonlinear effects. 5.5 Numerical results We perform the tests proposed by Pailha and Pouliquen for low and high viscosity (figures 2 and 3 respectively of (Pailha and Pouliquen 2009)). The initial conditions are ux(t= 0) = 0 m/s, vx(t= 0) = 0 m/s,¯ϕ(t= 0) = ¯ϕ0, hm(t= 0) = h0 m.(5.19) For the Iverson-George model, ps|b(t= 0) is taken hydrostatic. The common data for the tests are ρs= 2500 kg/m3,¯ϕstat c= 0.582,tan δ= 0.415, d = 160µm, K= 4.09, K1= 90.5, K2= 25,(5.20) and the specific ones are •Low viscosity: ηf= 9.8×10−3Pa ·s, ρf= 1026 kg/m3,|θ|= 28o,h0 m= 6.1 mm. For the loose and dense cases, we choose the values ¯ϕ0= 0.576 and ¯ϕ0= 0.592, respectively. •High viscosity: ηf= 96 ×10−3Pa ·s, ρf= 1041 kg/m3,|θ|= 25o,h0 m= 4.9 mm. For the loose and dense cases, we choose the values ¯ϕ0= 0.562 and ¯ϕ0= 0.588, respectively. We solve the three above uniform models, that we call below −“Proposed simplified model” (proposed), given by equations (5.1)-(5.9). −“Pailha-Pouliquen model” (PP), described in Subsection 5.2. −“Iverson-George type model” (IG), given by equations (5.10)-(5.14), (5.6), (5.8), (5.9). In (Pailha and Pouliquen 2009) the simulations are shown up to final times corresponding to the duration of measurements in the laboratory experiments. Here we compute the solution until the steady state is reached. The steady state attained in large time can be computed explicitly for the three models. One can checck that it is given by the following formulas valid for |θ|> δ, ¯ϕ∞= ¯ϕstat c−K2 K1|tan θ|−tan δ, ηf¯ ˙γ∞ ps∞ |b =1 K1|tan θ|−tan δ. (5.21) Note that equation (5.21) tells that the same viscous inertial number is reached at steady-state for all models. These are completed by the steady thickness, pressure and velocity that may differ for each model, namely •Steady states for the Pailha-Pouliquen model: h∞ m=h0 m, ps∞ |b= ¯ϕstat ch∞ m(ρs−ρf)gcos θ, vx∞=1 3h∞ m¯ ˙γ∞. (5.22) 36 •Steady states for the simplified proposed model and for the Iverson-George type model: h∞ m=¯ϕ0 ¯ϕ∞h0 m, ps∞ |b= ¯ϕ∞h∞ m(ρs−ρf)gcos θ, vx∞=1 3h∞ m¯ ˙γ∞. (5.23) We have for the three models ux∞=vx∞, since kb= 0. This is due to the fact that here in the immersed configuration the gravity force term does not appear in the fluid equation (5.3), because it is balanced by the gradient of the upper fluid layer thickness hf. On the contrary, in the thin fluid layer case the gravity term is present, and (4.39) shows that the relative velocity does not become very small since there ∇xhfcannot balance the friction term. There is one main difference between the above two sets of steady states: the simplification of the mass conservation equation made by Pailha-Pouliquen leads to a constant height hmequal to the initial value, while in our case the conservation of solid mass ¯ϕhmgives a steady height depending on the initial height but also on the ratio between the initial and the steady volume fractions. It implies that in the Pailha-Pouliquen model the steady pressure and velocity depend only on the initial height, while in our model they depend also on the initial volume fraction (indeed they depend only on the initial solid mass). As a consequence, for the Pailha-Pouliquen model the same steady pressure and velocity are obtained for dense and loose initial volume fraction, as far as the initial height is fixed. On the contrary, different steady values are obtained for the proposed model. However, if the initial mass were fixed, varying initial volume fractions would yield the same steady pressure and velocity in our model, but different ones for the Pailha-Pouliquen model. 5.5.1 Comparison of our model to the Pailha-Pouliquen model We show on figure 3 the solutions vx,ps|band ¯ϕfor both models in loose and dense configurations for the high viscosity case, together with the steady states. Figure 3d shows the total solid mass hm¯ϕ. According to the formulas above, the unique steady state for the solid volume fraction ¯ϕis reached for both models in all cases, while different limits are obtained for velocity, pressure and height. The fact that our model ensures the conservation of solid mass allows to see differences between the loose and dense behaviours. With any of the two models we observe a sharp initial transition during which some components of the solution vary rapidly. The orders of magnitude of characteristic times are written on table 2. The characteristic time tcrit is the largest one, and figure 3 shows that there is a delay before that vxbecomes sufficiently far from zero so that convergence to the steady state can be achieved. The characteristic time trel is much smaller than the others, and the insert on the first figure shows that the relative velocity is quickly damped. Note however that the time scale related to this relative velocity looks rather correlated to tvisc, which indicates the feedback by coupling terms in (5.2), (5.3) mentioned above. Figure 4 helps in understanding the dilation/contraction behaviour. It shows the volume fraction ¯ϕas well as the critical volume fraction ¯ϕeq c, the height hmand the mass transfer term Vf. Dilation occurs when ¯ϕ > ¯ϕeq c, inducing a positive dilation angle tan ψ > 0. Contraction is found in the opposite case ¯ϕ < ¯ϕeq cwhich gives a negative dilation angle tan ψ < 0. For the loose initial configuration ¯ϕ0<( ¯ϕeq c)0= ¯ϕstat cthere is contraction (note that the lines corresponding to ¯ϕeq care above the curves representing ¯ϕ) and the fluid is expelled from the mixture. It implies a reduction of the height and a positive Vf. In the opposite case, for the dense initial configuration ¯ϕ0>¯ϕstat cthe solid assembly dilates (note that the lines corresponding to ¯ϕeq c 37 tcrit tvisc trel low viscosity 0.18 s to ∞3.7×10−2s 5.4×10−5s high viscosity 4.7 s to ∞2.3×10−3s 5.5×10−6s Table 2: Orders of magnitude of the characteristic times defined by (5.16), (5.17) for low or high viscosity data sets. The values of tcrit vary because they are inversely proportional to vx, that vanishes initially and tends to vx∞for large time. are under the curves representing ¯ϕ), and the fluid is sucked into the mixture. It implies an increase of the height and a negative Vf. The mass transfer Vfremains extremely small, less than 0.5×10−3mm/s. Our model and the Pailha-Pouliquen model give extremely close results, because the mass transfer Vfremains very small, as well as the relative velocity ux−vx. Results obtained with low viscosity are shown on figures 5 and 6. In this case we do not have clear loose/contraction and dense/dilation situations as for the high viscosity case. We can see on figure 5 that the solutions for both loose and dense initial mixtures look as having the same behaviour. Indeed the initially dense configuration behaves normally as in the high viscosity case, but the initially loose configuration very quickly changes to a dense state. If we look carefully at the value of the volume fraction on figure 6b we can see that for the loose case we find contraction until t= 2.7×10−3s (critical state lines are above solutions lines). At that time the critical state is passed over, and later on the mixture begins to dilate (critical state lines are under solutions lines). After a very short initial duration both initially loose or dense configurations lead to a dense/dilating configuration. The mass transfer term Vfbecomes negative and the height increases. Interestingly, the laboratory experiments plotted on figure 9b show the same qualitative behaviour. The convergence to the steady state occurs faster than in the high viscosity case, because the characteristic time tcrit is smaller, see table 2. Our model and the Pailha-Pouliquen model give again extremely close results, the mass transfer Vfremains very small (10−2mm/s), as well as the relative velocity ux−vx. 5.5.2 Comparison to the Iverson-George type model For the comparison to the Iverson-George type model we consider the two values of the compressibility αdiscussed in Subsubsection 4.2.2, namely 5 ×10−10 Pa−1and 5 ×10−5Pa−1, in order to evaluate the influence of the characteristic timescale tIG defined in (5.18). We only consider in details the case of high viscosity with ¯ϕ0= 0.562, while the other cases are shown in Subsubection 5.5.3. The characteristic time scales tcrit,tvisc,trel are given in table 2. The coefficient involved in the formula (5.18) of tIG takes the value ¯ βh2 m/(2(1 −¯ϕ)2) = 3 ×104Pa ·s for high viscosity, thus tIG = 1.5×10−5s for α= 5 ×10−10 Pa−1, tIG = 1.5 s for α= 5 ×10−5Pa−1(5.24) (for low viscosity, the values are respectively tIG = 2 ×10−6s and tIG = 2 ×10−1s). For the first value α= 5 ×10−10 Pa−1,tIG is even smaller than tvisc, thus the relaxation to the relation (5.7) is extremely fast and we do not observe a difference with our model (figure 7c,d). For the second value α= 5 ×10−5Pa−1,tIG is still much less than tcrit which is large initially. If follows that still there is no significant difference with our model (figure 7a,b). 38 0 200 400 600 800 1000 1200 1400 0 0.1 0.2 0.3 0.4 0.5 t(s) vx(mm/s) Loose: ¯ϕ0= 0.562 proposed PP Dense: ¯ϕ0= 0.588 proposed PP Steady states proposed (Loose) proposed (Dense) PP 012 x 10−3 0 0.5 1 1.5 x 10−3 (a) 0 200 400 600 800 1000 1200 1400 32 34 36 38 40 42 t(s) ps|b(Pa) (b) 0 200 400 600 800 1000 1200 1400 0.565 0.57 0.575 0.58 0.585 0.59 0.595 t(s) ¯ϕ (c) 0 200 400 600 800 1000 1200 1400 2.75 2.8 2.85 2.9 t(s) hm¯ϕ(mm) (d) Figure 3: High viscosity case: simulated (a) solid velocity vx, (b) solid pressure ps|b, (c) volume fraction ¯ϕ, and (d) hm¯ϕ. Horizontal dotted and dash-dotted lines represent the limit steady states. The relative velocity vx−uxis plotted in the insert of figure (a) (also in mm/s). 39 0 200 400 600 800 1000 1200 1400 0.565 0.57 0.575 0.58 0.585 0.59 0.595 0.6 0.605 0.61 t(s) ¯ϕand ¯ϕeq c Loose: ¯ϕ0= 0.562 ¯ϕproposed ¯ϕPP ¯ϕeq cproposed ¯ϕeq cPP Dense: ¯ϕ0= 0.588 ¯ϕproposed ¯ϕPP ¯ϕeq cproposed ¯ϕeq cPP (a) 0 0.2 0.4 0.6 0.8 1 x 10−3 0.565 0.57 0.575 0.58 0.585 0.59 0.595 0.6 0.605 0.61 t(s) ¯ϕand ¯ϕeq c (b) 0 200 400 600 800 1000 1200 1400 4.85 4.9 4.95 5 5.05 5.1 t(s) hm(mm) (c) 0 200 400 600 800 1000 1200 1400 −1 −0.5 0 0.5 1x 10−3 t(s) Vf(mm/s) (d) Figure 4: High viscosity case: simulated (a-b) volume fraction ¯ϕand critical volume fraction ¯ϕeq cgiven by (5.8), (b) represents a zoom of figure 4a at the very beginning of the simulation t∈[0,0.001s], (c) mixture thickness hm, and (d) mass transfer Vf. 40 0 50 100 150 200 0 2 4 6 8 10 12 14 t(s) vx(mm/s) Loose: ¯ϕ0= 0.576 proposed PP Dense: ¯ϕ0= 0.592 proposed PP Steady states proposed (Loose) proposed (Dense) PP 0 1 2 x 10−3 0 0.5 1 1.5 x 10−3 (a) 0 50 100 150 200 40 45 50 55 60 t(s) ps|b(Pa) (b) 0 50 100 150 200 0.54 0.55 0.56 0.57 0.58 0.59 0.6 0.61 t(s) ¯ϕ (c) 0 50 100 150 200 3.35 3.4 3.45 3.5 3.55 3.6 3.65 3.7 3.75 t(s) hm¯ϕ(mm) (d) Figure 5: Low viscosity case: simulated (a) solid velocity vx, (b) solid pressure ps|b, (c) volume fraction ¯ϕ, and (d) hm¯ϕ. Horizontal dotted and dash-dotted lines represent the limit steady states. The relative velocity vx−uxis plotted in the insert of figure (a) (also in mm/s). 41 Therefore, taking the scalar product of (2.38) with ufand subtracting the result to (A.4) yields ρf|u|2 2−|uf|2 2(˜ Nt+uf·˜ NX)−pfm(u−uf)·˜ NX −p∗ sVf−ρfVf(u−uf)·uf+Tfm˜ NX·(u−uf)≥0. Using the formula for Vfin (2.37), this can be written ρf 2|u−uf|2−p∗ sVf+(Tfm−pfmId) ˜ NX·(u−uf)≥0.(A.7) Next, we write the tangential part of (2.38) using (A.5), that gives (Tfm−Tf)˜ NXτ=ρfVf(uf−u)τ.(A.8) Together with (2.40) it yields Tfm˜ NXτ=1 2ρfVf(uf−u)τ−ki(uf−u)τ, Tf˜ NXτ=−1 2ρfVf(uf−u)τ−ki(uf−u)τ, (A.9) or equivalently Tfm˜ NX= (Tfm˜ NX)·˜ NX ˜ NX |˜ NX|2+1 2ρfVf(uf−u)τ−ki(uf−u)τ, Tf˜ NX= (Tf˜ NX)·˜ NX ˜ NX |˜ NX|2−1 2ρfVf(uf−u)τ−ki(uf−u)τ. (A.10) Plugging this in (A.7) we get the energy dissipation condition ρf 2(u−uf)·˜ NX |˜ NX|2 −p∗ s!Vf+(Tfm˜ NX)·˜ NX |˜ NX|2−pfm˜ NX·(u−uf) +ki|(uf−u)τ|2≥0. (A.11) Now, the term proportional to ki≥0 is a dissipation. The other terms are proportional to Vf because of (A.3). Since Vfcan be positive or negative, we write that its factor vanishes, that is ρf 2(u−uf)·˜ NX |˜ NX|2 −p∗ s+(Tfm˜ NX)·˜ NX |˜ NX|2−pfmϕ∗ 1−ϕ∗= 0.(A.12) This gives the value of p∗ s, and reporting this in (A.5) we finally obtain the interface condition (2.39). B Asymptotic expansion of the pressures and averaged momentum equations In this Appendix B we give the details of the calculations that justify the formulas of Subsection 3.4. 48 •Asymptotic expansion of the pressures. For the fluid-only layer we integrate the normal momentum equation (3.7b) with respect to z and use (3.14), (3.36) to get for b+hm< z < b +hm+hf pf=Tzz f+O(ǫ2) = ρfgcos θ(b+hm+hf−z) + O(ǫ2),(B.1) establishing (3.43). In the mixture, the normal fluid momentum equation (3.4b) gives with (3.36) ∂zpfm=−ρfgcos θ−¯ β 1−¯ϕ(uz−vz) + O(ǫ).(B.2) Integrating with respect to z, we obtain for b < z < b +hm pfm=pfm|b+hm+ρfgcos θ(b+hm−z) + ¯ β 1−¯ϕZb+hm z (uz−vz)(z′)dz′+O(ǫ2),(B.3) where the notation |b+hmmeans that the quantity is evaluated at z=b+hm. From (3.19), we have pfm|b+hm=pf|b+hm−ps|b+hm+O(ǫ2). Also from (3.20b) we have ps|b+hm=p∗ s+O(ǫ2), with according to (3.21), p∗ s=O(ǫ2). Thus p∗ s=O(ǫ2), ps|b+hm=O(ǫ2), pfm|b+hm=pf|b+hm+O(ǫ2).(B.4) Then from (B.1) we obtain the pressure for the fluid in the mixture at the interface, pfm|b+hm=ρfgcos θhf+O(ǫ2).(B.5) Finally with (B.3) we deduce the fluid pressure for the mixture layer (3.44), (3.45). Then, the solid normal momentum equation (3.3b) gives ∂zps=−¯ϕ∂zpfm−¯ϕρsgcos θ+¯ β(uz−vz) + O(ǫ).(B.6) Integrating with respect to zgives the expression of the solid pressure, ps=ps|b+hm−¯ϕ(pfm−pfm|b+hm)+ ¯ϕρsgcos θ(b+hm−z)−¯ βZb+hm z (uz−vz)(z′)dz′+O(ǫ2).(B.7) Using (B.4), (B.5), (3.44) and the notation (3.45), we finally obtain (3.46). •Tangential components of averaged momentum equations. For the fluid-only layer (3.7a), taking into account (B.1), we have ρf(∂tux f+ux f·∇xux f) = −ρfgcos θ∇x(b+hm+hf)−∂zTxz f−ρfgsin θ(1,0)t+O(ǫ2).(B.8) Next we write the mixture tangential fluid momentum equation (3.4a), using (3.44), ρf(1 −¯ϕ)∂tux+ux·∇xux=−(1 −¯ϕ)ρfgcos θ∇x(b+hm+hf)−(1 −¯ϕ)∇xpe fm −∂zTf xz m−¯ β(ux−vx)−(1 −¯ϕ)ρfgsin θ(1,0)t+O(ǫ2).(B.9) Similarly, the tangential solid momentum equation (3.3a) gives with (3.46) ρs¯ϕ(∂tvx+vx·∇xvx) = −(ρs−ρf)gcos θ∇x¯ϕ(b+hm−z)+ (1 −¯ϕ)∇xpe fm −¯ϕρfgcos θ∇x(b+hm+hf)−∂zTxz s +¯ β(ux−vx)−¯ϕρsgsin θ(1,0)t+O(ǫ2). (B.10) 49 We are going now to average (B.8) over the fluid layer, and (B.9), (B.10) over the mixture layer, so that the effects of the rheology are only taken into account by the boundary values of Txz f, Tf xz m,Txz s. According to (3.34), (3.36), the equation (3.14) gives pf=O(ǫ2) at the free surface, and then (3.13) yields Txz f=O(ǫ3) at z=b+hm+hf.(B.11) Next, using (B.4), the equation (3.20a) gives Txz s=O(ǫ3) at the interface, while (3.18) gives Txz f−Tf xz m=ρfVf(ux−ux f) + O(ǫ3) at the interface, where we recall that Vfis the fluid mass exchange between the layers from (2.37). But (3.22) gives Txz f+Tf xz m=−2ki(ux f−ux) + O(ǫ3) at the interface. We conclude that Txz s=O(ǫ3) at z=b+hm, Txz f=−(ki+1 2ρfVf)(ux f−ux) + O(ǫ3) at z=b+hm, Tf xz m=−(ki−1 2ρfVf)(ux f−ux) + O(ǫ3) at z=b+hm. (B.12) Finally, the conditions (3.10), (3.12) at the bottom give Txz s=−tan δeff vx |vx|Tzz s−2Txz s·∇xb+O(ǫ3) at z=b, Tf xz m=−kbux+O(ǫ3) at z=b. (B.13) Now, to go further, one would need an information on Tzz s−psup to O(ǫ3) error terms, i.e. on the rheology, and an expansion of psup to O(ǫ3), that we do not have in (3.46). This should lead to complementary terms in the expression of ps|b, as for example the term proportional to the curvature of the bottom and quadratic in vxas in (Bouchut and Westdickenberg 2004). We prefer here to avoid further technical expansions, and to simplify the presentation by dropping the term Txz s· ∇xbin (B.13), replacing Tzz sby ps, and considering that (3.46) is enough accurate. The additional terms that arise rigorously will eventually be written in future work. Alternatively, even if it is not physically relevant, an assumption of small solid friction tan δeff =O(ǫ) as in (Bouchut and Westdickenberg 2004) would justify avoiding extra terms. Then, when averaging the mixture momentum equations (B.9), (B.10), one comes up with the average excess pore pressure, that we can express with (3.45) as pe fm≡1 hmZb+hm b pe fm(z)dz =¯ β 1−¯ϕZb+hm b z′−b hm (uz−vz)(z′)dz′.(B.14) Then one computes Zb+hm b∇xpe fmdz =∇xZb+hm b pe fmdz −(pe fm)|b+hm∇x(b+hm) + (pe fm)|b∇xb. (B.15) Since (pe fm)|b+hm= 0, we deduce the expression of the average excess pore pressure force, ∇xpe fm≡1 hmZb+hm b∇xpe fmdz =1 hm∇x(hmpe fm) + (pe fm)|b∇xb.(B.16) We have also to average in (B.10) the term ∇x¯ϕ(b+hm−z), which gives 1 hmZb+hm b∇x¯ϕ(b+hm−z)dz =1 hmZb+hm b¯ϕ∇x(b+hm) + (b+hm−z)∇x¯ϕdz = ¯ϕ∇x(b+hm) + hm 2∇x¯ϕ. (B.17) 50 Therefore, averaging (B.8) over the fluid layer and using (B.11), (B.12), we obtain the momentum equation for the fluid-only layer (3.47). For the fluid phase in the mixture, averaging (B.9) and using (B.12), (B.13), we obtain (3.48). For the solid phase, averaging (B.10) with again (B.12), (B.13), and (B.17), we get (3.49). C Local energy identity This Appendix C is devoted to computations of local energy identities related to our models and stated in Subsection 4.5 as (4.41) and (4.46), and to comments on well-posedness related to these energy identities. •In order to establish (4.41), we proceed as in (Bouchut et al. 2015), by first performing the linear combination 1 2ρs|vx|2×(4.10a)+hmvx·(4.10b) +1 2ρf|ux|2×(4.11a)+hmux·(4.11b) +1 2ρf|ux f|2×(4.12a) +hfux f·(4.12b). Noticing that the terms in Vfcancel out, we obtain ∂tρs¯ϕhm|vx|2 2+ρf(1 −¯ϕ)hm|ux|2 2+ρfhf |ux f|2 2 +∇x·ρs¯ϕhm|vx|2 2vx+ρf(1 −¯ϕ)hm|ux|2 2ux+ρfhf |ux f|2 2ux f =−ρs¯ϕg cos θhmvx·∇x(b+˜ b+hm)−ρf¯ϕg cos θhmvx·∇xhf −(ρs−ρf)gcos θh2 m 2vx·∇x¯ϕ−ρfgcos θ(1 −¯ϕ)hmux+hfux f·∇x(b+˜ b+hm+hf) −(1 −¯ϕ)hm(ux−vx)·∇xpe fm−¯ βhm|ux−vx|2 −|vx|tan δeff ¯ϕ(ρs−ρf)gcos θhm−(pe fm)|b+−ki|ux f−ux|2−kb|ux|2. (C.1) The terms in ∇x(b+˜ b+hm) in the right-hand side of (C.1) are written using the mass equations (4.1), (4.4) as −gcos θρs¯ϕhmvx+ρf(1 −¯ϕ)hmux+hfux f·∇x(b+˜ b+hm) =−∇x·gcos θρs¯ϕhmvx+ρf(1 −¯ϕ)hmux+hfux f(b+˜ b+hm) +gcos θ(b+˜ b+hm)∇x·ρs¯ϕhmvx+ρf(1 −¯ϕ)hmux+hfux f =−∇x·gcos θρs¯ϕhmvx+ρf(1 −¯ϕ)hmux+hfux f(b+˜ b+hm) −gcos θ(b+˜ b+hm)∂tρs¯ϕhm+ρf(1 −¯ϕ)hm+hf =−∇x·gcos θρs¯ϕhmvx+ρf(1 −¯ϕ)hmux+hfux f(b+˜ b+hm) −gcos θ∂tρs¯ϕhm+ρf(1 −¯ϕ)hm+hf(b+˜ b) −ρfgcos θhm∂t(hm+hf)−(ρs−ρf)gcos θhm∂t( ¯ϕhm). (C.2) Similarly, the terms in ∇xhfin the right-hand side of (C.1) are written −ρfgcos θ¯ϕhmvx+ (1 −¯ϕ)hmux+hfux f·∇xhf =−∇x·ρfgcos θ¯ϕhmvx+ (1 −¯ϕ)hmux+hfux fhf −ρfgcos θhf∂t¯ϕhm+ (1 −¯ϕ)hm+hf. (C.3) 51 Then, the last term in the right-hand side of (C.2) is combined with the term in ∇x¯ϕin the right-hand side of (C.1), according to the identity −hm∂t( ¯ϕhm)−h2 m 2vx·∇x¯ϕ=−∂t¯ϕh2 m 2−h2 m 2∂t¯ϕ+vx·∇x¯ϕ,(C.4) where the right-hand side can be expressed with (4.6). Next, the excess pore pressure term in (C.1) can be written using (4.14) −(1 −¯ϕ)hm(ux−vx)·∇xpe fm =−∇x·(1 −¯ϕ)hmpe fm(ux−vx) +hmpe fm∇x·(1 −¯ϕ)(ux−vx)−(1 −¯ϕ)(pe fm)|b(ux−vx)·∇xb. (C.5) Finally, multiplying (4.6) by deeq c dϕ ( ¯ϕ) and using (2.21), we obtain ∂teeq c( ¯ϕ)+vx·∇xeeq c( ¯ϕ)=−peq c( ¯ϕ) ρs¯ϕ¯ Φ.(C.6) Multiplying this by hm¯ϕand using (4.10a), we deduce ∂thm¯ϕeeq c( ¯ϕ)+∇x·hm¯ϕeeq c( ¯ϕ)vx=−hm peq c( ¯ϕ) ρs ¯ Φ.(C.7) Using the formulas (C.2)-(C.5), (C.7) in (C.1), this yields the energy balance equation (4.41). •Next we prove the identity (4.46) for case (II). The value (4.44) of uz−vzgives Zb+hm b (uz−vz)2dz =h3 m 3A2−h2 mAB +hmB2.(C.8) Then, using the definition of pe fmin (4.16), −hmpe fm ¯ Φ−hmpe fm∇x·(1 −¯ϕ)(ux−vx)=−hmpe fm(1 −¯ϕ)A=¯ βhmAh2 m 3A−hm 2B.(C.9) With the definition of (pe fm)|bin (4.16) it yields −hmpe fm ¯ Φ−hmpe fm∇x·(1 −¯ϕ)(ux−vx)+ (1 −¯ϕ)B(pe fm)|b =¯ βhmAh2 m 3A−hm 2B−¯ βBh2 m 2A−hmB =¯ βh3 m 3A2−h2 mAB +hmB2=¯ βZb+hm b (uz−vz)2dz. (C.10) According to the definition (4.42) of Re, this can be written equivalently (4.46), proving the claimed identity. •Concerning well-posedness of the model, case (I) shows a quasilinear system with an energy identity that has a formally small right-hand side, but that contains derivatives because of the term ∇x·(1 −¯ϕ)(ux−vx)in Rein (4.42). Thus we do not obtain a mathematical entropy for the system. In case (II) the system contains second-order terms because of the term ∇xpe fmin (4.10b) and (4.11b), defined by (4.14) where pe fminvolves ∇x·(1−¯ϕ)(ux−vx)in (4.16). However, the right-hand side Rof the energy equation (4.41) contains (4.49) that can be expressed with (C.8) as a positive definite quadratic form in terms of A,Bdefined in (4.45). We therefore have bounds on A,B(especially in case of closure (C)) and thus also on pe fm,¯ Φ, and on ∇x·(1−¯ϕ)(ux−vx) since Ais bounded. Thus we have a hyperbolic/parabolic system for which we can expect wellposedness. Note that the term ∇x·(1 −¯ϕ)(ux−vx)is the depth-averaged counterpart of ∇·v in (2.24), (2.26). 52 References T. B. Anderson and R. Jackson. A fluid mechanical description of fluidized beds. Ind. Eng. Chem. Fundam., 6:527–539, 1967. N. Andreini, C. Ancey, and G. Epely-Chauvin. Granular suspension avalanches II. Plastic regime. Phys. Fluids, 25:033302, 2013. B. Andreotti, Y. Forterre, and O. Pouliquen. Les milieux granulaires. Physique Savoirs Actuels. EDP Sciences, 2011. G.S. Beavers and D.D. Joseph. Boundary conditions at a naturally permeable wall. J. Fluid Mech., 30:197–207, 1967. M.D. Bolton. The strength and dilatancy of sands. Geotechnique, 36:65–78, 1986. F. Bouchut and S. Boyaval. Unified derivation of thin-layer reduced models for shallow freesurface gravity flows of viscous fluids. European J. Mech.-B/Fluids, 55:116–131, 2016. F. Bouchut and M. Westdickenberg. Gravity driven shallow water models for arbitrary topography. Comm. in Math. Sci., 2:359–389, 2004. 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´erie I, 336:531–536, 2003. F. Bouchut, E. D. Fern´andez-Nieto, A. Mangeney, and G. Narbona-Reina. A two-phase shallow debris flow model with energy balance. ESAIM: Math. Modelling Num. Anal., 49:101–140, 2015. H. Brenner. Bi-velocity hydrodynamics. multicomponent fluids. Intern. J. Engin. Sci., 47: 902–929, 2009. H. Brenner. Diffuse volume transport in fluids. Phys. A, 389:4026–4045, 2010. C. Cassar, M. Nicolas, and O. Pouliquen. Submarine granular flows down inclined planes. Phys. Fluids, 17:103301, 2005. F. Da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir. Rheophysics of dense granular materials: Discrete simulation of plane shear flows. Phys. Rev. E, 72:021309, 2005. E. D. Fern´andez-Nieto, F. Bouchut, D. Bresch, M. J. Castro D´ı`az, and A. Mangeney. A new savage-hutter type model for submarine avalanches and generated tsunami. J. Comput. Phys., 227:7720–7754, 2008. Y. Forterre and O. Pouliquen. Flows of dense granular media. Annual Rev. Fluid Mech., 40: 1–24, 2008. GDR MiDi group. On dense granular flows. Eur. Phys. J. E, 14:341–365, 2004. 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. Italian Journal of Engineering Geology and Environment, pages 415–424, 2011. URL DOI: http://dx.doi.org/10.4408/IJEGE.2011-03.B-047. 53 D. L. George and R. M. Iverson. A depth-averaged debris-flow model that includes the effects of evolving dilatancy. II. Numerical predictions and experimental tests. Proc. R. Soc. A, 470: 20130820, 2014. J. M. N. T. Gray and A. N. Edwards. A depth-averaged µ(I)-rheology for shallow granular free-surface flows. J. Fluid Mech., 755:503–534, 2014. R. M. Iverson. The physics of debris flows. Rev. Geophys., 35:245–296, 1997. R. M. Iverson. Landslide triggering by rain infiltration. Water Resour. Res., 36:1897–1910, 2000. R. M. Iverson. Regulation of landslide motion by dilatancy and pore pressure feedback. J. Geophys. Res., 110:F02015, 2005. R. M. Iverson. Elements of an improved model of debris-flow motion. Powders and Grains, American Institute of Physics, Proceedings, 1145:9–16, 2009. R. M. Iverson and D. L. George. Modelling landslide liquefaction, mobility bifurcation and the dynamics of the 2014 oso disaster. Geotechnique, 66:175–187, 2016. R. M. Iverson and D. L. George. A depth-averaged debris-flow model that includes the effects of evolving dilatancy. I. Physical basis. Proc. R. Soc. A, 470:20130819, 2014. 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:F03005, 2010. R. Jackson. Some mathematical and physical aspects of continuum models for the motion of the granular materials, volume 714 of in Theory of Dispersed Multiphase Flow. R. E. Meyer, 1983. R. Jackson. The Dynamics of Fluidized Particles. Cambridges Monographs on Mechanics. Cambridge University Press, 2000. J. Kowalski and J. N. McElwaine. Shallow two-component gravity-driven flows with vertical variation. J. Fluid Mech., 714:434–462, 2013. C. H. Lee, C. J. Huang, and Y.M. Chiew. A three-dimensional continuum model incorporating static and kinetic effects for granular flows with applications to collapse of a two-dimensional granular column. Phys. Fluids, 27:113303, 2015. D. Lhuillier. Migration of rigid particles in non-brownian viscous suspensions. Phys. Fluids, 21: 023302, 2009. J. K. Mitchell. Fundamentals of soil behaviours. New York. Wiley, 1993. S. Montserrat, A. Tamburrino, O. Roche, and Y. Ni˜no. Pore fluid pressure diffusion in defluidizing granular columns. J. Geophys. Res., 117:F02034, 2012. T. Morales de Luna. A saint venant model for gravity driven shallow water flows with variable density and compressibility effects. Math. Comput. Modelling, 47:436–444, 2008. J. F. Morris and F. Boulay. Curvilinear flows of non-colloidal suspensions: the role of normal stresses. J. Rheol., 43:1213–1237, 1999. P. R. Nott, E. Guazzelli, and O. Pouliquen. The suspension balance model revisited. Phys. Fluids, 23:043304, 2011. 54 M. Pailha and O. Pouliquen. A two-phase flow description of the initiation of underwater granular avalanches. J. Fluid Mech., 633:115–135, 2009. M. Pailha, M. Nicolas, and O. Pouliquen. Initiation of underwater granular avalanches: influence of the initial volume fraction. Phys. Fluids, 20:111701, 2008. M. Pelanti, F. Bouchut, and A. Mangeney. A roe-type scheme for two-phase shallow granular flows over variable topography. ESAIM: Math. Modelling Num. Anal., 42:851–885, 2008. Y. Penel, S. Dellacherie, and B. Despr´es. Coupling strategies for compressible-low mach number flows. Math. Models Meth. Appl. Sci., 25:1045–1089, 2015. E.B. Pitman and L. Le. A two-fluid model for avalanche and debris flows. Phil.Trans. R. Soc. A, 363:1573–1601, 2005. O. Reynolds. On the dilatancy of media composed of rigid particles in contact. Phil. Mag. Ser. 5, 20:469–481, 1885. J. F. Richardson and W. N. Zaki. Sedimentation and fluidization: part I. Trans. Inst. Chem. Eng., 32:35–53, 1954. L. Rondon, O. Pouliquen, and P. Aussillous. Granular collapse in a fluid: role of the initial volume fraction. Phys. Fluids, 23:073301, 2011. S. Roux and F. Radjai. Texture-dependent rigid plastic behaviour. In H. J. Herrmann et al., editor, Physics of Dry Granular Media, pages 305–311. Kluwer, 1998. D. G. Schaeffer and R. Iverson. Steady and intermittent slipping in a model of landslide motion regulated by pore-pressure feedback. SIAM J. Appl. Math., 69:768–786, 2008. A. N. Schofield and C. P. Wroth. Critical State Soil Mechanics. McGraw-Hill Inc., 1968. I. Vardoulakis. Dynamic stability analysis of undrained simple shear on water-saturated granular soils. Int. J. Num. Anal. Mech. Geomechanics, 10:177–190, 1986. D. M. Wood. Soil behavior and critical state soil mechanics. Cambridge, UK. Cambridge University Press, 1990. 55