J Garres-Díaz

F Bouchut

E D Fernández-Nieto

A Mangeney
email: mangeney@ipgp.fr

G Narbona-Reina

Multilayer models for shallow two-phase debris flows with dilatancy effects

. In the numerical results, we observe that the proposed model with a two layer description of the mixture accurately represents the velocity measured at the surface of the mixture in the laboratory experiments. This is obviously poorly represented by the depth-averaged velocity in single-layer models while the other quantities (solid volume fraction, basal excess pore fluid pressure) are similar to those obtained with single-layer models. Our numerical results show a significant impact of the parameters involved in dilatancy law, in particular on the calculation of the time evolution of the excess pore fluid pressure.

Introduction

Many efforts have been devoted to the study of granular flows (aerial, sub-aerial, debris flows,...) in recent years. Gravity driven flows such as landslides, submarine avalanches, or rock avalanches are important natural hazards. One of the challenges of those studies is to predict the velocity and the runout distance in rapid landslides in order to provide new tools for prediction and prevention systems. The physical understanding of these flows, and their theoretical description is still a challenge from the geophysical, mathematical and numerical point of view. In last years, mathematical modelling has contributed to investigate granular flows. However, many questions remain open regarding the internal behavior of these flows, as the rheology and the fluid-grain interaction in fluidized flows.

On the one hand, rheological laws able to explain the complex behavior of granular flows have been widely investigated. Currently, the most accepted law is the so-called µ(I)-rheology, introduced in [START_REF] Jop | A constitutive law for dense granular flows[END_REF], and other equivalent laws that consider a variable friction coefficient in the definition of the stress tensor (see Capart et al. [START_REF] Capart | Depth-integrated equations for entraining granular flows in narrow channels[END_REF]). On the other hand, most of natural landslides involve a fluid (commonly water) mixed with the granular material. Interactions between the fluid and solid phases play a key role in the dynamics of these debris flows. This interaction basically depends on the pore fluid pressure, which determines the effective friction between the fluid and the grains. Many studies have been devoted to the study of dilatancy effects in granular materials. Compression/dilatation implies a decrease/increase of the fluid pressure that has a strong effects on the dynamics of fluidized flows (see e.g. [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF][START_REF] Bouchut | A two-phase solid-fluid model for dense granular flows including dilatancy effects: comparison with submarine granular collapse experiments[END_REF][START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF]).

Most of the models in the literature describing two-phase (grain-fluid) flows are based on the Jackson's model [START_REF] Jackson | The Dynamics of Fluidized Particles[END_REF], which considers the buoyancy force and the friction between the phases. This model has five unknowns: the solid (grain) volume fraction, the fluid and solid pressures, and solid and fluid velocities. However, only four equations are used: mass and momentum conservation equations for each phase. Therefore, the Jackson's model is underdetermined and a closure equation has to be added. It seems convenient to formulate this closure equation in term of the contraction/dilatation effect as discussed in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF].

As it is well known, the computational cost of solving 3D models is huge. For this reason, the thin-layer (i. e. shallow) approximation have been commonly used to obtain depth-averaged models from the Jackson's model. Pitman and Le [START_REF] Pitman | A two-fluid model for avalanche and debris flows[END_REF] and later Petanti et al. [START_REF] Pelanti | A Roe-type scheme for two-phase shallow granular flows over variable topography[END_REF] proposed a depth-averaged model where no closure relation was included but used an extra boundary condition instead. This approach makes it impossible to obtain a dissipative energy balance, leading therefore to a model that is physically meaningless. This is widely discussed in [START_REF] Bouchut | A two-phase shallow debris flow model with energy balance[END_REF], where a model with a consistent energy balance is proposed based on the incompressibility of the solid phase. However, the granular phase is not incompressible and dilatancy effects have to be accounted for.

Roux and Radjai introduced a dilatancy model in [START_REF] Roux | Texture-Dependent Rigid-Plastic Behavior[END_REF], namely the rate of the solid volume change is defined in terms of the shear rate of the granular material γ and a dilatation angle ψ as -γ tan ψ. Actually, they write the evolution of the solid volume fraction ϕ as

∂ t ϕ + v • ∇ϕ = -ϕ γ tan ψ,
where v is the solid velocity. Pailha and Pouliquen [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] used this relation as a closure of the Jackson's model to obtain a depth-averaged model, but they also assume an extra boundary condition and therefore their system is overdetermined. When the previous equation is combined with the mass conservation equation, it leads to

∇ • v = γ tan ψ, (1)
which is the closure relation considered by Bouchut et al. [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF][START_REF] Bouchut | A two-phase solid-fluid model for dense granular flows including dilatancy effects: comparison with submarine granular collapse experiments[END_REF]. A key point in that model is to consider a thin layer of fluid over the mixture layer. Then, the contraction/dilatation of the mixture is allowed through a fluid transference between the two layers. As in the Pailha and Pouliquen model [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF], the dilatancy effects appear in this model in particular through an excess pore fluid pressure term representing a deviation from the hydrostatic solid and fluid pressures.

A different approach to two-phase models is to assume a single-phase flow (e.g. Iverson [START_REF] Iverson | The physics of debris flows[END_REF], George and Iverson [START_REF] George | A two-phase debris-flow model that includes coupled evolution of volume fractions, granular dilatancy, and pore-fluid pressure[END_REF][START_REF] George | A depth-averaged debris-flow model that includes the effects of evolving dilatancy. II. Numerical predictions and experimental tests[END_REF], Iverson and George [START_REF] Iverson | A depth-averaged debris-flow model that includes the effects of evolving dilatancy. I. Physical basis[END_REF]) where the mass and velocity of the mixture, which are defined as an average in terms of the solid volume fraction, are used instead of the masses and velocities of each phase. Therefore, the relative motion between the fluid and solid is not explicit in the model. In [START_REF] Iverson | A depth-averaged debris-flow model that includes the effects of evolving dilatancy. I. Physical basis[END_REF] they consider dilatancy effects through a modification of the dilatancy law of Roux and Radjai [START_REF] Roux | Texture-Dependent Rigid-Plastic Behavior[END_REF] and took into account the compressibility of the granular material.

The previous depth-averaged models share an important limitation. The dilatancy effect is written in terms of the shear rate of the solid phase γ, which strongly depends on the variations in the direction normal to the topography. Therefore, these "single-layer" models give a poor approximation of this term. Moreover, contraction/dilatation produces a relative motion of the fluid and solid phases in the normal direction that these models hardly reproduce. Multilayer models were introduced as an intermediate step between depth-averaged and 3D Navier-Stokes models (e.g. [START_REF] Audusse | A multilayer Saint-Venant model: derivation and numerical validation[END_REF][START_REF] Fernández-Nieto | A Multilayer Method for the Hydrostatic Navier-Stokes Equations: A Particular Weak Solution[END_REF]). These models allow us to recover the normal structure of the flow, and consequently to get a better approximation of the terms that depend on the normal variations. Fernández-Nieto et al. [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF] introduced a multilayer model for dry granular flows with the µ(I) rheology. They showed that the rheological terms are better approximated thanks to this approach, making it possible to numerically reproduce relevant behaviour of dry granular flows. In particular, they recover the position and evolution of the flow/no-flow transition. Furthermore, their model were able to simulate naturally the change of velocity profiles from Bagnold to S-shaped. This was deeply investigated by Fernández-Nieto et al. [START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: a well-balanced multilayer discretization[END_REF], where side walls friction was added, changing significantly the flow dynamics.

The goal of this work is to propose a multilayer extension of model [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] for fluidized granular flows with dilatancy effects. To this aim, we start from the Jackson's model with the closure relation [START_REF] Abgrall | Two-Layer Shallow Water System: A Relaxation Approach[END_REF], and the final model is obtained from an asymptotic analysis and the multilayer approach. Then, the velocity and solid volume fraction is no more constant along the normal direction. Actually, both are piecewise constant functions, with discontinuities at the internal interfaces of the multilayer domain. One of the key point of this work is to highlight the main numerical difficulty in such multilayer (or 3D) models for debris flows: to deal with the strongly non-linear and abrupt changes of the excess pore fluid pressure when starting from non-equilibrium conditions. We propose an approximation of this excess pore fluid pressure in simple uniform flows and compare our numerical results with previous depth-averaged (i. e. single-layer) models and with analytical solutions to quantify the errors made in the proposed numerical approximation and how much adding layers in the normal direction improves the comparison with laboratory experiments compared to single-layer models.

The paper is organized as follows: in section 2 the governing equations and the rheology are presented. Section 3 is devoted to the details of the derivation of the multilayer model, from the dimensional analysis of the Jackson's model to the final two-phase multilayer model. The numerical tests are presented in section 4 and some conclusions are introduced in section 5. The explicit expression of the excess pore pressure is done in Appendix A. Appendix B collects all the details of the special case of the 2 layer model preserving the granular mass, named PGM-2 model. This model is the multilayer model with a smaller computational cost and which we compare with depth-averaged models. We show that it reproduces accurately experimental data, mainly the velocity measured at the surface of the mixture, as can be seen in section 4.

The initial system

In this section we present the governing equations together with the closure relations and the rheology describing fluidized granular flows, as well as the boundary conditions of the system. As in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF], the domain consists of a grain-fluid mixture layer and a fluid layer on its top. The role of this upper layer is to allow the exchange of fluid V f between the two layers (Figure 1) that is a key point on the dilatancy effect. This exchange is given through the kinematic boundary condition at the mixture interface.

Governing equations

As mentioned before, the goal of this paper is to derive a new multilayer model for fluidized granular flows that can be seen as an extension of the work developed in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF]. We focus on the modelling of the mixture layer and the model for the single upper layer is adopted from [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] since it does not change. Thus, we focus on the detailed modelling of the mixture layer in the following.

We consider a granular material and a fluid with constant densities ρ s , ρ f . The governing equations are given by the Jackson's model for fluidized granular materials. The system consists of the concatenation of the mass and the momentum equations for the granular and fluid phases. The equations for the solid phase are considered only in the mixture domain and the fluid equations are considered in both domains. We denote u f the velocity of the fluid in the top layer and v, u ∈ R 2 the solid and fluid velocities respectively in the mixture layer. Thus, for the upper layer the mass and momentum equations read:

∇ • u f = 0, ρ f (∂ t u f + (u f • ∇)u f ) = -∇ • σ f + ρ f g,
whereas in the mixture domain, the Jackson's complete system reads:

             ρ s ∂ t ϕ + ∇ • (ϕv) = 0, ρ f ∂ t (1 -ϕ) + ∇ • ((1 -ϕ) u) = 0, ρ s ϕ ∂ t v + (v • ∇) v = ∇ • σ s -ϕ∇p fm + f + ρ s ϕg, ρ f (1 -ϕ) ∂ t u + (u • ∇) u = ∇ • σ fm + ϕ∇p fm -f + ρ f (1 -ϕ) g, (2)
where g is the gravity force and 0 < ϕ < 1 is the solid volume fraction. The total stress tensors are

σ f = -p f I + τ f , σ s = -p s I + τ s , σ fm = -p fm I + τ fm ,
with p ∈ R the pressure for the solid and fluid phases in the mixture (subscript s and f m respectively), and for the fluid in the upper layer (subscript f ). I is the 2D identity tensor and the deviatoric stress tensors are given by

τ f = 2η f D(u f ), τ s = 2η s D(v), τ fm = 2η fm D(u),
where η f ∈ R is the constant viscosity of the fluid phase, and the viscosity of the solid phase η s ∈ R is defined by the rheological law for the granular phase. The strain-rate tensors are defined as usual by

D(u f ) = 1 2 (∇u f + (∇u f ) ), D(v) = 1 2 (∇v + (∇v) ), D(u) = 1 2 (∇u + (∇u) ).
In system (2), the terms ϕ∇p fm and f are the buoyancy and the drag forces between the phases. We consider here that

f = β (u -v) ,
where β is the drag coefficient. Following precedent works, this coefficient is defined as

β = (1 -ϕ) 2 η f κ , with κ = d 2 s (1 -ϕ) 3 150ϕ 2 , ( 3
)
where κ is the hydraulic permeability of the granular mass and d s the grain diameter [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF][START_REF] Iverson | A depth-averaged debris-flow model that includes the effects of evolving dilatancy. I. Physical basis[END_REF][START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF]. The Jackson's model is underdetermined and a closure equation must be added. In this work we assume the closure introduced in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] including the dilatancy effects

∇ • v = Φ, (4)
where the dilatancy function Φ is defined as in [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] and it is detailed in the subsection 2.2.

Boundary conditions

In this work we consider analogous boundary condition to [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] that we summarize here, we refer the reader to this reference for more details.

• At the bottom we consider the non penetration conditions

u • n = 0, v • n = 0 at the bottom,
where n is the space unit normal pointing out of the domain. This is completed with friction conditions for solid and fluid phases. At first, a Coulomb friction law,

(σ s n) τ = -tan µ eff v |v| (σ s n) • n at the bottom,
where µ eff is the effective friction coefficient given in terms of the rheology introduced in the next subsection and the subscript τ denotes the tangential projection. For the fluid phase we consider a Navier friction condition,

(σ fm n) τ = -k b u at the bottom,
for some coefficient k b ≥ 0.

• At the free surface we assume no tension for the fluid and the kinematic condition

σ f N X = 0; N t + u f • N X = 0 at the free surface,
where N = (N t , N X ) is a time-space normal to the free surface.

• At the interface, we consider the kinematic condition for the solid phase

Ñt + v • ÑX = 0 at the interface,
where Ñ = ( Ñt , ÑX ) is a time-space downward normal to the interface. Since we assume that there is an exchange of fluid mass between the layers, it must be conserved across the interface. Thus the fluid mass leaving the upper layer is the same than the fluid entering the mixture. If this quantity is denoted by V f we have

V f ≡ Ñt + u f • Ñx = (1 -ϕ * ) Ñt + (1 -ϕ * )u • Ñx at the interface,
where we write the Rankine-Hugoniot condition for equations

∇ • u f = 0 and ∂ t (1 -ϕ) + ∇ • ((1 -ϕ)u) = 0.
That can be also written in terms of the relative velocity as

- ϕ * 1 -ϕ * V f = (u f -u) • Ñx = -ϕ * ( Ñt + u • Ñx ).
Notice that the value ϕ * is introduced to denote the value of the solid volume fraction at the interface, which is discontinuous since there is no solid phase in the upper layer. The conservation of the total momentum gives,

ρ f V f (u -u f ) + (σ s + σ fm ) ÑX = σ f ÑX at the interface.
The energy balance through the interface yields the stress transfer condition

σ s ÑX = ρ f 2 (u -u f ) • ÑX | ÑX | 2 + (σ fm ÑX ) • ÑX | ÑX | 2 -p fm ϕ * 1 -ϕ * ÑX .
Finally a Navier fluid friction condition is considered at this level

σ fm + σ f 2 ÑX τ = -k i (u f -u) τ at the interface,
where k i ≥ 0 is a friction coefficient.

Rheology for fluidized granular flows

A rheology accounting for the behavior of immersed granular flows must be considered. For dry granular flows the µ(I)-rheology [START_REF] Jop | A constitutive law for dense granular flows[END_REF] is usually considered where the viscosity of the granular material is given by µ(I)ps D(v) , with µ(I) being the friction coefficient and I the inertial number. Dilatancy affects the friction law and can be introduced as follows (see [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF][START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF]):

η s = (µ(I) + tan ψ) p s D(v) ,
where D = √ 0.5 D : D and ψ is the dilatation angle as is defined by equation [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF] (see [START_REF] Roux | Texture-Dependent Rigid-Plastic Behavior[END_REF][START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF][START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF]). We introduce µ eff = µ(I) + tan ψ the effective friction coefficient, appearing also in the Coulomb friction law. The friction coefficient is

µ(I) = µ s + µ 2 -µ s I 0 + I I, (5)
which depends on the inertial number I, where I 0 , µ 2 > µ s are constant values depending on the material. The main difference between the dry and the fluidized case is the definition of the inertial number I (see [START_REF] Andreotti | Granular Media: Between Fluid and Solid[END_REF]). In the case of immersed granular flows in the viscous regime, this dimensionless number is defined as

I = 2 η f D(v) p s , (6)
in contrast to

I dry = 2d D(v) √ ps/ρs
that is the inertial number for the dry case with d being the particle diameter.

Following this idea, we can also consider an alternative simplified rheology where the friction coefficient is defined as

µ(I) = µ s + K 1 I, (7)
where K 1 is a constant value. In this case, the viscosity coefficient can be written as follows,

η s = (µ s + tan ψ) p s D(v) + 2K 1 η f .
Note that the viscosity coefficient is not defined in the case D(v) = 0. As in previous works for dry granular flows [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF][START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: a well-balanced multilayer discretization[END_REF] we consider a regularization method (e. g. [START_REF] Bercovier | A finite-element method for incompressible non-Newtonian flows[END_REF][START_REF] Lusso | Two-dimensional simulation by regularization of free surface viscoplastic flows with Drucker-Prager yield stress and application to granular collapse[END_REF], since the results are good enough (although zero velocity is not obtained in static solutions, but very small velocities) and they are cheaper computationally than using a duality method (e. g. [START_REF] Ionescu | Viscoplastic modeling of granular column collapse with pressure-dependent rheology[END_REF][START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF])). Then, the regularized viscosity coefficient is

η s = (µ(I) + tan ψ) p s D(v) 2 + δ 2 , ( 8
)
with δ > 0 the regularization parameter. Remark that other regularization techniques might be used (see e.g. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes with a µ(I)-rheology[END_REF][START_REF] Papanastasiou | Flows of Materials with Yield[END_REF]).

Finally, the dilatancy function is adopted from [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] and it is given by:

Φ = γ tan ψ (9)
with γ = 2 D(v) the shear rate, ψ the dilatation angle and

tan ψ = K (ϕ -ϕ eq c ) , ϕ eq c = ϕ stat c -K 2 I, (10)
which allows us to rewrite (4) as

∇ • v = K ϕ -ϕ stat c + K 2 I ,
where K is the dilatation constant and K 2 is a constant value, ϕ stat c is a constant volume fraction corresponding to a static equilibrium, and ϕ eq c the critical-state equilibrium compacity. This value determines if the current state of the granular medium is dilatation, contraction or equilibrium. Assuming a deformation occurs, γ > 0, the solid dilates if ϕ > ϕ eq c and contracts if ϕ < ϕ eq c . In the dilatation case, the fluid is sucked into the fluid-solid mixture and the pore pressure decreases. In the other case the fluid is expelled from the mixture and the pore pressure increases. If ϕ = ϕ eq c then there is no dilatation nor contraction.

Local coordinates

We write the system for the mixture in tilted coordinates. Let b(x) be an inclined fixed plane of constant angle θ with respect to the horizontal axis, we define the coordinates (x, z) ∈ Ω×R + ⊂ R 2 . The x (respectively z) axis is measured along the inclined plane (respectively the normal direction) as depicted in Figure 2. In this reference frame the gravity force is written as g = (-g sin θ, -g cos θ) .

We also set b(x) an arbitrary bottom topography and the layer of the fluidized material over it with thickness h(t, x), which are measured in the normal direction to the inclined plane b(x). The fluid and granular velocities are u = (u, u z ) and v = (v, v z ) where u, v denotes the downslope velocities of the fluid and granular phase, and u z , v z the normal components of u, v. Finally, we set ∇ = (∂ x , ∂ z ), the usual differential operator in the space variables.

We denote the components of the total stress tensor for the fluid phase as:

τ fm = τ xx fm τ xz fm τ xz fm τ zz fm and D(u) = 1 2 2∂ x u ∂ z u + ∂ x u z ∂ z u + ∂ x u z 2∂ z u z
, and analogously defined for the solid phase. With these definitions and embedding the mass equations into the momentum equations, system (2) is written as

                                               ρ s ∂ t ϕ + ∂ x (ϕv) + ∂ z (ϕv z ) = 0, ρ f ∂ t (1 -ϕ) + ∂ x ((1 -ϕ) u) + ∂ z ((1 -ϕ) u z ) = 0, ρ s ∂ t (ϕv) + ∂ x (ϕv 2 ) + ∂ z (ϕvv z ) + ∂ x p s = ∂ x τ xx s + ∂ z τ xz s -ϕ∂ x p fm + β (u -v) -ρ s ϕg sin θ, ρ f ∂ t ((1 -ϕ) u) + ∂ x ((1 -ϕ) u 2 ) + ∂ z ((1 -ϕ) uu z ) + ∂ x p fm = ∂ x τ xx fm + ∂ z τ xz fm +ϕ∂ x p fm -β (u -v) -ρ f (1 -ϕ) g sin θ, ρ s ∂ t (ϕv z ) + ∂ x (ϕvv z ) + ∂ z (ϕ(v z ) 2 ) + ∂ z p s = ∂ x τ xz s + ∂ z τ zz s -ϕ∂ z p fm + β (u z -v z ) -ρ s ϕg cos θ, ρ f ∂ t ((1 -ϕ) u z ) + ∂ x ((1 -ϕ) uu z ) + ∂ z ((1 -ϕ) (u z ) 2 ) + ∂ z p fm = ∂ x τ xz fm + ∂ z τ zz fm +ϕ∂ z p fm -β (u z -v z ) -ρ f (1 -ϕ) g cos θ, (11)
and the dilatancy closure equation is

∂ x v + ∂ z v z = Φ.
The boundary conditions are detailed directly in section 3 after the asymptotic analysis.

3 Two-phase multilayer models with dilatancy

In this section we deduce the multilayer model for fluidized granular flows focusing on the points that are different from previous multilayer models (e.g. [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF]). The model is obtained from a dimensional analysis and considering a multilayer approach of the domain in the vertical direction. We focus on the approximation of the deviatoric tensor at the internal interfaces, as well as the definition of the vertical velocity and the stress for both, the solid and the fluid phases.

A difficulty of this model is that the solid pressure is not purely hydrostatic: an excess pore fluid pressure appears in the equations as a consequence of the contraction/dilatation of the granular material (see ( 29)-( 30)). This term also appears in the vertical velocities of both phases [START_REF] Pelanti | A Roe-type scheme for two-phase shallow granular flows over variable topography[END_REF]. Finally, the model and the explicit expressions for the mass transference terms are presented.

Dimensional analysis

We perform here a dimensional analysis of the system [START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: a well-balanced multilayer discretization[END_REF]. We consider a shallow domain by assuming that the dimensionless parameter ε = H/L is small, where H and L are the characteristic height and length of the flowing mass, respectively. We denote as U the characteristic velocity and the characteristic time is then given by L/U . We define the dimensionless variables, denoted with the tilde symbol (.), as follows:

(x, z, t) = (L x, H z, (L/U ) t), h = H h, ρ s,f = ρ 0 ρ s,f , ϕ = ϕ 0 ϕ (u, u z ) = (U u, εU u z ), (v, v z ) = (U v, εU v z ), p s,f = ρ 0 U 2 p s,f , η s,f = ερ 0 U H η s,f , κ i = U κ i , k b = U k b τ xx s,f , τ xz s,f , τ zz s,f = ϕ 0 ρ 0 U 2 ε 2 τ xx s,f , ε τ xz s,f , ε 2 τ zz s,f , V f = εU V f .
Note that since

D(u) = U H 1 2   2ε 2 ∂ x u ε∂ z u + ε 3 ∂ x u z ε∂ z u + ε 3 ∂ x u z 2ε 2 ∂ z u z   ,
we obtain that andD(v) and the components of τ s are analogously written. In this work we have considered that the characteristic velocity for the solid phase is equal to the one of the fluid, and we also take the characteristic solid volume fraction is ϕ 0 = 1. We also define the Froude number

τ xx f = η∂ x u, τ xz f = η 2 ∂ z u + ε 2 ∂ x u z , τ zz f = η∂ z u z ,
F r = U √ gH cos θ ,
and the nondimensional drag coefficient

β = ρ 0 U L ε k β,
where k = -1 strong friction; 0 moderate friction, [START_REF] Fernández-Nieto | A Multilayer Method for the Hydrostatic Navier-Stokes Equations: A Particular Weak Solution[END_REF] and the dilatancy function

Φ = U H Φ.
Thus, the system (11) can be rewritten using dimensionless variables as (tildes have been dropped for simplicity):

                                                                         ρ s ∂ t ϕ + ∂ x (ϕv) + ∂ z (ϕv z ) = 0, ρ f ∂ t (1 -ϕ) + ∂ x ((1 -ϕ) u) + ∂ z ((1 -ϕ) u z ) = 0, ρ s ∂ t (ϕv) + ∂ x ϕv 2 + ∂ z (ϕvv z ) + ∂ x p s = ε 2 ∂ x τ xx s + ∂ z τ xz s -ϕ∂ x p fm + ε k β (u -v) - ρ s ϕ εF r 2 tan θ, ρ f ∂ t ((1 -ϕ) u) + ∂ x (1 -ϕ) u 2 + ∂ z ((1 -ϕ) uu z ) + ∂ x p fm = ε 2 ∂ x τ xx fm + ∂ z τ xz fm + ϕ∂ x p fm -ε k β (u -v) - ρ f (1 -ϕ) εF r 2 tan θ, ε 2 ρ s ∂ t (ϕv z ) + ∂ x (ϕvv z ) + ∂ z ϕ(v z ) 2 + ∂ z p s = ε 2 ∂ x τ xz s + ε 2 ∂ z τ zz s -ϕ∂ z p fm + ε k+2 β (u z -v z ) - ρ s ϕ F r 2 , ε 2 ρ f ∂ t ((1 -ϕ) u z ) + ∂ x ((1 -ϕ) uu z ) + ∂ z (1 -ϕ) (u z ) 2 + ∂ z p fm = ε 2 ∂ x τ xz fm + ε 2 ∂ z τ zz fm + ϕ∂ z p fm -ε k+2 β (u z -v z ) - ρ f (1 -ϕ) F r 2 (13a) (13b) (13c) (13d) (13e) (13f)
and

∂ x v + ∂ z v z = 1 ε Φ. ( 13g
)
Before applying the multilayer approach it is suitable to write the system in matrix notation. To this aim we collapse all terms coming from the stress tensors under the following notation for the fluid phase (analogously for the solid phase) σ ε,fm = -p fm E + ετ ε,fm where the subindex ε marks the dependence of these terms on ε, as in [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF], with

τ ε,fm = ηD ε,fm (u), D ε,fm (u) := 1 2 2ε 2 ∂ x u ∂ z u + ε 2 ∂ x u z ∂ z u + ε 2 ∂ x u z 2 ∂ z u z
, and

E = ε 0 0 1/ε . Then, defining F = tan θ F r 2 ,
1 εF r 2 , the system (13) is finally written as

                                         ρ s ∂ t ϕ + ∇ • (ϕv) = 0, ρ f ∂ t (1 -ϕ) + ∇ • ((1 -ϕ) u) = 0, ρ s ∂ t (ϕv) + ∇ • (ϕv ⊗ v) - 1 ε ∇ • σ ε,s + 1 ε ϕ∇ • (p fm E) = ε k β (u -v) - 1 ε ρ s ϕF , ρ f ∂ t ((1 -ϕ) u) + ∇ • ((1 -ϕ) u ⊗ u) - 1 ε ∇ • σ ε,fm - 1 ε ϕ∇ • (p fm E) = -ε k β (u -v) - 1 ε ρ f (1 -ϕ) F , ∇ • v = 1 ε Φ. ( 14
)
Next we specify the asymptotic expressions of the boundary conditions described in section 2.1. Taking into account the asymptotic analysis up to second order performed in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] (see its Appendixes A and B for the details), these conditions read

(1 -ϕ * ) ∂ t h + u | z=b+h ∂ x (b + h) -u z | z=b+h = -V f ; ∂ t h + v | z=b+h ∂ x (b + h) -v z | z=b+h = 0; p fm| z=b+h = ρ f g cos θh f , p s| z=b+h = 0, (15)
The energy balance gives

η f 2 ∂ z u | z=b+h = 1 ε ε 2 ρ f V f -κ i (u f -u) | z=b+h
where u f is the first component of the velocity of the water in the upper layer. At the bottom z = b the no penetration and the friction conditions read

u | z=b ∂ x b = u z | z=b ; (η f ∂ z u) | z=b = 1 ε k b u | z=b , v | z=b ∂ x b = v z | z=b ; η s 2 ∂ z v | z=b = 1 ε (µ(I) + tan ψ) p s v |v| | z=b .
Moreover we consider that there is no transference of mass at the bottom, neither solid nor fluid.

A multilayer approach

We remind the notation of the multilayer domain in order to apply this approach to system [START_REF] George | A depth-averaged debris-flow model that includes the effects of evolving dilatancy. II. Numerical predictions and experimental tests[END_REF]. The mixture domain is denoted by Ω F (t) and

I F (t) = x ∈ R; (x, z) ∈ Ω F (t) ,
is its projection on the reference inclined plane, for any time t > 0. The idea of this approach is to divide the domain along the direction normal to the plane in N ∈ N \ {0} shallow layers with thickness h α (t, x), such that the total height is h = N α=1 h α (see Figure 2). In practice, the normal partition is previously fixed by using the positive coefficient l α satisfying

h α = l α h for α = 1, ..., N ; N α=1 l α = 1.
The N + 1 interfaces separating the vertical layers are denoted by Γ α+ 1 2 (t), and are described by z = z α+ 1 2 (t, x) for α = 0, 1, .., N , x ∈ I F (t). These interfaces, which are assumed to be smooth enough, shut the subdomain Ω α (t) up, i.e., Ω α (t) = (x, z); x ∈ I F (t) and z α-

1 2 < z < z α+ 1 2 .
A key point of this approach is the approximation of the variables at the interfaces. Then, for an arbitrary function f and for α = 0, 1, ..., N , we set

f - α+ 1 2 := (f | Ωα(t) ) | Γ α+ 1 2 (t) and f + α+ 1 2 := (f | Ω α+1 (t) ) | Γ α+ 1 2 (t) ,
where if the function f is continuous, then

f α+ 1 2 := f | Γ α+ 1 2 (t) = f + α+ 1 2 = f - α+ 1 2 .
Finally, the normal vectors at the interfaces are defined. For a time t ≥ 0,

n T,α+ 1 2 = ∂ t z α+ 1 2 , ∂ x z α+ 1 2 , - 1
1 + ∂ x z α+ 1 2 2 + ∂ t z α+ 1 2 2 and n α+ 1 2 = ∂ x z α+ 1 2 , - 1
1 + ∂ x z α+ 1 2 2
denote the space-time unit normal vector and the space unit normal vector to the interface Γ α+ 1 2 (t) outward to the layer Ω α+1 (t) for α = 0, ..., N (see detail in Figure 2).

Weak solution with discontinuities

In order to apply the multilayer approach, we look for a particular weak solution (u, v, p fm , p s , ρ f , ρ s , ϕ) of ( 14) (see [START_REF] Fernández-Nieto | A Multilayer Method for the Hydrostatic Navier-Stokes Equations: A Particular Weak Solution[END_REF][START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF]) satisfying that:

(i) (u, v, p fm , p s , ρ f , ρ s , ϕ) is a standard weak solution of ( 14) in each layer Ω α (t), (ii) (u, v, p fm , p s , ρ f , ρ s , ϕ) satisfies the normal flux jump condition for mass and momentum equations for each phase, at the interfaces Γ α+ 1 2 (t), namely:

[(ρ s ϕ ; ρ s ϕv)] α+ 1 2 n T,α+ 1 2 = 0, [(ρ f (1 -ϕ) ; ρ f (1 -ϕ) u)] α+ 1 2 n T,α+ 1 2 = 0, (16)
[(ρ s ϕv; ρ s ϕv ⊗ v - 1 ε σ ε,s ) α+ 1 2 n T,α+ 1 2 = 0, (ρ f (1 -ϕ) u; ρ f (1 -ϕ) u ⊗ u - 1 ε σ ε,fm ) α+ 1 2 n T,α+ 1 2 = 0, (17)
where [(a; b)] α+ 1 2 denotes the jump of (a; b) across the interface Γ α+ 1 2 (t).

In the previous normal flux jump condition for the momentum equation, we have used that the fluid pressure p fm is continuous and the gravitational force and friction between the phases are source terms, therefore there is no jump in the bouyancy term and all the terms at the interfaces. Thus, there is no significant differences between the jump condition in [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF] and our case. A particular family of test function is considered where the horizontal component of the velocity does not depend on z, and the vertical one is linear on z. Then, the velocities of the granular and the fluid phase in each layer are

u | Ωα(t) := u α := (u α , u z α ) , v | Ωα(t) := v α := (v α , v z α )
, where u α ,v α and u z α , v z α are the horizontal and normal velocities, respectively, on layer α. Then,

∂ z u α = 0, ∂ z v α = 0; and ∂ z u z α = d f,α (t, x), ∂ z v z α = d s,α (t,
x) for some smooth functions d f,α (t, x), d s,α (t, x). We also consider that the solid volume fraction is constant in each layer

ϕ | Ωα(t) := ϕ α = ϕ α (t, x).
From the jump condition for the mass equations ( 16) we can easily obtain the definition of the solid mass flux (G s,α+ 1 2 ) and the fluid mass flux (G f,α+ 1 2 ):

G s,α+ 1 2 := G + s,α+ 1 2 = G - s,α+ 1 2 and G f,α+ 1 2 := G + f,α+ 1 2 = G - f,α+ 1 2 with G ± s,α+ 1 2 = ϕ ± α+ 1 2 ∂ t z α+ 1 2 + v ± α+ 1 2 ∂ x z α+ 1 2 -(v z α+ 1 2 ) ± , (18a)
G ± f,α+ 1 2 = (1 -ϕ) ± α+ 1 2 ∂ t z α+ 1 2 + u ± α+ 1 2 ∂ x z α+ 1 2 -(u z α+ 1 2 ) ± . ( 18b
)
Note that condition [START_REF] Ionescu | Viscoplastic modeling of granular column collapse with pressure-dependent rheology[END_REF] gives the mass transference at the top of the mixture layer z = z N +1/2 (t, x), i.e. the exchange with the upper fluid domain,

G f,N +1/2 = -V f and G s,N +1/2 = 0. ( 19
)
Now, after some straightforward calculations we get the approximation of the deviatoric tensor at the interface Γ α+ 1 2 using (17) as in [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF]:

τ ± ε,s,α+ 1 2 n α+ 1 2 = τ ε,s,α+ 1 2 n α+ 1 2 ± 1 2 ρ s G s,α+ 1 2 1 + ∂ x z α+ 1 2 2 [v] | Γ s,α+ 1 2
,

where τ ε,s,α+ 1 2 is an approximation of η s D ε (v α ) |Γ α+ 1 2
, defined by

τ ε,s,α+ 1 2 = η s,α+ 1 2 D ε,s,α+ 1 2 = 1 2 η s,α+ 1 2       2ε 2 ∂ x v + α+ 1 2 + v - α+ 1 2 2 D xz ε,s,α+ 1 2 D xz ε,s,α+ 1 2 2 V V Z ,α+ 1 2       , with D xz ε,s,α+ 1 2 = ε 2 ∂ x (v z α+ 1
2

) + + (v z a+ 1 2 ) - 2 + V H Z ,α+ 1 2 .
As in previous works,

V Z = (V H Z , V V Z ) is an approximation of the derivatives in z of the velocity V Z -∂ z v = 0, at the interfaces.
Analogously, we can define for the fluid phase τ ε,fm,α+

1 2 , D ε,fm,α+ 1 2 and U Z ,α+ 1 2
and get the approximation of the deviatoric tensor of the fluid phase:

τ ± ε,fm,α+ 1 2 n α+ 1 2 = τ ε,fm,α+ 1 2 n α+ 1 2 ± 1 2 ρ f G f,α+ 1 2 1 + ∂ x z α+ 1 2 2 [u] | Γ f,α+ 1 2
.

Finally, we need to approximate η s,α+ 1 2 , the viscosity coefficient [START_REF] Capart | Depth-integrated equations for entraining granular flows in narrow channels[END_REF], up to first order in ε. We consider

D(v) α+ 1 2 ≈ 1 2 V H Z ,α+ 1 2 ,
then, the viscosity coefficient is

η s,α+ 1 2 = µ(I α+ 1 2 ) + (tan ψ) α+ 1 2 p s,α+ 1 2 V H Z ,α+ 1 2 2 /4 + δ 2 , ( 20
) with V H Z ,α+ 1 2 = v α+1 -v α h α+ 1 2 , for α = 1, . . . , N -1,
and h α+ 1 2 the distance between the midpoints of layers α and α + 1. For the particular case α = 0 we take

V H Z , 1 2 = λv 1 h 1 , (21)
with λ depending on the friction condition at the bottom, λ = 1 (friction), λ = 2 (no slip). p s,α+ 1 2 is the solid pressure at the interface z α+ 1

2

, whose expression will be given later, and

I α+ 1 2 = η f V H Z ,α+ 1 2 p s,α+ 1 2 , for α = 0, . . . , N -1. (22)
Finally, thanks to [START_REF] Ionescu | Viscoplastic modeling of granular column collapse with pressure-dependent rheology[END_REF] we trivially have η s,N +1/2 = 0.

Vertical velocities

The procedure to recover the vertical velocities is the same as in [START_REF] Fernández-Nieto | A Multilayer Method for the Hydrostatic Navier-Stokes Equations: A Particular Weak Solution[END_REF]. In order to compute the vertical velocity for the solid phase, we use both the solid mass transference term and the dilatancy closure equation. Firstly, since there is not transference of solid mass at the bottom level, we have G s,1/2 = 0, i.e., we get (v z 1 2

) + = v 1 ∂ x b + ∂ t b. Now, integrating the dilatancy equation (13g) between z α-1 2 and z ∈ z α-1 2 , z α+ 1 2
we obtain that

v z α (z) = (v z α-1 2 ) + -(z -z α-1 2 ) ∂ x v α - 1 ε Φ α ,
and from the solid mass transference term (18a) we get

(v z α+ 1
2

) + = 1 ϕ α+1 (ϕ α+1 -ϕ α ) ∂ t z α+ 1 2 + (ϕ α+1 u α+1 -ϕ α u α ) ∂ x z α+ 1 2 + ϕ α (v z α+ 1 2 ) -, with (v z α+ 1
2

) -= (v z α-1 2 ) + -h α ∂ x v α - 1 ε Φ α .
Computing the vertical velocity of the fluid phase is possible by using the incompressibility of the mixture

∇ • (ϕv + (1 -ϕ) u) = 0. ( 24
)
Analogously to the procedure above, we can obtain the vertical velocity for the fluid phase using ( 24), (18b) and (25b), resulting:

(u z 1 2
)

+ = u 1 ∂ x b + ∂ t b, (25a)
u z α (z) = (u z α-1 2 ) + - z -z α-1 2 (1 -ϕ α ) ∂ x (ϕ α v α + (1 -ϕ α ) u α ) -ϕ α ∂ x v α + 1 ε ϕ α Φ α , (25b)
(u z α+ 1
2

) + = 1 (1 -ϕ α+1 ) (ϕ α -ϕ α+1 ) ∂ t z α+ 1 2 + (1 -ϕ α+1 ) u α+1 -(1 -ϕ α ) u α ∂ x z α+ 1 2 + (1 -ϕ α ) (u z α+ 1 2 ) -, (25c)
and

(u z α+ 1
2

) -= (u z α-1 2 ) + - h α (1 -ϕ α ) ∂ x (ϕ α v α + (1 -ϕ α ) u α ) -ϕ α ∂ x v α + 1 ε ϕ α Φ α . (25d)
Explicit expression for fluid and solid pressures: p fm , p s

From the nondimensional vertical momentum equation (13f) we obtain up to first order that

∂ z p fm = ϕ∂ z p fm -ε k+2 β(u z -v z ) - ρ f (1 -ϕ) F r 2 . ( 26
)
By integrating this expression from z to b + h (thanks to the continuity of the dynamic fluid pressure), we obtain

p fm,α (z) = ρ f F r 2 (b + h + h f -z) + p e f,α , (27a)
where

p e f,α = p e f,α+ 1 2 + ε k+2 z α+ 1 2 z β (1 -ϕ) (u z -v z ) dz , (27b)
with

p e f,α+ 1 2 = p e f,α+1 z = z α+ 1 2 = ε k+2 N γ=α+1 z γ+ 1 2 z γ-1 2 β (1 -ϕ) (u z -v z ) dz . ( 27c
)
Once we have the fluid pressure, we can obtain an explicit expression for the solid one. From the vertical momentum equation (13e) for the solid phase we have 26) in the previous equation and integrating from z to z α+ 1 2 we get

∂ z p s = -ϕ∂ z p fm + ε k+2 β(u z -v z ) - ρ s ϕ F r 2 . Using (
p s,α (z) = p s,α+ 1 2 + ϕ α (ρ s -ρ f ) F r 2 z α+ 1 2 -z -ε k+2 z α+ 1 2 z β (1 -ϕ) (u z -v z ) dz ,
where

p s,α+ 1 2 = (ρ s -ρ f ) F r 2 N γ=α+1 h γ ϕ γ -p e f,α+ 1 2 . ( 28
)
We can write

p s,α (z) = (ρ s -ρ f ) F r 2 N γ=α+1 ϕ γ h γ + ϕ α (ρ s -ρ f ) F r 2 z α+ 1 2 -z -p e f,α , (29)
where p e f,α is given by (27b). The computation of the excess pore pressure is subtle. It is carried out using the dilatancy equation, the mass transference terms for the solid and the fluid, and the fact that we obtain the incompressibility of the mixture with velocity (ϕv + (1 -ϕ) u) from the mass equations in [START_REF] George | A depth-averaged debris-flow model that includes the effects of evolving dilatancy. II. Numerical predictions and experimental tests[END_REF]. We also assume that we are in the case where the friction between the solid and the fluid phase is strong (k = -1 in ( 12)). This leads to (details in appendix A)

p e f,α (z) = p e f,α+ 1 2 -β α h 2 α -z -z α-1 2 2 2 (1 -ϕ α ) 2 Φ α - β α ϕ α (1 -ϕ α ) 2 z α+ 1 2 -z α-1 γ=1 ϕ γ h γ Φ γ , (30a)
with

p e f,α+ 1 2 = p e f,α+1 z = z α+ 1 2 = N ξ=α+1 -β ξ h ξ ϕ ξ (1 -ϕ ξ ) 2   ξ-1 γ=1 ϕ γ h γ Φ γ + ϕ ξ h ξ Φ ξ 2   . ( 30b
)

Final model

We derive the final model for the mixture by looking for a particular weak solution of system [START_REF] George | A depth-averaged debris-flow model that includes the effects of evolving dilatancy. II. Numerical predictions and experimental tests[END_REF] with the procedure introduced in [START_REF] Fernández-Nieto | A Multilayer Method for the Hydrostatic Navier-Stokes Equations: A Particular Weak Solution[END_REF], which was applied to the case of a variable viscosity in [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF].

As we mentioned before, for the upper layer we adopt the model derived in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF]. Finally, the resulting system must be written in dimensional variables. As conclusion, the final model reads, for α = 1, ..., N ,

                                                                                                                     ∂ t h f + ∂ x (h f u f ) = -G f,N + 1 2 , ρ f ∂ t (h f u f ) + ∂ x h f u 2 f + g cos θh f ∂ x (b + b + h + h f ) = - 1 2 ρ f G f,N + 1 2 (u f + u) -k i (u f -u), l α ϕ α ∂ t h + ∂ x (hv α ) = G s,α+ 1 2 -G s,α-1 2 + ϕ α l α hΦ α , l α ∂ t (hϕ α ) + ∂ x (hϕ α v α ) = G s,α+ 1 2 -G s,α-1 2 , l α ∂ t (h (1 -ϕ α )) + ∂ x (h (1 -ϕ α ) u α ) = G f,α+ 1 2 -G f,α-1 2 , l α ρ s ∂ t (hϕ α v α ) + ρ s ∂ x hϕ α v 2 α + g cos θhϕ α ∂ x ρ s b + b + h + ρ f h f + g cos θh (ρ s -ρ f ) N γ=α+1 h γ ∂ x ϕ γ + h α 2 ∂ x ϕ α + g cos θh (ρ s -ρ f ) N γ=α+1 (ϕ γ -ϕ α ) ∂ x h γ -h (1 -ϕ α ) ∂ x p e f,α = l α hβ α (u α -v α ) + K s,α-1 2 -K s,α+ 1 2 + 1 2 ρ s G s,α+ 1 2 (v α+1 + v α ) - 1 2 ρ s G s,α-1 2 (v α + v α-1 ) , l α ρ f ∂ t (h (1 -ϕ α ) u α ) + ρ f ∂ x h (1 -ϕ α ) u 2 α + ρ f g cos θh (1 -ϕ α ) ∂ x b + b + h + h f + h (1 -ϕ α ) ∂ x p e f,α = -l α hβ α (u α -v α ) + K f,α-1 2 -K f,α+ 1 2 + 1 2 ρ f G f,α+ 1 2 (u α+1 + u α ) - 1 2 ρ f G f,α-1 2 (u α + u α-1 ) , (31a) (31b) (31c) (31d) (31e) (31f) (31g)
where

β α = 1 h α z α-1 2 z α-1 2 β α dz and ∂ x p e f,α = 1 h α z α-1 2 z α-1 2 ∂ x p e f,α dz, (31h)
where β α and p e f,α are given by ( 3) and (30), respectively. The mass transference terms G s,α+ 1 2 , G f,α+ 1 2 are given by ( 18), and

K s,α+ 1 2 = - 1 2 η s,α+ 1 2 (V H Z ,α+ 1 2 ) V H Z ,α+ 1 2 , K s, 1 2 = -ρ s g cos θ h µ I 1 2 + (tan ψ) 1/2 v 1 |v 1 | , K f,α+ 1 2 = -η f U H Z ,α+ 1 2 , K f, 1 2 = -k b u 1 , (31i) for V H Z ,α+ 1 2 , η s,α+ 1 2
defined in ( 20)- [START_REF] Lusso | Two-dimensional simulation by regularization of free surface viscoplastic flows with Drucker-Prager yield stress and application to granular collapse[END_REF]. For the specific case of the bottom and mixture interface, we introduce the following definitions:

u 0 = 0; v 0 = (1 -λ)v 1 ; u N +1 = u f ; v N +1 = 0;
notice that the definition of v 0 comes from [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes with a µ(I)-rheology[END_REF]. Moreover, the frictions between the mixture and the upper fluid layer are defined by

K s,N + 1 2 = 0, and K f,N + 1 2 = k i (u f -u).
At the mixture interface we also remind that the transference terms are given in [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF]. Note that although the previous system has 5N equations, only 5N -1 are linearly independent, therefore we have a system with 5N -1 equations and unknowns, which are:

(h, {ϕ α , hu α , hv α } α=1,...,N , {G f,α+ 1 2 , G s,α+ 1 2 } α=1,...,N -1 ).
However, an explicit expression can be found for the mass transference terms, as detailed below.

Explicit mass transference terms

In order to obtain explicit expressions for the solid mass transference term, we start by combining the averaged mass equations for the solid and fluid phase (31d)-(31e), giving

l α ∂ t h + ∂ x h ϕ α v α + (1 -ϕ α ) u α = G s,α+ 1 2 + G f,α+ 1 2 -G s,α-1 2 + G f,α-1 2 . ( 32
)
Next, denoting

U α = ϕ α v α + (1 -ϕ α ) u α ,
and summing up the previous equation from α = 1, . . . , N we obtain

∂ t h + ∂ x N α=1 l α h U α = G s,N + 1 2 + G f,N + 1 2 . ( 33
)
Here we have used that there is no mass transference at the bottom, i.e. G s,

1 2 = G f, 1 2 = 0. The transference terms G s,N + 1 2 and G f,N + 1 2
are defined in [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF]. Using previous equation in the averaged dilatancy equation (31c), it leads to

l α ϕ α -∂ x N γ=1 l γ h U γ + ∂ x (hv α ) -hΦ α + l α ϕ α G s,N + 1 2 + G f,N + 1 2 = G s,α+ 1 2 -G s,α-1 2 ,
and summing up the previous equation from α = 1, . . . , N , it gives

G s,α+ 1 2 = α β=1 l β ϕ β ∂ x hv β - N γ=1 l γ h U γ -hΦ β + G s,N + 1 2 + G f,N + 1 2 α β=1 l β ϕ β . (34a)
Once we have an explicit expression for the solid mass transference term, we can use it in order to obtain the fluid mass transference term. We embbed (33) in (32), and sum it up from α = 1, . . . , N , it gives

G f,α+ 1 2 = α β=1 l β ∂ x h U β - N γ=1 l γ h U γ -G s,α+ 1 2 + G s,N + 1 2 + G f,N + 1 2 α β=1 l β . ( 34b
)
Taking into account those expressions for the mass transference terms, system (31) can be rewritten as a system with 3N +1 equations and unknowns. The unknowns are now the total height (h), the solid volume fraction ({ϕ α } N α=1 ) and the solid and fluid velocities ({v α } N α=1 , {u α } N α=1 ). Thus, the equations of the final system are (31d), (31f), (31g) and (33).

Note that we need a closure relation for the fluid and solid mass transference at the free surface introduced in [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF]. Moreover, some important properties of the model, as the mass conservation, depend on how these terms are defined. This is discussed in subsection 3.4.1.

Since the proposed model is a two-phase multilayer model, its mathematical nature would be non-hyperbolic [START_REF] Abgrall | Two-Layer Shallow Water System: A Relaxation Approach[END_REF][START_REF] Chiapolino | Models and methods for two-layer shallow water flows[END_REF]. Nevertheless, the friction between phases is considered strong which leads to obtain similar velocities. Thus we expect that the system is not far from being hyperbolic, which is reinforced by the fact that we do not find numerical difficulties in the approximation of the model. Even in simple configurations, the numerical approximation of our model involves important difficulties, namely when approximating the pressure. For this reason, in the following we study the particular case of immersed uniform flows, for which numerical tests will be presented in section 4.

Particular configuration: immersed uniform flows

We consider the analogous case to [START_REF] Bouchut | A two-phase solid-fluid model for dense granular flows including dilatancy effects: comparison with submarine granular collapse experiments[END_REF] for uniform flows, i.e. an immersed flow with rigid top boundary at the surface (see Figure 3), in order to compare with their results. To this aim, we rewrite the derivative of the fluid depth taking into account that we have a rigid lid boundary (b = 0). This implies that

∂ x (x tan θ + h + h f ) = 0.
Since we consider the case of a uniform flow we also have that ∂ x h = 0. Then, from the previous equation we deduce that

∂ x h f = -tan θ.
Taking these hypotheses into account in (31c)-(31i) with u f = 0, k i = 0, we obtain the following Figure 3: Sketch of the domain in the immersed configuration and flat free surface.

simplified system

                                       ∂ t h = G s,N + 1 2 + G f,N + 1 2 , ∂ t ϕ α = -ϕ α Φ α , l α ρ s ∂ t (hϕ α v α ) = l α h β α (u α -v α ) -(ρ s -ρ f ) g sin θl α hϕ α + K s,α-1 2 -K s,α+ 1 2 + 1 2 ρ s G s,α+ 1 2 (v α+1 + v α ) - 1 2 ρ s G s,α-1 2 (v α + v α-1 ) , l α ρ f ∂ t (h (1 -ϕ α ) u α ) = -l α h β α (u α -v α ) + K f,α-1 2 -K f,α+ 1 2 + 1 2 ρ f G f,α+ 1 2 (u α+1 + u α ) - 1 2 ρ f G f,α-1 2 (u α + u α-1 ) , (35a) (35b) (35c) (35d)
Note that (35b) is deduced from (31c) and (31d). We also consider the following definitions: for the drag coefficient (see [START_REF] Audusse | A multilayer Saint-Venant model: derivation and numerical validation[END_REF])

β α = 150ϕ 2 α d 2 s (1 -ϕ α ) η f ,
and from ( 9)-( 10) the dilatancy function is

Φ α = γα K ϕ α -ϕ eq c,α , with ϕ eq c,α = ϕ stat c -K 2 I α , with γα = |v α -v α-1 | h α .
Note that in the definition of γα we set a backward approximation because this is the one that gives us better results, namely at the bottom where the friction boundary condition must be taken into account. In particular, the central approximation produced spurious oscillations.

The system (35) is written in terms of the mass transference at the mixture interface, appearing in the first equation (35a). An important property of this system that is shared with the system proposed in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] is that it preserves the solid mass in the mixture domain under the definitions given in [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF], that is,

G s,N + 1 2 = 0 and G f,N + 1 2
. In this case we get

∂ t h = G f,N + 1 2 ,
hence, thanks to (33) the total solid mass is preserved:

∂ t N β=1 h β ϕ β = 0.
In order to give an expression for the mass transference terms in the multilayer domain, we first use that G s,N + 1 2 = 0 in (34a) to obtain the fluid mass transference at the free surface, so finally we have

G s,N + 1 2 = 0, G f,N + 1 2 = N γ=1 l γ ϕ γ hΦ γ N γ=1 l γ ϕ γ . ( 36
)
With this definition, the solid and fluid mass transference terms are written as

G s,α+ 1 2 = α β=1 l β ϕ β G f,N + 1 2 -hΦ β and G f,α+ 1 2 = G f,N + 1 2 α β=1 l β -G s,α+ 1 2 , (37)
for α = 1, . . . , N -1.

From now on we refer to model ( 35), ( 36), (37) as PGM model (Preserving Granular Mass). Actually, this model is the multilayer extension of the model proposed by Bouchut et al. [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF], where the thickness of the mixture (h m in their notation) is not preserved. In fact, their case is the single-layer case of the proposed model. In that case, as authors assume there is no solid mass transference between these layers, therefore G s,3/2 = 0 and from (34a) the fluid mass transference at the free surface remains

G f,3/2 = h 1 Φ 1 + ∂ x h 1 (1 -ϕ 1 ) (u 1 -v 1 ) ,
which matches with -V f as defined in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF], where the minus sign comes from the fact that they consider the upward normal vectors while we work here with downward normal vectors. This section is organized as follows: in subsection 3.4.1 we propose a multilayer extension preserving the mixture height (as the Pailha and Pouliquen model [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF]). An analytical solution is deduced in subsection 3.4.2, by neglecting the friction between the phases. In subsection 3.4.3 we deduce an analytical solution that includes the effect of the side walls friction for confined flows, which plays an important role in the case of granular flows, namely in laboratory experiments. The numerical approximation of the resulting system (35) is a difficult task even in this simple configuration of a uniform flow. The main difficulty is the approximation of the excess pore fluid pressure term p e f,α , which is necessary to define I α . Its numerical approximation is a key point of the numerical scheme, which is detailed in subsection 3.4.4.

A multilayer model preserving the mixture height

As discussed in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] one of the main difference between their model and the one of Pailha and Pouliquen is that this last one preserves the total height of the mixture and not the total mass. This implies different steady state solutions for each model. In order to perform a similar comparison in the numerical tests we show a multilayer model that preserves the height instead of the granular mass in the mixture.

In order to preserve the total height of the mixture layer it is required that the solid mass leaving the mixture domain equals the fluid mass entering the domain, i.e., G s,N + 1 2 = -G f,N + 1 2 in the model (35). Note that in this case, the model in the upper layer should also take into account the granular phase.

From (35a) we obtain that ∂ t h = 0, and (31d) and (34a), under the uniform flow assumption, lead to

∂ t N β=1 h β ϕ β = G s,N + 1 2 , with G s,N + 1 2 = - N β=1 l β hϕ β Φ β . (38)
Therefore, the solid mass of the considered mixture may be not preserved at the domain. The solid and fluid mass transference terms (34) are

G s,α+ 1 2 = - α β=1 l β h ϕ β Φ β and G f,α+ 1 2 = -G s,α+ 1 2 . ( 39
)
for α = 1, . . . , N . From now on we refer to model ( 35), ( 38), (39) as the PH model (Preserving Height).

Analytical solution

The analytical solution for a uniform immersed flow can be obtained with an analogous procedure to the case of a granular flow, i.e, Bagnold flow (see e.g. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes with a µ(I)-rheology[END_REF]). To this aim, we consider the particular case where no friction effects appear in the model (even β(u -v)). This allows us to obtain an explicit expression for analytical solid velocity, pressure and solid volume fraction. The effects of the term β(u -v) will be studied in the numerical tests (see section 4.1.1).

Firstly, since we focus on steady uniform flows, where

G s,N + 1 2 = G f,N + 1 2
= 0, once the final height is computed (which will be different in each model), the analytical solution will be the same for both models analyzed in previous section. Starting from system [START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: a well-balanced multilayer discretization[END_REF] in the immersed configuration (∂ x h f = -tan θ), and assuming a steady uniform flow (tan ψ = 0), we obtain

∂ z τ xz s = (ρ s -ρ f ) ϕ(z) g sin θ, ∂ z p s = -(ρ s -ρ f ) ϕ(z) g cos θ.
(

) 40
Moreover we consider

τ s (z) = µ(I(z)) p s (z), µ(I(z)) = µ s + K 1 I(z), and ϕ(z) = ϕ eq c (z) = ϕ stat c -K 2 I(z).
Integrating from z > b to the free surface b + h we obtain

p s (z) = (ρ s -ρ f ) g cos θ b+h z ϕ(s)ds, (41)
τ xz s (z) = -(ρ s -ρ f ) g sin θ b+h z ϕ(s)ds.
Using that

τ xz s (z) = µ(I(z)) p s (z), we get -(ρ s -ρ f ) g sin θ b+h z ϕ(s)ds = µ(I(z)) (ρ s -ρ f ) g cos θ b+h z ϕ(s)ds, so (µ(I(z)) + tan θ) b+h z ϕ(s)ds = 0.
This leads to a constant friction coefficient µ(I) = -tan θ. Thus, we obtain that the inertial number

I(µ s , ∆µ = µ 2 -µ s , tan θ, K 1 ) is constant I = -tan θ -µ s ∆µ + µ s + tan θ I 0 , or I = -tan θ -µ s K 1 ,
depending of the considered rheology ( 5) or [START_REF] Bouchut | A two-phase solid-fluid model for dense granular flows including dilatancy effects: comparison with submarine granular collapse experiments[END_REF], and therefore the solid volume fraction

ϕ = ϕ eq c = ϕ stat c -K 2 I,
is also constant. In addition, using the definition of the inertial number (6) we get an expression for ∂ z v. Integrating from z = b to 0 < z < b + h and using the no-slip condition at the bottom, the velocity profile reads:

v(z) = (ρ s -ρ f ) I ϕ cos θ 2η f (b + h) 2 -(b + h -z) 2 .

Side walls friction for confined flows: analytical solution

To model side wall friction the following additional viscous term is considered at the momentum solid phase (see [START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: a well-balanced multilayer discretization[END_REF] and [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF]),

- 2 W µ w p s (z) v |v| , ( 42
)
where W is the width of the domain and µ w the friction coefficient at the lateral walls.

In order to obtain an analytical solution in this case we start again from (40), where the side walls friction term is added

∂ z τ xz s (z) - 2 W µ w p s (z) = (ρ s -ρ f ) ϕ(z) g sin θ, ∂ z p s (z) = -(ρ s -ρ s ) ϕ(z) g cos θ.
where we have also assumed that sign(∂ z v) = sign(v) = 1. Following the same procedure as in the previous section we get

τ xz s (z) + 2 W µ w (ρ s -ρ f ) g cos θ b+h z b+h s ϕ(ξ)dξ ds = -(ρ s -ρ f ) g sin θ b+h z ϕ(s)ds. ( 43
)
and writing (43) in terms of ϕ(z), we obtain

µ s + K 1 ϕ stat c -ϕ(z) K 2 b+h z ϕ(s)ds + 2 W µ w b+h z b+h s ϕ(ξ)dξ ds = -tan θ b+h z ϕ(s)ds. (44)
Now, deriving the previous equation, it leads to the integro-differential equation

ϕ (z) K 1 K 2 b+h z ϕ(s)ds + µ s + K 1 ϕ stat c -ϕ(z) K 2 ϕ(z) + 2 W µ w b+h z ϕ(s)ds = -tan θϕ(z).
Considering the variable

F (z) = b+h z ϕ(s)ds,
we obtain an initial value problem with two differential equations

F (z) ϕ(z) =   -ϕ(z) - K 2 K 1 F (z) µ s + K 1 ϕ stat c -ϕ(z) K 2 ϕ(z) + 2 W µ w F (z) + tan θϕ(z)   ,
where

F (β + h) = , 1, ϕ(b + h) = ϕ stat c + K 1 K 2 (tan θ + µ s ) .
This system is solved using a numerical ODE solver, obtaining the concentration ϕ(z) and F (z).

The pressure is computed using (41) and the inertial number is

I(z) = (ϕ stat c -ϕ(z)) /K 2 .
Using the definition of I we find the analytical solid velocity solving the differential equation

v (z) = p s (z)I(z) η f , with v(0) = 0.
Remark 1. Note that if we take (µ w = 0) in (44) we recover the previous case without side walls friction (subsection 3.4.2). In particular, it leads to µ(I) = -tan θ, which implies that both, the inertial number I and the solid volume fraction ϕ, are constant. Remark 2. In the case of uniform flows for dry granular flows and a hydrostatic pressure it is equivalent to consider the additional viscous term (42) or the following modification of the friction coefficient (see [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF])

µ(I) = µ(I) + µ w b + h -z W .
We show in what follows that for the case considered in this paper, where we deal with a two-phase flow, this modification of the friction coefficient, under the same hypothesis, is not equivalent but it is a second order approximation of it.

It is enough to use the trapezoidal rule to approximate the double integral in (44). We obtain, up to second order,

µ(I) + µ w (b + h -z) W + tan θ b+h z ϕ(s)ds = 0.
Note that as a consequence we can deduce an explicit expression that is a second order approximation of the analytical solution. We obtain the following approximation,

ϕ = ϕ eq c = ϕ stat c + C 2 (b + h -z) -(b + h -h * ) , v(z) C 1 = C 2 4 (b + h -z) 4 -(b + h -h * ) 4 - 2C 2 (b + h -h * ) -ϕ stat c 3 (b + h -z) 3 -(b + h -h * ) 3 - ϕ eq c -C 2 (b + h -h * ) 2 (b + h -h * ) (b + h -z) 2 -(b + h -h * ) 2 ,
for z > h * , where

h * = (b + h) - tan θ -µ s µ w W, C 1 = µ w (ρ s -ρ f ) g cos θ K 1 η f W , and C 2 = K 2 µ w K 1 W .
Note that in this case, h * is the position of the static/flowing interface.

Numerical approximation: computation of the linearized solid pressure

We approximate the system of ODEs for uniform flows (35) (for the two cases described in section 3.4.1) using a standard first order semi-implicit scheme to go from a time t n to t n+1 = t n + ∆t as follows:

• The first equation is trivially discretized as

h n+1 = h n + ∆tG n f,N + 1 2 .
• For the equation (35b) for the solid volume fraction, we consider an explicit approximation

ϕ n+1 α = ϕ n α -∆tϕ n α Φ n a .
We can also use an implicit discretization or the exact solution of (35b). We have chosen the explicit one for the sake of simplicity, since similar results are obtained with all of them.

• For the momentum equations (35c) and (35d) we use a semi-implicit scheme in order to avoid the restriction on the time step that the viscous terms usually involve:

ρ s l α h n+1 ϕ n+1 α v n+1 α = ρ s l α h n ϕ n α v n α + l α h n β α n (u n α -v n α ) -(ρ s -ρ f ) l α h n ϕ n α g sin θ + ρ s G n s,α+ 1 2 v n α+ 1 2 -ρ s G n s,α-1 2 v n α-1 2 + 1 2 η n s,α+ 1 2 v n+1 α+1 -v n+1 α l α+ 1 2 h n - 1 2 η n s,α-1 2 v n+1 α -v n+1 α-1 l α-1 2 h n , ρ f l α h n+1 (1 -ϕ α ) n+1 u n+1 α = ρ f l α h n (1 -ϕ α ) n u n α -l α h n β α n (u n α -v n α ) + ρ f G n f,α+ 1 2 u n α+ 1 2 -ρ f G n s,α-1 2 u n α-1 2 + η f u n+1 α+1 -u n+1 α l α+ 1 2 h n -η f u n+1 α -u n+1 α-1 l α-1 2 h n .
For α = 1 and α = N the definitions of K α+ 1 2 must be taken into account. The new values are obtained solving a tridiagonal linear system with N equations and unknowns.

If the side walls friction (42) is considered, the term

- 2 W µ w h n α p n s,α v n+1 α |v n α | 2 + δ 2
must be added to the right hand side of the momentum equation for the solid phase, where p s,α = p s,α+

1 2 + p s,α-1 2
/2 is the pressure in the midpoint of layer α.

Once we have h n+1 , ϕ n+1 α , (hϕ α v α ) n+1 , (h (1 -ϕ α ) v α ) n+1 we find u n+1 α , v n+1 α . The difficulty of the scheme comes when computing the solid pressure at the interfaces p n s,α-1 2 at time t n (see [START_REF] Roux | Texture-Dependent Rigid-Plastic Behavior[END_REF] and (30b)):

p n s,α-1 2 = (ρ s -ρ f ) g cos θ N γ=α h n γ ϕ n γ + N ξ=α β n ξ h n ξ ϕ n ξ 1 -ϕ n ξ 2   ξ-1 γ=1 ϕ n γ h n γ Φ n γ + ϕ n ξ h n ξ Φ n ξ 2   , (45a)
depending on the dilatation function

Φ n β = γn β K ϕ n β -ϕ stat c + K 2 I n β , for β = 1, . . . , N. (45b)
In previous equation the inertial number must be approximated in the middle of the layer. For the sake of simplicity in the computation of I n β in (45b), we approximate p n s,α by the value in the interface below. Thus we take

I n β = η f γn β p n s,β-1 2 , (45c)
Note that this approximation follows the same idea of what it is made in single-layer models, see for example [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF]. In this paper the basal pressure is used to define the dilatancy relation arguing that the dilatancy at the bottom gives the right order of magnitude of the dilatancy inside the layer. Here, this hypothesis is adopted for each layer β, thus being p n s,β-

. Therefore, we obtain a nonlinear system with N equations and unknowns (the solid pressures), where the main difficulty is the fact that the system is fully coupled, i.e., each equation depends on all the rest, and this is due to the excess pore pressure.

The first attempt that we carried out was to linearise the problem embedding into (45b) the pressure at the previous time step p n-1 s,β-1 2 . However, due to the abrupt time variations of the excess fluid pressure, namely at short times, this method is inadequate. We propose to consider a numerical method that is based on an approximation of the function f (x) = 1/x by the first-order term of its Taylor polynomial. Thus, in (45b) we write

1 p n s,β-1 2 = 2 p n-1 s,β-1 2 - 1 p n-1 s,β-1 2 2 p n s,β-1 2 + O(p n s,β-1 2 -p n-1 s,β- 1 2
) 2 , for β = 1, . . . , N.

Thanks to this approximation, and defining the coefficients

A n α = 150ϕ n α h n α η f d 2 s (1 -ϕ n α ) 3 , B n α = ϕ n α h n α γn α K ϕ n α -ϕ stat c , C n α = 2 ϕ n α h n α ( γn α ) 2 KK 2 η f p n-1 s,α-1 2 , D n α = C n α 2p n-1 s,α-1 2
, the previous nonlinear system becomes, up to second order, the following N × N linear system

   a n 11 a n 12 • • • a n 1N a n 21 a n 22 • • • a n 2N • • • • • • • • • • • • a n N 1 a n N 2 • • • a n N N         p n s, 1 2 p n s, 3 2 • • • p n s,N -1 2      =    b n 1 b n 2 • • • b n N    ,
where

a ij = D j N k=j+i A k + A j 2 , if i < j, a jj = 1 + D j N k=j+i A k + A j 2 , a ij = D j N k=i A k , if i > j,
and

b j = (ρ s -ρ f ) g cos θ N ξ=j h ξ ϕ ξ + N ξ=j A ξ ξ-1 γ=1 (B γ + C γ ) + A ξ 2 (B ξ + C ξ ) ,
that can be solve using a direct method.

Numerical tests for uniform flows

In this section we present some numerical results in order to validate our model. In this work we only deal with uniform immersed flows and consider dilatancy effects by starting with initially loose or initially dense configurations. With the purpose of comparing the results of the proposed model with model in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF], we consider the rheology given by ( 7), i.e., µ(I) = µ s + K 1 I, in the immersed configuration. First, we compare the results obtained for the proposed model with the analytical solution with and without side walls friction. For sake of simplicity, we only use in this case the model preserving the total height (PH model). The reason is that it is more simple since the final height is known and therefore we can directly compute the velocity and the solid volume fraction profiles. In a second step, we show that the method proposed to approximate the pressure in the previous section 3.4.4 makes it possible to recover the appropriate root of the non-linear system defined by (45), for the particular cases N = 1, 2. However, the error made when approximating the pressure increases with the number of layers, so we decide to use the PGM model with only 2 layers to compare with experimental data (see section 4.2.1). Then, we briefly analyse the influence of the dilatation constant K (see [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF]) and compare the results of the 2 layer PGM model with previous depth-averaged single-layer models in the literature, namely the models introduced in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF][START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] and with laboratory experiments in [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF]. For this comparison we use the model preserving the total solid mass because this is the only physically relevant model.

For the tests presented in this section we set the same values for the parameters than in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF], i.e., the ones proposed in [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] for low and high fluid viscosity. For each case a dense and a loose initial configuration is simulated. The physical and rheological parameters that are common for these tests are detailed in Table 1. The specific parameters in the high or low viscosity cases and for the loose or dense initial conditions are shown in Table 2. In all the presented tests, we consider no-friction with the bottom for the fluid phase, k b = 0, and a no-slip condition (see [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes with a µ(I)-rheology[END_REF]) for the solid one. At the mixture interface, no friction between the fluid in the mixture layer and the fluid in the upper layer is also assumed, i.e., k i = 0.

ρ s (Kg/m 3 ) d s (µm) µ s ϕ stat c K K 1 K

Comparison with analytical solution

In this test we focus on steady uniform flows, where

G f,N + 1 2 = G s,N + 1 2
= 0. Therefore, once the final height is computed (which will be different from the proposed PGM model and from the PH model described in subsection 3.4.1), the analytical solution is the same for both models PH model and PGM model. For the sake of simplicity, we only consider in the comparison with the analytical solution and the PH model, which is simpler than the PGM model since the height is prescribed by the initial condition.

We simulate a uniform flow initially at rest (u α = v α = 0, for α = 1, . . . , N ), with height h 0 , and initial solid volume fraction ϕ α = ϕ 0 for α = 1, . . . , N , until a steady state is reached.

In the low viscosity case we use 50 vertical layers and a time step ∆t L = 10 -5 while 20 layers are used in the high viscosity case and ∆t H = 10 -6 . Note that ∆t H = ∆t L /10, i.e., we need a smaller time step for the high viscosity case, and therefore 20 layers are used in that case because the time needed to complete the simulation is huge. Concretely, we need 35.6 minutes to simulate t f = 300 s for the low viscosity case with 20 layers, whereas the high viscosity case takes 28.9 hours to complete t f = 1500 s of simulation. Finally, the low viscosity case with 50 layers takes 5.1 hours. These CPU times corresponds to an Euler semi-implicit time discretization (see section 3.4.4) and have been measured on a Mac Mini with Intel R Core TM i7-4578U and 16 GB of RAM. Remark that although some differences can be found in the transient regime, the steady solution is the same in both cases (20 or more layers).

No side walls friction

Firstly, we compare our results with the analytical solution for a flow without side walls friction. In Figure 4 we show that the model approximate properly, with loose and dense initial configuration, the analytical maximum velocity (the one on the top of the flow) and the averaged velocity v = N α=1 l α v α . We can see that the convergence in time is much slower in the high viscosity case than in the low viscosity case. This fact, together with the use of a smaller time step justify the use of 20 vertical layers instead of 50 as in the low viscosity case. The main advantage of multilayer models is that the vertical structure of the flow can be recovered. We show these profiles and the analytical profile of velocity in Figure 5, where the steady states reached for the loose/dense initial configuration are the same. In these figures we only show the velocity of the solid phase because we obtain that both, the solid and the fluid velocities, are equal with a relative error of order 10 -4 .

In Figure 6(a) we show the convergence in time of the averaged solid volume fraction ( φ = N α=1 l α ϕ α ) and the equilibrium concentration ( φeq c = N α=1 l α ϕ eq c,α ) towards its steady state. In the inset in Figure 6(a.2) we also show that the solid volume fraction go to the equilibrium concentration quickly (we don't show the loose case because the dynamics towards equilibrium is even quicker in that case). Figure 6(b) shows the normal profile of solid volume fraction. In this case with no lateral friction, when the friction between the phases β α (u α -v α ) is not considered in system (35), and the analytical profile of the solid volume fraction is constant. This is properly reproduced by our model except for the first layer (the closest to the bottom). Moreover, when increasing the number of vertical layers, we obtain the same behaviour. This loss of accuracy in is due to the friction between the solid and fluid phases. Actually, we have checked that if the term βα (u α -v α ) is neglected in the momentum equations then the constant profile of solid volume fraction is obtained, (red-crossed solution in 6(b)). In addition, a convergence test in the number of vertical layers has been performed in this case. We see that the method is second order accurate for the velocity (see Table 3 and Figure 8). For the solid volume fraction, we always obtain an estimated error of 4.63 × 10 -8 , which confirms the fact that the friction between the solid and the fluid phases is the responsible for not obtaining the constant profile in Figure 6(b). Figure 7 shows the pressure and the excess pore pressure at the bottom. We see that the numerical solutions go to the analytical one in all the configurations, and that there is no excess pore pressure in the steady state.

N (∆z = 1/N ) L 1 -Error L 1 -Order L 2 -Error L 2 -Order L ∞ -Error L ∞ -Order 2 1.25×10 -1 - 1.17×10 -1 - 9.09×10 -2 - 4 3.12×10 -
Figures 9 and10 show the normal profiles of excess pore pressure at different times. We see that the profiles are approximately linear at short times, when the excess pore pressure is larger, except near the bottom.

Side walls friction

In subsection 3.4.3 the analytical solution is obtained in the case of a flow that is confined in a channel whose width is W . In that case, the side walls friction has an important effect on the dynamics of the flow. Mainly, the normal profile can become S-shaped instead of Bagnold type (see [START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: a well-balanced multilayer discretization[END_REF]) leading to the appearance of a flow/no-flow transition. We consider a flow with the same physical and rheological properties as in previous subsection in the low viscosity regime and loose initial condition, and we add the side walls friction term with the friction coefficient µ w = tan (13.1 • ). In this case we use 20 layers in the multilayer system since the computational cost, until reaching the stationary solution, is very high.

Figure 11 shows the normal profiles of velocity and solid volume fraction in the presence of lateral walls, for several channel widths W = 0.005, 0.0075, 0.01, 0.012, 0.015, 0.02, 0.03, 0.05, 0.1, 1, 10 m. These values correspond to α layers of granular particles, W = αd s , with α = 31, 47, . . . , 62500, approximately. For the velocity profiles, the dynamics is similar to the case of a dry granular flow (see [START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: a well-balanced multilayer discretization[END_REF]). In the case of the solid volume fraction, a constant profile is obtained for large widths, similar to the case of no lateral wall friction.

However, looking at the analytical solution (see remark 2) we see that the analytical profile of the volume fraction is always linear, starting from the equilibrium volume fraction ϕ eq c at the top of the mixture layer. Actually, the linear counterpart is K 2 µ w z/(K 1 W ). Therefore, when W is large, this contribution is nearly zero and we obtain an (almost) constant profile. As observed on Figure 11(b), wall friction makes the solid volume fraction vary linearly with a slope increasing when the channel width W decreases. An interesting remark is that when a flow/no-flow transition appears in the velocity profile, the solid volume fraction becomes a constant, whose value is ϕ stat c = 0.582, from the flow/no-flow position to the bottom representing a region of rigid material (zero velocity). It decreases from the flow/no-flow transition up to the upper surface of the mixture. In Figure 11 we also see the influence of the friction parameter µ w . For widths W = 31d s , 47d s we have simulated the same flow but taking as friction coefficient at the lateral walls µ w = 2µ w . It makes the friction force increase (similar effect is obtained by decreasing the channel width) and the solid volume profile becomes constant almost everywhere in the normal direction, decreasing just close to the mixture upper surface (see blue symbols and red dash-dotted lines in Figure 11(b)).

In Figure 11 we remark that errors between the analytical and the numerical solutions are bigger than the expected ones, specially for the solid volume fraction and a small channel width (see e.g. red dash-dotted lines in Figure 11(b)). This is due to the fact that we have used 20 layers to reproduce steady state solutions with a strong variations in the normal direction and the multilayer approach is based on a constant profile in each layer. In order to clarify that, we have performed a convergence test for W = 94d s , also neglecting the term β α (u α -v α ), since it is not considered for the analytical solution. We see that the method is first order accurate in this case, for both, the velocity and the solid volume fraction (see Table 4 and Figure 12). In addition, for the two configurations for which the side wall friction is stronger, namely W = 31d s , 47d s , we have also added the corresponding solutions using 160 layers in Figure 11(b) (cyan dash-dotted lines). The analytical solution is properly reproduced in that case. (a) Velocity

N (∆z = 1/N ) L 1 -Error L 1 -Order L 2 -Error L 2 -Order L ∞ -Error L ∞ -Order 2 1.94×10 -1 - 1.84×10 -1 - 1.74×10 -1 - 4 1.45×10 -1 0.47 1.42×10 -1 0.37 1.34×10 -
N (∆z = 1/N ) L 1 -Error L 1 -Order L 2 -Error L 2 -Order L ∞ -Error L ∞ -Order 2 1.57×10 -2 - 1.63×10 -2 - 1.97×10 -2 - 4
9.79×10 -3 0.68 1.00×10 In Figure 11 we have used 20 layers in the simulation in order to properly reproduce complex vertical profiles of velocity and solid volume fraction. However, if the number of layers is reduced, although the profile is not exactly reproduced, its shape and mean variation is well ap- proximated. In Figure 13 we show these vertical profiles for two representative channel widths W = 62d s , 62500d s , computed with the multilayer model with 2, 3 and 20 vertical layers. In the first case, the side walls friction has a strong effect on the vertical profiles, leading to a rigid zone with zero velocity close to the bottom. We see that this solution is well reproduced with 20 layers. The models with 2 and 3 layers give also a good approximation. Even the position of the flow/no-flow interface is also reasonably well approximated. The mean variation of the profile is also even better reproduced if the side walls friction is smaller, as shown for W = 62500d s .

Comparison with other models and laboratory experiments

In this section we compare the results of the PGM model with the results of the model in Bouchut et al. (2016) [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] and Pailha and Pouliquen (2009) [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] (hereafter B-2016 and PP-2009 respectively). However, before that, we evaluate in the next subsection the appropriate model that we propose to use in practice, since the approximation of the solid pressure entrains serious difficulties in the multilayer framework, as we mentioned before.

Choice of the model

In this section the simulated cases, as well as the physical and rheological parameters, take the same values as in the previous section. We focus on the dense case for both the high and low viscosity cases because it is the more complex case, since the convergence to the steady solution is slower and it differs far more from other analyzed models, as it is shown in following subsections. Firstly, we check that the proposed linearization to compute the pressure gives good results. Note that equations (45) lead to a non-linear system with N equations and unknowns, where each of them is a quadratic equation, and the system is fully coupled due to the excess pore pressure. Therefore, for a configuration with N vertical layers, there exist 2 N vector solutions of the system. In particular, for N = 1, at each time step, once the rest of variables are computed, the pressure at the bottom, p n s, 1 2 , is a root of the quadratic equation resulting from (45) with N = 1. That equation can be exactly solved, obtaining two different roots for p n s, 1 2 . However, one of them becomes negative and therefore there is only one solution with a physical meaning. In Figure 14(a)(b) we show that the proposed method in 3.4.4 gives as result an approximation to this appropriate solution. Analogously, for N = 2, system (45) has 4 different solutions for the vector p n s, 1 2 , p n s, 3 2 at each time step, corresponding to the pressure at the interfaces z 1 2 and z3

2

. Again, the system obtained by (45) with N = 2 can be exactly solved, where just one of the roots is relevant, since the others become negative. Figure 14(c)(d) shows that the proposed method approximates this solution. Nevertheless, we did not managed to calculate an accurate approximation of this solution at short times where the variation and effect of the excess pore pressure are strong.

Our results show that it is very difficult to calculate the excess pore pressure when we have a large number of layers because the system is fully coupled, i.e. each pressure depends on all the pressures in the other layers. In figures 15 and 16 we see a comparison of the velocity at short times with the multilayer model with 2, 3, 4, and 10 layers, and the velocity measured at the surface of the mixture in the laboratory experiments in [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF], for the low and high viscosity cases and the dense initial configuration.

We observe that the model with 2 or 3 layers better reproduces the experimental results than the model with 10 layers, in particular at the beginning, when the velocity is very small. The 10 layers model predicts that the velocity grows up too fast in both cases, with low and high viscosity. We see that the single-layer model B-2016 better predicts the time when the velocity starts to grow up, while the model with 2, 3 layers captures the maximum (i. e. surface) velocity much better than the single layer model.

We conclude that the lack of accuracy at short times when using a large number of layers is due to the approximation of the pressure detailed in subsection 3.4.4 and the fact that the non-linear system to solve is strongly coupled, i.e., the errors approximating the pressure in each layer accumulates within the whole domain. A more accurate method to approximate the pressure should be developed if more than 2-3 layers are used. As a result, we will mainly use the 2-layer model in the following. It is also justified looking at Figure 13, where we have shown that, despite of the fact that we need many layers to give an accurate approximation of the vertical profiles of velocity and solid volume fraction, the models with 2 and 3 layers are able to capture surprisingly well the mean variation of these profiles.

Figures 15 and16 show the influence of the dilatation constant K (see [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF]), which has been calibrated for single-layer models. If we compare the results obtained with two values of K by showing also the results obtained with K = 2K (dashed blue lines), we observe a strong sensitivity of the results to this parameter. This is the case in particular for the high viscosity system (see the case with 10 layers in Figure 16) for the velocity behaviour in particular at the beginning when the velocity starts to grow up. These results suggest that this constant should be calibrated for 3D models (or at least for multi-layer cases).

Therefore, in the following we keep only the model with two layers (Preserving Granular Mass two-layers model, denoted PGM-2 model hereafter) because it is the simplest case giving an appropriate solution, although we could choose also the one with 3 layers. In appendix B, we write this particular case for the reader that is interested on using this particular model, and not on the general multilayer case. In the following subsection we finally compare this two-layer model with previous models in the literature and experimental data.

PGM-2 model vs other models

We compare here the results of the proposed two-layer model (PGM-2 model) with the results of models B-2016 and PP-2009. These depth-averaged models only compute the averaged velocity and solid volume fraction, while the PGM-2 model allows us to compute two values of these quantities in the normal direction. In particular, in the configuration considered here, the velocity in the bottom layer is lower than the velocity in the top layer, further called the maximum velocity. As previously, the physical and rheological parameter are fixed in tables 1 and 2.

In figures 17 and 18 we see the time evolution of the maximum and the averaged velocities, the pressure at the bottom, the concentration, the solid mass and the height for all the models in both, the loose and dense initial configuration, for the high and low viscosity cases. All the models give similar results at the steady state although the time evolution can differ. We see that the results of the proposed PGM-2 model and the B-2016 model are very close for the pressure, averaged solid volume fraction and height, whereas the velocities are slightly different. The behaviour is very similar, for example, the steady velocity for the dense initial configuration is greater that the one for the loose configuration. However, with the 2-layers model the obtained maximum velocity in the top layer is much closer to the velocity measured at the surface of the mixture in the experiments (see Figure 21). Therefore, the main advantage of this model with respect to the previous depth-averaged models B-2016 and PP-2009 is the approximation of the velocity field.

Both PGM-2 and B-2016 models preserve the granular mass Figures 19 and20 show the convergence of the solid volume fraction to the equilibrium concentration in the high and low viscosity cases. This convergence is much slower in the dense configuration than in the loose one, as commented before. In fact, we can see in Figure 19b that for the loose configuration the solid volume fraction φ goes to φeq c at very short times (t ∼ 10 -3 ). In Figure 21 we show the velocity and excess pore pressure with models PP-2009, B-2016 and PGM-2, and also the laboratory experiments in Pailha et al. 2009 [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF], in both the high and low viscosity case. Looking at the averaged velocities computed with the depth-averaged models, we see that they cannot predict the velocity observed at the mixture surface in the experiments. However, the maximum velocity computed by the proposed two-layers model has the order of magnitude of the surface velocity in the experiments. On the other hand, having two layers instead one layer does not significantly change the calculated excess pore pressure. As expected, this excess pore pressure goes to zero after some time in all models.

Finally, in Figure 22 we see the evolution of the fluid transference, G f,N + 1 2 , at the interface between the mixture layer and the upper fluid layer. We recall that G f,N + 1 2 have the opposite sign to V f in previous work [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] (see [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF]). For the high viscosity case ( Figure 19(b)) we see very clearly the dynamics. When the mixtures starts to flow in the initially dense case, we have ϕ > ϕ eq c , leading to a positive dilatation angle tan ψ = K (ϕ -ϕ eq c ) > 0 and therefore to dilatation. As a consequence the fluid is sucked into the mixture (V f = -G f,N + 1 2 < 0) and the height of the mixture increases (see inset in Figure 17c). On the contrary, in the initially loose case, ϕ < ϕ eq c , leading to tan ψ = K (ϕ -ϕ eq c ) < 0 and therefore to contraction. The fluid is thus expelled from the mixture (V f = -G f,N + 1 2 > 0) and the mixture height decreases. Focusing now on the low viscosity case, we observe the same behaviour for the initially dense case, but not for the loose case. Figure 23 shows the loose case at very short times, where the material contracts at the beginning as expected, i.e. the fluid is expelled (-G f,N + 1 2 > 0) and the height decreases. However we see that after some time ϕ 2 < ϕ eq c,2 (contraction) in the top layer (layer 2) whereas ϕ 1 > ϕ eq c,1 (dilatation) in the bottom layer (layer 1). Moreover, it makes that the fluid transference at the top changes its sign and becomes negative at time t ≈ 5.1 × 10 -3 , and Physically, it could mean that when the fluid is expelled at the beginning, the lower part is the first one that remain with less fluid (because each layer is expelling fluid toward the upper layer, but receiving it from the layer below), and therefore it becomes to a dense situation. After that, the lower part dilates and the upper part contracts, going to a globally dilatation situation, as we can see in Figure 23.

Conclusions

In this work we deal with dilatancy effects in granular models of fluidized flows. A multilayer extension of the two-phase model proposed by [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] has been presented. It is derived from a dimensional analysis and the multilayer approach following an analogous procedure to [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF]. This model is in principle able to recover the structure of the fluid in the direction normal to the topography. This is essential to improve the description of friction and viscous terms when considering a complex rheology compared to depth-averaged single-layer models, as showed in [START_REF] Fernández-Nieto | A multilayer shallow model for dry granular flows with the µ(I) rheology: Application to granular collapse on erodible beds[END_REF][START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: a well-balanced multilayer discretization[END_REF]. In particular, it allows to recover the normal velocity profiles, but also the normal profiles of the solid volume fraction and pore fluid pressure.

The proposed two-phase model behaves as the B-2016 [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] model, preserving the total solid mass. In addition, a model that behaves as the PP-2009 model [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] -preserving the total mixture height-has been also obtained for a particular closure relation for the solid mass transference at the mixture upper interface. That closure relation allows us to obtain either a model (PH) preserving the total height (and maybe loosing solid mass) or a model (PGM) preserving the total solid mass (the height vary depending on the contraction/dilatation of the solid). This last model is the physically meaningful model that we further compare with laboratory experiments.

The main difficulty of the proposed two-phase model is the excess pore fluid pressure that appears in particular in the solid pressure, which is no more hydrostatic even with the thin layer approximation. This pressure varies strongly in time when starting from non-equilibrium conditions due to dilatation or contraction of the granular medium. It further goes to zero when reaching steady states. Dealing with this strongly varying, non-linear term is shown to be a real numerical challenge to approximate the model when variations in the normal direction are considered. For this reason, we only deal with uniform flows in the numerical tests in this work, although the model is deduced for general flows, and simulation of x-dependent flows could be part of a forthcoming work. Even in the case of uniform flows, where the model is reduced to a system of ODEs, the numerical approximation is not easy and a specific numerical scheme is necessary. As consequence of this pore excess pressure term, the pressure cannot be easily computed, even in the single-layer case where the pressure at the bottom is one of the two roots of a quadratic equation. A numerical scheme based on the first-order Taylor polynomial has been proposed to approximate the pressure.

We have validated the model by comparing with the analytical solution for two-phase uniform flows. In this case we have chosen the PH model because computing the analytical solution in that case is simpler than for the PGM model, since the height is directly obtained. We have seen that the steady solution does not depend on the dense/loose initial value of the solid volume fraction for the PH model (also PP-2009 model), where the solid mass varies in time, in contrast with the results obtained with the PGM model (also B-2016 model), where the solid mass is preserved and the steady states for the velocity and pressure depend on the initial solid volume fraction. For these tests, a constant profile of solid volume fraction has been obtained. Furthermore, an analytical solution for confined flows, where side wall friction notably changes the dynamics of the flow, has been deduced from the resolution of a integro-differential equation. In that case, the velocity profile goes from a Bagnold profile to a S-shaped profile, very similar to the case of dry flows (see [START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: a well-balanced multilayer discretization[END_REF]). It is interesting that the profile of solid volume fraction is not constant in this case, but linear, and when a flow/no-flow transition is observed in the velocity profile, then the solid volume fraction becomes constant in the no-flow (i. e. rigid) zone. In addition, we see that including side wall friction as in Jop et al. [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF] is a good approximation for shallow flows. An interesting feature of this approach is that is allows to give an approximation of flow/no-flow transition which is obviously not the case for depth-averaged single-layer models. Furthermore, for the analytical solutions (with and without side wall friction) we have also shown convergence tests depending on the number of vertical layers. In particular, for solutions with a strong vertical structure (the case of a narrow channel width) the analytical solution is perfectly recovered when increasing the number of layers.

We have shown that the approximation of the pressure with the proposed method is acceptable, in the sense that it approximates the good solution from the 2 N possible solution vectors for the pressure. However, this approximation is not accurate enough at short times and the error notably increases with the number of layers in the normal direction. We have seen that the model with 2 or 3 layers gives reasonable results when compared with experimental data, allowing to make the quantities (velocity, volume fraction, pore pressure, etc.) vary in the normal direction. The pore pressure estimates however is getting less and less accurate when increasing the number of layers due to the accumulation of the numerical error for this full coupled models. This problem will be encountered also if solving the complete Jackson's equation for any 3D method. Another conclusion of this work is that the dilatation constant should be fitted for 3D models, since it has a strong influence on the results, in particular when using a large number of layers.

The model with two layers (PGM-2 model) has been chosen for further comparison with experimental results. We conclude that a more accurate method to compute the pressure is necessary in order to approximate the solution at short times. An interesting alternative could be to use an adaptive multilayer method, which uses few layers at short times giving a good approximation of the solid pressure, and increases the number of layers when the excess pore pressure becomes small. That would allow us to get a better approximation of the velocity profile at latter times and will be investigated in the future.

We compare the results of the PGM-2 model with the model of [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] [6] (B-2016) and [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] (PP-2009). We see that our results are similar to the B-2016 model, except for the velocity. Actually, these previous depth-averaged models were not able to predict the velocity measured at the mixture surface, while the maximum velocity calculated at the top layer with the PGM-2 model is closer of this observed velocity. More differences are found when comparing the proposed model with the PP-2009 model as expected, since this model preserves the total height instead of the mass and therefore the steady states of the loose/dense initial condition are the same, in contrast with models B-2016 and PGM-2 models. and using the expression for the solid and fluid mass transference terms (34), we get

ε β α (1 -ϕ α ) z α+ 1 2 -z (u z α-1 2 ) + -(v z s,α-1 2 ) + = -β α ϕ α (1 -ϕ α ) 2 z α+ 1 2 -z α-1 γ=1 ϕ γ h γ Φ γ + O(ε).
Note that here we have used the nondimensional expression of the mass transference terms, whose leading order is 1/ε, corresponding to the term involving the dilatancy function. As a result, we get the following explicit expression for the excess pore pressure p e f,α = p e f,α+ In this case, we consider the moderate friction regime by fixing k = 0 in (50), the excess pore pressure results p e f,α = p e f,α+

1 2 -β α h 2 α -z -z α-1 2 2 2 (1 -ϕ α ) 2 Φ α - β α ϕ α (1 -ϕ α ) 2 z α+ 1 2 -z α-1 γ=1 ϕ γ h γ Φ γ + O(ε),
1 2 + ε 2 β α (1 -ϕ α ) z α+ 1 2 -z G s,α-1 2 ϕ α - G f,α-1 2 (1 -ϕ α ) + (u α -v α ) ∂ x z α-1 2 -ε 2 β α (1 -ϕ α ) h 2 α -z -z α-1 2 2 2 (1 -ϕ α ) ∂ x (1 -ϕ α ) (u α -v α ) + 1 ε Φ α ,
where p e f,α+ 1 2 = p e f,α+1 z = z α+ 1 2 , and G s , G f are given by (34). Note that this expression is more complicated because it involves second order derivatives in the momentum equations.

B Preserving Granular Mass two-layers model (PGM-2 model)

This appendix is devoted to presenting the two-layer model, which we have used in section 4.2.2 to compare our results with previous models in the literature: Bouchut et al. (2016) [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] and Pailha and Pouliquen (2009) [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF], for readers that are interested on using this model without look at the detail of the multilayer case.

Figure 1 :

 1 Figure 1: Sketch of the domain. The term V f , which may be no-zero, represents the fluid transference at the mixture interface with the upper fluid domain.

Figure 2 :

 2 Figure 2: Sketch of domain with the particular multilayer division for the mixture grain-fluid layer and its notation.

  Low vis cos ity (b) High vis cos ity

Figure 4 :

 4 Figure 4: Time evolution of the maximum (v max , brown lines) and averaged (v, green lines) velocities in the (a) low and (b) high viscosity case. Dashed (resp. solid) lines are the solution starting from a loose (resp. dense) configuration, and the dash-dotted (resp. dashed) grey line is the analytical averaged (resp. maximum) velocity.

Figure 5 :

 5 Figure 5: Normal profile of downslope velocity in the (a) low and (b) high viscosity case. Green crosses (resp. brown squares) are the solution in the dense (resp. loose) initial configuration, and the dashed grey line is the analytical profile of velocity.

Figure 6 :

 6 Figure 6: (a) Time evolution of the averaged solid volume fraction (green lines) and equilibrium concentration (brown lines). Dashed (resp. solid) lines are the solution starting from a loose (resp. dense) configuration, and the dashed grey line is the analytical averaged solid volume fraction. The inner figure in (a.2) is the difference | φ -φeq c | in the dense initial configuration. (b) Normal profiles of solid volume fraction in the (b.1) low and (b.2) high viscosity case. The dash-dotted grey line is the analytical normal profile of the solid volume fraction. Red-crossed are the simulation neglecting the friction term β α (u α -v α ) in the momentum equation.

Figure 7 :

 7 Figure 7: Time evolution of (a) the pressure and (b) the excess pore pressure at the bottom in the (x.1) low and (x.2) high viscosity case. Dashed (resp. solid) lines are the solution starting from a loose (resp. dense) configuration, and the dashed grey lines are the analytical steady states.

Figure 6 (

 6 Figure 6(b) is due to the friction between the solid and fluid phases. Actually, we have checked that if the term βα (u α -v α ) is neglected in the momentum equations then the constant profile of solid volume fraction is obtained, (red-crossed solution in 6(b)). In addition, a convergence test in the number of vertical layers has been performed in this case. We see that the method is second order accurate for the velocity (see Table3and Figure8). For the solid volume fraction, we always

Figure 8 :

 8 Figure 8: Relative errors of the velocity for the analytical solution without side walls friction and neglecting the friction term β α (u α -v α ).

Figure 9 :Figure 10 :

 910 Figure 9: Normal profiles of excess pore pressure for the low viscosity case, in the (a) loose and (b) dense configurations, at times t = 0.1, 1, 5, 10, 20, 30, 50, 70, 100, 400 s.

Figure 11 :

 11 Figure 11: Normal profiles of velocity (a) and solid volume fraction (b) for different widths W . Dashed lines correspond to the simulations and black crosses to the analytical solutions. Blue symbols and red dash-dotted lines in figure (b) are the simulations for widths W = 31d s , 47d s and µ w = 2µ w , additionally cyan dash-dotted lines are the simulations for this case with 160 vertical layers.

Figure 12 :

 12 Figure 12: Relative errors of (a) the velocity and (b) the solid volume fraction for the analytical solution with side walls friction and W = 94d s .

  (a) the velocity and (b) the solid volume fraction of the analytical solution with side walls friction, W = 94d s , and neglecting the friction term β α (u α -v α ).

Figure 13 :

 13 Figure 13: Normal profiles of velocity and solid volume fraction for two representative widths (a) W = 62d s and (b) W = 62500d s , for different number of vertical layers. Dashed blue lines correspond to the simulations with 20 layers, red circled dash-dotted lines with 3 layers, green square dotted lines with 2 layers and black crosses are the analytical solutions. Notice that symbols for 2 and 3 layers correspond to the values in the middle of each layer.

Figure 14 :

 14 Figure 14: Time evolution of exact (solid lines) and computed (blue symbols) solid pressures for the case of a single layer ((a) and (b)) and two layers ((c) and (d)) in the multilayer system, in the low ((a) and (c)) and high ((b) and (d)) viscosity case with dense initial configuration. Dotted lines in the two layer case are the pressures at the internal interface z = z 3/2 , while solid lines are the pressures at the bottom z = z 1/2 . Inset figures are the long time exact solutions.

Figure 15 :

 15 Figure 15: Low viscosity case with dense initial condition. Time evolution of the velocity of each layer for the multilayer model with 2, 3, 4 and 10 layers. Solid red lines are computed with dilatation constant K = 4.09 and dashed blue lines are with K = 8.18. The black-squared solid line is the solution of the single-layer model B-2016, and grey lines is the velocity measured at the surface of the mixture in laboratory experiments.

N α=1 h α

  ϕ α , contrary to the PP-2009 model for which the total height (h) is preserved (see figures 17(c)(d) and 18(c)(d)). Note that a multilayer model preserving the total height has been also obtained in subsection 3.4.1. Imposing mass conservation leads to different velocity, pressure, and solid mass in steady state for initially loose and dense cases (PGM-2 and B-2016 models), while the steady state reached with PP-2009 model does not depend on the initial dense/loose configuration. These results are the

Figure 16 :

 16 Figure 16: High viscosity case with dense initial condition. Time evolution of the velocity of each layer for the multilayer model with 2, 3, 4 and 10 layers. Solid red lines are computed with dilatation constant K = 4.09 and dashed blue lines are with K = 8.18. The black-squared solid line is the solution of the single-layer model B-2016, and grey lines is the velocity measured at the surface of the mixture in laboratory experiments.

Figure 17 :Figure 18 :

 1718 Figure 17: Time evolution of (a) velocity, (b) pressure at the bottom, (c) averaged solid volume fraction, (d) solid granular mass in the high viscosity case.

Figure 19 :Figure 20 :

 1920 Figure 19: (a) Time evolution of the averaged volume solid fraction and equilibrium concentration for the high viscosity case and (b) time evolution at short times.

Figure 21 :

 21 Figure 21: (a),(c) Time evolution of the averaged velocity; (b),(d) the excess pore pressure at the bottom in the high viscosity case ((a),(b)) and the low viscosity case ((c),(d)). Grey lines are the experiment data.

Figure 22 :

 22 Figure 22: Time evolution of the fluid transference at the interface between the mixture and the upper fluid layer in the (a) low and (b) high viscosity cases, with dense/loose initial condition.

Figure 23 :

 23 Figure 23: Time evolution of the (a) fluid transference at the interface between the mixture and the upper fluid layer, and (b) the solid volume fraction at short times for the low viscosity case and loose initial condition.

2 

 2 h ξ ϕ ξ (1 -ϕ ξ ) 2   ξ-1 γ=1 ϕ γ h γ Φ γ + ϕ ξ h ξ Φ ξ + O(ε).

Table 1 :

 1 Common physical and rheological parameter for the two-phase simulations.η f case ρ f (Kg/m 3 ) η f (Pa/s) θ ( • ) h 0 (mm) ϕ 0 (loose) ϕ 0 (dense)

	2

Table 2 :

 2 Parameter for the two-phase simulations depending on the high/low fluid viscosity and loose/dense initial configuration.

Table 3 :

 3 Order of the error for the velocity of the analytical solution without side walls friction and neglecting the friction term β α (u α -v α ).

		2	2.00	2.87×10 -2	2.02	2.12×10 -2	2.09
	8	7.81×10 -3	2.00	7.14×10 -3	2.00	5.23×10 -3	2.02
	16	1.95×10 -3	2.00	1.78×10 -3	2.00	1.30×10 -3	2.00
	32	4.88×10 -4	1.99	4.45×10 -4	2.00	3.25×10 -4	2.00
	64	1.22×10 -4	1.99	1.11×10 -4	1.99	8.15×10 -5	1.99
	128	3.05×10 -5	1.99	2.79×10 -5	1.99	2.04×10 -5	1.99
	256	7.67×10 -6	1.99	7.00×10 -6	1.99	5.21×10 -6	1.97

Table 4 :

 4 Order of the error for

	-2	0.69	1.18×10 -2	0.74

Acknowledgements

This research has been partially supported by the Spanish Government under grants MTM2015-70490-C2-2-R and RTI2018-096064-B-C22 with the participation of FEDER, by the ERC contract ERC-CG-2013-PE10-617472 SLIDEQUAKES.

A Explicit expression of the excess pore pressure

In this appendix we detail the computation of the excess pore pressure in equation ( 27)

with

Firstly, we integrate the dilatancy equation (13g) between z α- as it is made for the calculation of the vertical velocity, resulting

We can also integrate equation [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF] regarding the incompressibility of the mixture, obtaining

Now, by considering equation (49) minus equation ( 48) and dividing the result by (1 -ϕ α ), we get

Once we have the explicit expression for the difference of the vertical velocities, we have to compute the integrals in ( 46)-( 47). Tanking into account that

we obtain that

In this work we are interested on the case of a strong friction between the fluid and granular phases. This case is more suitable for natural context applications as discussed in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF], and moreover it is simpler than the moderate friction choice, involving only first-order derivatives in the momentum equations. Therefore we take k = -1 in previous equation, leading to

Using the mass transference terms [START_REF] Jackson | The Dynamics of Fluidized Particles[END_REF] we obtain

2