François Bouchut

Enrique D Fernández-Nieto

El Hadji Koné

Anne Mangeney

Gladys Narbona-Reina

Dilatancy in dry granular flows with a compressible µ(I) rheology

Keywords: Dry granular materials, compressible rheology, dilatancy, depth-averaged model, granular collapse, well-balanced scheme

Dilatancy plays a key role in mixtures of grains and fluid but is poorly investigated in dry granular flows. These flows may however dilate by more than 10% in granular column collapses. We investigate here dilatancy effects in dry flows with a shallow depth-averaged model designed to be further applied to simulate natural landslides. We use a compressible µ(I), φ eq (I) rheology with a dilatancy law, where φ eq (I) is the volume fraction at the equilibrium (i.e. critical) state and I the so-called inertial number. This law is obtained by simply removing the fluid phase in the solid/fluid model of our previous work [J. Fluid Mech, 801, 166-221 (2016)] and derived from critical state theory. A numerical method is proposed to solve the equations, that have however singularities that are rather difficult to handle.

Simulations of granular collapses on horizontal and sloping beds show that the maximum height of the deposits changes as a function of the initial volume fraction with higher (lower) deposits for initially denser (looser) granular masses, as observed with Discrete Element simulations. The front position and the deposit shape behind it are on the contrary poorly affected by the initial volume fraction, as if the flow had almost forgotten its initial state. However subtle effects can be observed with the occurrence of low velocity regimes on steep slopes that strongly depend on the initial volume fraction. Simulations show complex compression/dilation effects during the flow, in particular with front dilation (compression) during the acceleration (deceleration) phases. These effects may dramatically change the effective friction that is observed to decrease at the front in some situations, while the µ(I) rheology without dilatancy would have predicted an increasing friction. The model predicts an increasing dilation of the mass for increasing slopes by up to 10% in the studied configurations, in agreement with laboratory experiments. Our results suggest that this compressible model contains key features to describe granular dilatancy.

Introduction

Dilatancy plays a key role in mixtures of grains and fluid, making it possible to dramatically change their flow dynamics and deposit [START_REF] Iverson | The perfect debris flow? Aggregated results from 28 large-scale experiments[END_REF][START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF][START_REF] Bouchut | A twophase solid-fluid model for dense granular flows including dilatancy effects: comparison with submarine granular collapse experiments[END_REF][START_REF] Delannay | Granular and particleladen flows: from laboratory experiments to field observations[END_REF]. This process is an essential ingredient in landslide numerical models, increasingly used to build hazard maps for risk assessment [START_REF] Kean | Debris-flow hazards mitigation: mechanics, monitoring, modeling, and assesssment[END_REF].

Dilatancy can be described using the critical state theory [START_REF] Schofield | Critical state soil mechanics[END_REF]. The essential feature is that a granular material, when sheared, reaches an equilibrium state called the critical state. When the volume fraction of the granular mass is initially smaller than the equilibrium volume fraction (initially loose material), the granular material compresses under shear. This compression expels the interstitial fluid, increasing the pore pressure and therefore decreasing the effective friction experienced by the grains. This effect may possibly lead to the so-called liquefaction of granular mass. The contrary is observed for initially dense materials that dilate under shear, leading to a decrease of fluid pressure and an increase of effective friction.

Dilatancy has been introduced in models of grain/fluid mixtures by using the law proposed by [START_REF] Roux | Texture-dependent rigid plastic behavior[END_REF] for dry granular material [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][START_REF] Meng | Modeling dynamic flows of grain-fluid mixtures by coupling the mixture theory with a dilatancy law[END_REF]. Depth-averaged models based on the shallow approximation (the depth of the landslide is assumed to be much smaller than its extension along the slope) are extensively used to simulate granular flows in natural geological context. Indeed, these approximations reduce the high computational cost related to the need of accurately describing the complex topography (e. g. [START_REF] Pirulli | Result of back-analysis of the propagation of rock avalanches as a function of the assumed rheology[END_REF][START_REF] Moretti | Numerical modeling of the mount Meager landslide constrained by its force history derived from seismic data[END_REF][START_REF] Lucas | Frictional velocity-weakening in landslides on earth and on other planetary bodies[END_REF][START_REF] Mcdougall | A model for the analysis of rapid landslide motion across three-dimensional terrain[END_REF][START_REF] Iverson | Modelling landslide liquefaction, mobility bifurcation and the dynamics of the 2014 oso disaster[END_REF]). Solving numerically two-phase equations with dilatancy in such models however involves several difficulties: shear rate and pore fluid pressure dependency in the dilatancy law that induces several singularities at low depth, tendency to produce oscillations on the volume fraction due to the coupling between convection and dilatancy, nontrivial steady states to preserve [START_REF] Bouchut | A twophase solid-fluid model for dense granular flows including dilatancy effects: comparison with submarine granular collapse experiments[END_REF].

Even though dilatancy is expected to have a much smaller impact on dry granular flows than in grain/fluid mixtures, the volume of the flowing mass may increase by more than 10% during granular column collapses over sloping beds [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF]. The effect of dilatancy in dry granular flows has been poorly investigated. We can found some works where dilatancy has been tackled in dry granular flows [START_REF] Lee | A three-dimensional continuum model incorporating static and kinetic effects for granular flows with applications to collapse of a twodimensional granular column[END_REF][START_REF] Si | A general frictional-collisional model for dense granular flows[END_REF]. In particular, dilatancy is considered in the closure of the stress tensor as constitutive relations. Nevertheless, these authors did not investigated the role of dilatancy in their results. The question is how and how much dilatancy changes flow dynamics, effective friction and deposit shape of dry flows. Recent simulations based on Discrete Element Methods (DEM) showed that the deposit of dry granular collapses on an horizontal plane changes as a function of the initial state (loose or dense) of the released mass [START_REF] Kermani | Simulation of Collapse of Granular Columns Using the Discrete Element Method[END_REF]. In particular the maximum thickness of the deposit is significantly smaller for initially loose than for initially dense granular columns. On the contrary, the front position and the deposit shape near the front are poorly affected by the initial volume fraction as if the material had forgotten the memory of its initial state in the front zone. Other DEM simulations and laboratory experiments on sloping topography suggest that the front travels further for initially dense masses [START_REF] Cagnoli | Grain size and flow volume effects on granular flow mobility in numerical simulations: 3-D discrete element modeling of flows of angular rock fragments[END_REF]. Although physically relevant, Discrete Element Methods are hardly applicable to natural flows due to their huge computational cost and to the unknown particle size distribution in real landslides, that has been shown to significantly affect simulation results [START_REF] Cagnoli | Grain size and flow volume effects on granular flow mobility in numerical simulations: 3-D discrete element modeling of flows of angular rock fragments[END_REF][START_REF] Borykov | Investigation of frictional weakening of martian and terrestrial landslides using discrete element simulations[END_REF].

We propose here to investigate these effects by including dilatancy in a shallow depthaveraged model of dry granular flows with the ultimate objective to further apply it to simulate real landslides. One of the question is as to whether the simple law proposed by [START_REF] Roux | Texture-dependent rigid plastic behavior[END_REF] makes it possible to qualitatively and possible quantitatively reproduce these effects. The frictional behavior in the model will be described by the µ(I) rheology where I is the so-called inertial number. This widely used rheology has been shown to reproduce a wide range of laboratory experiments (see references in [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Delannay | Granular and particleladen flows: from laboratory experiments to field observations[END_REF]) even though it is ill-posed for some values of I [START_REF] Barker | Partial regularisation of the incompressible µ(I)-rheology for granular flow[END_REF][START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF]. It has recently been shown that compressible effects may resolve the ill-posedness of the initially incompressible µ(I) rheology [START_REF] Barker | Partial regularisation of the incompressible µ(I)-rheology for granular flow[END_REF][START_REF] Heyman | Compressibility regularizes the µ(I)rheology for dense granular flows[END_REF]. However, when describing compressibility only by adding a variable volume fraction φ(I) (so-called µ(I), φ(I) rheology), the model is even more prone to ill-posedness in time-dependent calculations than the incompressible µ(I) rheology [START_REF] Heyman | Compressibility regularizes the µ(I)rheology for dense granular flows[END_REF][START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF]. A more realistic way to introduce compressible effects such as dilatancy is to use, in addition to the µ(I), φ(I) rheology, flow rules derived from critical state theory as that proposed by [START_REF] Roux | Texture-dependent rigid plastic behavior[END_REF]. Then dilatancy occurs when the volume fraction φ deviates from the equilibrium (i.e. critical) volume fraction φ eq (I) [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF], this is the formulation that we take. The authors of [START_REF] Barker | Partial regularisation of the incompressible µ(I)-rheology for granular flow[END_REF][START_REF] Schaeffer | Constitutive relations for compressible granular flow in the inertial regime[END_REF] propose relations to be satisfied by the yield function (related to the frictional behavior) and the flow rule (related to compressible effects) to insure wellposedness.

The relevance of the model should be tested against laboratory experiments. Quantitative comparison between laboratory experiments and 2D granular flow simulations is however challenging in particular because most experimental flows are 3D, within a channel, and therefore significantly affected by lateral walls [START_REF] Taberlet | Super stable granular heap in thin channel[END_REF][START_REF] Jop | Crucial role of sidewalls in dense granular flows: consequences for the 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]. Furthermore, introducing wall effects in shallow depth-averaged models may lead to large errors because these effects depend on the flowing depth and not on the whole depth as assumed in depth-averaged models [START_REF] Fernández-Nieto | Narbona-Reina, 2D granular flows with the mu(I) rheology and side walls friction: a well balanced multilayer discretization[END_REF]. As a result, we will mainly focus here on the qualitative effect of dilatancy on the flow and deposit even though the results will be compared to laboratory experiments.

We present in Section 2 the compressible (further called dilatant) µ(I) model which is obtained by simply removing the fluid in the two-phase model derived in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF]. We then propose in Section 3 a numerical method to solve the equations, pointing out the numerical difficulties. We then investigate in Section 4 compression/dilation effects in dry granular column collapse showing the subtle role of dilatancy on flow dynamics, effective friction and final deposit.

Dry granular model with dilatancy

We consider a one-dimensional version of the 2D equations presented in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF], for the case of a dry granular material with dilatancy effects.

To obtain this model we start from the Jackson's model -originally defined for two-phase systems made of fluid and grains-(cf. [START_REF] Jackson | The Dynamics of Fluidized Particles. Cambridges Monographs on Mechanics[END_REF]) where the fluid is replaced by the vacuum. Thus, the conservation of mass and momentum read:

∂ t (ρ s ϕ) + ∇ • (ρ s ϕv) = 0, ρ s ϕ(∂ t v + (v • ∇)v) = -∇ • T + ρ s ϕg.
The velocity field is denoted by v, T denotes the (symmetric) stress tensor and the constant solid density is denoted by ρ s . The solid volume fraction is ϕ and g is the gravity. In this system neither buoyancy nor friction with air are accounted for. We consider the dilatancy law

∇ • v = γ tan ψ (2.1)
as the closure for the system above (see also equation (2.8) below).

The derivation of the model follows the classical depth-averaging process under the shallow flow asymptotic regime (see Appendix A for details). We obtain a closed system with the unknowns ϕ, h and v that are averaged quantities representing the volume fraction, the thickness of the granular layer in the direction normal to a fixed inclined plane, and its velocity in the slope-parallel direction, respectively (see Figure 1). Notice that equation (2.1) is used to obtain the mass equation in the depth-averaging procedure, resulting in the source term of the solid volume fraction continuity equation (2.2). Moreover this equation is used to obtain an expression of the normal velocity v z , even if it does not appear in the equations. Thus, the original dilatancy law (2.1) is recovered at the correct order of accuracy with respect to the asymptotic approximation of the derived system (see Remark A.1 for more details).

The resulting model, with x defined as the slope aligned coordinate, reads

∂ t ϕ + v∂ x ϕ = -ϕ γ tan ψ, (2.2)
∂ t (ϕh) + ∂ x (ϕhv) = 0, (2.3)
ρ s ϕ(∂ t v + v∂ x v) = -ρ s ϕg c ∂ x ( b + h) -ρ s g c h 2 ∂ x ϕ - τ b h sgn(v), (2.4)
with g c = g cos θ, where g is acceleration due to gravity and θ the slope of the inclined plane. We set b(x) = b(x) + x tan θ with b denoting the bottom elevation in the direction normal to the inclined plane (see Figure 1). Note that equations (2.2) and (2.3) lead to the equation on the thickness h

∂ t h + ∂ x (hv) = h γ tan ψ. (2.5)
The term sgn(v) in (2.4) represents the solid friction, and is multivalued in the sense that its value can be any number in [-1, 1] when v = 0. The model satisfies the momentum equation ((4.17) in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF])

∂ t (ρ s ϕhv) + ∂ x (ρ s ϕhv 2 ) + ∂ x ρ s ϕg c h 2 2 = -ρ s ϕg c h∂ x b -τ b sgn(v), (2.6)
and the energy equation ((4.41) in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF])

∂ t ρ s ϕh v 2 2 + ρ s ϕg c h 2 2 + ρ s ϕg c h b + ∂ x ρ s ϕh v 2 2 v + ρ s ϕg c h 2 v + ρ s ϕg c h bv ≤ ρ s ϕg c h 2 2 γ tan ψ -τ b |v|, (2.7)
the inequality being strict only when discontinuities (shocks) appear. We use the dilatancy law tan ψ = K(ϕϕ eq ).

(2.8)

The quantity γ is the (averaged) shear rate and ϕ eq is the critical-state compacity given by

ϕ eq = ϕ stat -K I, (2.9)
where ϕ stat is a constant volume fraction corresponding to the static equilibrium, K, K are constants, and I is the inertial number. In our case of dry material it is defined as

I = γ d p s| b /ρ s , (2.10)
where d is the mean diameter of the solid particles, ρ s is the density of the solid material, and p s| b is the bottom solid pressure,

p s| b = ρ s ϕg c h. (2.11)
The velocity equation (2.4) involves the bottom solid pressure through the effective bottom solid friction

τ b = µ eff p s| b , (2.12)
where we define

µ eff = max(0, tan δ eff ) (2.13) with tan δ eff = µ(I) + tan ψ and µ(I) = µ 1 + µ 2 -µ 1 I + I 0 I. (2.14)
The constant values µ 1 , µ 2 are the tangents of the internal friction angles and I 0 is a constant parameter, all of them involved in the µ(I) rheology. The cutoff at 0 in the definition of µ eff avoids negative values, even if this possibility does not occur in practice. We compute the shear rate with the approximation

γ = 5 2 |v| h , (2.15)
compatible with the free fall or inertial regime (cf. [START_REF] Cassar | Submarine granular flows down inclined planes[END_REF]). In the formulation (2.7) of the energy equation, we do not include the dilation internal energy e eq c terms of [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] (equation (C7) in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF]) because, with our above defined dilatancy law, this would not give a nonpositive right-hand side. This is because ϕ eq does not depend only on p s| b , but also on γ.

The system has the steady state at rest (i.e. with vanishing velocity) solutions characterized by

v = 0, ∂ x ( b + h) + h 2ϕ ∂ x ϕ ≤ µ eff . (2.16)
Note that the volume fraction equation (2.2) and the definition (2.8) of γ tan ψ indicate that ϕ eq is an attractive value for ϕ. To obtain a relevant volume fraction ϕ ∈ [0, 1], we therefore need to have ϕ eq ∈ [0, 1]. This is not the case with the definition (2.9) of ϕ eq since I can take arbitrary large values, in particular I ∼ h -3 2 when h approaches zero. Therefore we consider several modifications of (2.9) that remain in [0, 1] while keeping the same behavior when I is small. We consider either (f0) ϕ eq = max(0, ϕ stat -K I), or the following approximations that allow a smooth decay (f1) ϕ eq = ϕ stat e

-K ϕ stat I , (f2) ϕ eq = ϕ stat 1 + K N ϕstat I N , for a given N > 0.
Figure 2 shows the effect of these approximations for dense initial packing at inclination angle 22 • and t = 1.38s. As the inertial number tends to infinity near the front (Figure 2(b)), the equilibrium volume fraction ϕ eq reaches unphysical negative values owing to equation (2.9), as shown by the black circles in Figure 2(a). The simplest regularization (f0) is very sharp and generates oscillation of the mass thickness near the front at θ = 0 • (Figure 2(c)).

Regularization (f1) is smoother but still gives rise to small instabilities near the front at θ = 0 • where the thickness h is small. The decreasing rate of change of ϕ eq towards zero can be controlled by the exponent N in regularization (f2). As N increases regularization (f2) tends to (f1) (N ≥ 5 is almost superimposed to (f1)). For N = 1 the decrease to zero is smoother than (f1) but still close to ϕ eq and no more oscillations occur (red line in Figure 2(c)). At higher slope angles (θ = 22 • ), all regularizations give a stable h (Figure 2(d)). This may be due to greater thicknesses at the front at higher slope angles whereas at small angles the front thickness decreases very slowly to zero. However, whatever the regularization, the system still has a singularity when h is small, because of (2.15) that implies that γ tan ψ in (2.8) tends to infinity like 1/h. This property induces some numerical difficulties close to the vacuum (i.e. close to h = 0). ϕ eq (f0) ϕ eq (f1) ϕ eq (f2), N=0.1 ϕ eq (f2), N=1 ϕ eq (f2), N=100 ϕ eq (a) ϕ eq at θ = 22 • and t = 1.38s and associated regularizations

Numerical scheme

The system (2.2), (2.3), (2.4) can be written as the following shallow water system with volume fraction and dilatancy, In practice, a value of α = 2 is used. Note that although classically the variable considered to write a two-phase system is (1ϕ)h (that would be related here to the vacuum part) instead of ϕ α h, we do not consider it because of the lack of convexity. Indeed it corresponds to the case α = 0 in (3.3), that combined with (3.1) gives

∂ t (ϕh) + ∂ x (ϕhv) = 0, (3.1)
∂ t (ϕhv) + ∂ x ϕhv 2 + 1 2 g c ϕh 2 = -g c ϕh∂ x b - τ b ρ s sgn(v), (3.2)
∂ t (ϕ α h) + ∂ x (ϕ α hv) = -(α -1)ϕ α h γ tan ψ, (3.3
∂ t ((1 -ϕ)h) + ∂ x ((1 -ϕ)hv) = h γ tan ψ.
The system (3.1), (3.2), (3.3) is discretized with a combination of finite volume and finite difference schemes, in a two-step splitting approach. The finite volume method is devoted to the system in which the dilatancy source term involving γ tan ψ is excluded, while the finite difference scheme handles only that source term.

Finite volume step

Denoting U = (ϕh, ϕhv, ϕ α h), the finite volume step to deal with the PDE part of (3.1)-(3.3) (i.e. with γ tan ψ = 0) is performed by the well-balanced formula

U n+1 i -U n i + ∆t ∆x F l (U n i , U n i+1 , ∆z i+1/2 ) -F r (U n i-1 , U n i , ∆z i-1/2 ) = 0, (3.5)
where the upper index n refers to time, the lower index i refers to the space location, and the left/right numerical fluxes F l (U l , U r , ∆z), F r (U l , U r , ∆z) are defined according to Appendix B. The "total topography" z is defined formally by ∂ x z = ∂ x b + µ eff sgn(v) to enter the framework of the apparent topography approach (see [4, section 4.12.1]). More explicitly, this means that at the interface i + 1/2, denoting

U l = U n i , U r = U n i+1 , ∆ b = bi+1 -bi , we define ∆z i+1/2 ≡ ∆z by ∆z = ∆ b + ∆b app , (3.6)
∆b app = ∆x Proj

µ interf eff X interf , v interf g c ∆t , (3.7)
where

v interf = ϕ l h l v l + ϕ r h r v r ϕ l h l + ϕ r h r , (3.8)
X interf is an approximation of

X = -∂ x (h + b) + h 2ϕ ∂ x ϕ , (3.9)
and where

Proj µ (X, Y ) = proj µ proj µ (X) + 2 1 + max(1, -X • Y /µ|Y |) Y , (3.10) with proj µ (X) =    X if |X| ≤ µ, µ X |X| if |X| > µ. (3.11)
The argument X in (3.7) ensures the well-balanced property for the steady states at rest (i.e. to exactly preserve discrete steady states at rest). If we want to approximate the stopping phase well when ϕ is not constant, we have to take into account the term in ∂ x ϕ from (3.9) in the computation of X. We take

X interf = 1 ∆x h l -h r -∆ b + 1 2(α -1) h l + h r ϕ α-1 l + ϕ α-1 r ϕ α-1 l -ϕ α-1 r . (3.12)
Considering that (∂ x ϕ)/ϕ = (∂ x ϕ α-1 )/((α -1)ϕ α-1 ), we can see that this definition is consistent with (3.9). Moreover because of the well-balanced property of the numerical flux (see Appendix B, Theorem B.1), the scheme preserves exactly the discrete data satisfying We thus obtain the well-balanced property, not for all the steady states at rest (2.16), but at least for those that satisfy either ϕ = cst or hϕ 1-α = cst. For ϕ = cst, these steady states at rest are characterized by the property |∂ x ( b + h)| ≤ µ eff . The steady states at rest verifying hϕ 1-α = cst are characterized by the property |∂ x ( b + α-1/2 α-1 h)| ≤ µ eff . In the latter case the free surface is not flat and the volume fraction ϕ is proportional to h 1/(α-1) . Setting α = 2, we obtain that the volume fraction ϕ is proportional to h.

v l = v r = 0, |X interf | ≤ µ interf
Note that a difficulty is encountered when applying the scheme (3.5) to a state with h = 0, because as seen in the formulas of Appendix B, the value of ϕ is really needed in this case, even if not determined by the value of U (which vanishes). For this case, we therefore set

ϕ i = max(ϕ i-1 , ϕ i+1 ) whenever h i = 0. (3.14)
A remark is that in this work the issue of numerical stability is important because of the stiff coupling with the dilatancy source terms. A scheme satisfying an entropy inequality is guaranteed to be stable because the integral of the entropy over the domain decreases, giving a bound on the computed solution and preventing from blow up. Here with such a scheme for the finite volume step (and with the stabilization step defined in Subsection 3.3) our computation is stable. It is indeed important to be sure that if there is an instability this is not due to the numerical flux that we use, but it is rather inherent to the nature of the system that we are solving. We would not be sure of that with another solver without entropy inequality. The steady states hϕ 1-α = cst in (3.13) have no clear physical interpretation, contrary to those for which ϕ = cst. Thus having a scheme that leaves them invariant should not be so important. Nevertheless the use of these steady states hϕ 1-α = cst is a way to build systematically the numerical method in a coherent way (formulas (B. [START_REF] Farin | Insights into the generated seismic signal and dynamics of granular flows on horizontal and sloping beds[END_REF])-(B.32)). These formulas could also be useful in other modeling situations where the system is of gas dynamics type, similar to (B.1)-(B.3).

Finite difference step

To complete the numerical procedure, we need to solve only the reaction term that involves the dilatancy effects, i.e.

∂ t (ϕh) = 0, ∂ t (ϕhv) = 0, ∂ t ϕ = -ϕ γ tan ψ. (3.15)
According to (2.8), γ tan ψ = γ K(ϕϕ eq ). We use a semi-implicit method here since the shear rate is involved in the expression of the dilation term through the equilibrium solid volume fraction ϕ eq . Let us denote the solution after the finite volume step by ϕ * , h * , v * and the final solution by ϕ, h, v. Then γ * = 5

2 |v * | h * and we write ϕ = ϕ * -∆tϕ * γ tan ψ = ϕ * 1 -∆t γ * K(ϕ -ϕ * eq
) . Thus we obtain ϕϕ * eq = ϕ *ϕ * eq -∆t γ * Kϕ * (ϕϕ * eq ), and we deduce the equality

(1 + ∆t γ * Kϕ * )(ϕ -ϕ * eq ) = ϕ * -ϕ * eq . (3.16)
We can deduce the solution for ϕ either directly from (3.16) by setting

ϕ = ϕ * eq + ϕ * -ϕ * eq 1 + ∆t γ * Kϕ * , (3.17)
or using the following approximation that mimics the analytical solution of the corresponding differential equation,

ϕ = ϕ * eq + (ϕ * -ϕ * eq ) e -∆t γ * Kϕ * . (3.18)
Next, we get the components h and v from

ϕh = ϕ * h * , v = v * . (3.19)
The value of ϕ * eq is computed using (f2) with N = 1 instead of (2.9), so that ϕ * eq decays slowly to 0 as I → ∞.

Stabilization step

As we are interested in simulating the spreading of a dry granular column, the scheme should be able to deal with the transition between h > 0 and h = 0 (the so-called dry/wet interface in numerical methods). This is a key and generally delicate issue in the design of numerical schemes solving shallow water type equations. Note that with the splitting numerical method described above, we obtain instabilities when using perfectly dry states in the simulation of the granular column collapse. Indeed when h approaches zero, the blow up concerns the shear rate γ and the inertial number I, that are involved in the computation of various other quantities. A classical approach is therefore to use a threshold on the thickness h and to consider a very thin layer of granular material instead of a perfectly dry state. We can choose a thickness threshold of about tens of particles, because for a thickness under this threshold the numerical results are physically meaningless.

In spite of the above procedure, the scheme still shows instabilities in some dry front configurations. These instabilities are actually due to the handling of the non-standard dilatancy source term that is coupled with the standard topography source term. Therefore, we stabilize the numerical method using a very classical approach that consists in including artificial viscosity. Hence the system we solve actually reads

       ∂ t (ϕh) + ∂ x (ϕhv) -∂ x [ι 1 ∂ x (ϕh)] = 0, ∂ t (ϕhv) + ∂ x (ϕhv 2 + 1 2 g c ϕh 2 ) + g c ϕh∂ x z -∂ x [ι 2 ∂ x (ϕhv)] = 0, ∂ t (ϕ α h) + ∂ x (ϕ α hv) -∂ x [ι 3 ∂ x (ϕ α h)] = -(α -1)ϕ α h γ tan ψ, (3.20
) where we denote formally ∂ x z = ∂ x b + µ eff sgn(v) as above in the finite volume step. The numerical viscosity coefficients ι i , for i = 1, 2, 3, are defined by

ι i = max(0, ξ i v∆z), (3.21)
with ξ i nonnegative constant parameters. For the numerical results, we take ξ 1 = ξ 2 = 2, ξ 3 = 5. These numerical viscosities are of the order of ∆x, thus system (3.20) remains close to (3.1)-(3.3). The choice (3.21) is made because we noticed that instabilities appear only when the quantity v∂ x z is positive. We apply a standard centered finite difference scheme to deal with viscous terms. In addition to the classical CFL condition in the finite volume framework, we require the parabolic CFL condition ∆t ∆x 2 max i=1,2,3

ι i ≤ 1 2 . (3.22)
Figure 3 shows that the viscosity has almost no impact on the calculated thickness profiles whatever the inclination angle. When zooming close to the front, we however see that using such small viscosity makes it possible to eliminate the small instabilities observed at the front (insets in Figures 3(a) and 3(b)). The front position changes more at θ = 0 • than at θ = 22 • where it changes by 0.43% of the total extent. Even though the front position difference is up to 9.1% of the total extent at θ = 0 • , it occurs for very small thicknesses (i.e. h ≤ 5×10 -4 m) that are not physically relevant for a continuum model representing granular materials with a grain diameter of about 5 × 10 -4 m. As a result, the small viscosities used in our numerical method are shown to have negligible effect on the calculated fields.

Application to granular collapse 4.1 Experimental setup

We simulate the experiments of [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF][START_REF] Farin | Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles: Insights from laboratory experiments[END_REF] as done in [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF]. The setup consists of a w = 20 cm wide inclined channel with a rectangular cross-section and plexiglasswalls. A rectangular mass of glass beads with height h 0 = 14 cm and length l 0 = 20 cm (i.e. aspect ratio h 0 /l 0 = 0.7), initially at rest, is released from a reservoir at time t = 0 s by lifting a gate (Figure 4). The glass beads are subspherical, cohesionless, and highly rigid with a diameter d = 0.7 ± 0.1 mm. They flow down a plane inclined at θ, roughened with a glued layer of the same glass beads. The particle density is ρ s = 2500 kg m -3 and the initial volume fraction is estimated at ϕ 0 = 0.62 giving an apparent density of ρ = 1550 kg m -3 . The thickness profiles of the flowing granular mass are recorded by a fast camera.

Note that most of the parameters involved in the µ(I)-rheology (2.14) (cf. [START_REF] Jop | A constitutive law for dense granular flows[END_REF]) and the dilatancy law (2.8) (cf. [START_REF] Roux | Texture-dependent rigid plastic behavior[END_REF]) come from empirical laws for granular material. The parameters involved in the µ(I)-rheology have been largely used but less literature can be found about the dilatancy law. The values used here have been taken from recognized works that we detail next. We use the same rheological parameters µ 1 , µ 2 and I 0 as [START_REF] Ionescu | Viscoplastic modeling of granular column collapse with pressure-dependent rheology[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], summarized in table 1. We also test how the classical approach to accounting for wall effects in shallow depth-averaged models (using a grain/wall friction µ w ) impacts our results. This will be presented in section 4.5. The parameters related to dilatancy effects are K = 4.09, K = 0.2, and ϕ stat = 0.585. The value of K is taken from [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF] even though they were for mixtures of grains and fluid while we are dealing with dry granular flows. The parameters involved

h 0 l 0 x z θ v front Figure 4: Experiment setup
in the definition of the solid volume fraction at equilibrium, K and ϕ stat are taken from [START_REF] Andreotti | Granular media: between fluid and solid[END_REF] (equation 6.7) obtained experimentally for dry granular flows. We performed simulations with two initial volume fractions: ϕ 0-dense = 0.62 as in the experiments, referred to as dense packing, and ϕ 0-loose = 0.55, hereafter referred to as loose packing.

µ 1 µ 2 µ w I 0 d(m) K K ϕ stat ϕ 0-dense ϕ 0-loose 0.

Intrinsic limitation of the model for granular collapse simulation

The current one-dimensional model does not take into account the lifting of the gate. The gate has been shown to impact the flow dynamics but has little effect on the deposit (Figures 14 and15 of [START_REF] Ionescu | Viscoplastic modeling of granular column collapse with pressure-dependent rheology[END_REF]). The presence of the gate reduces the mass velocity and leads to steeper flow front, at least during the first half of the collapse. Because it does not take gate lifting into account, our model overestimates the front position (distance from the gate to the front) during the flow by about 10%. Furthermore, the effect of the lateral walls of the channel could be significant, in particular at higher slope angles (Figures 3 and4 of [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF]). Wall effects globally straighten the granular mass. Because it does not take wall effects into account, our model underestimates the maximum thickness at the rear of the flow and overestimates the front position. At θ = 22 • , our simulation can overestimate the front position by 15% and underestimate the maximum thickness by 25%. However taking into account wall effects in shallow depth-averaged models is challenging and the classical way of doing that (see equation (4.1)) can lead to unrealistic effects (e.g. see Figure 11 of [START_REF] Fernández-Nieto | Narbona-Reina, 2D granular flows with the mu(I) rheology and side walls friction: a well balanced multilayer discretization[END_REF]).

Another limitation of the model is the thin-layer approximation (i. e. shallow approximation) that leads to the hydrostatic pressure when water is not present as it is the case here. Indeed, this approximation leads to overestimate the velocity in the beginning of the collapse where the granular mass has an aspect ratio close to 1 (i. e. it is not shallow). Such effect is observed for example when comparing the results of the hydrostatic multi-layer model 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 laboratory experiments of granular collapse (see their figure 15).

Finally, the depth-averaged nature of the model makes it impossible to take into account any motion or inhomogeneity (such as static or flowing grains) in the direction normal to the slope. Neglecting these effects by depth-averaging the equations may lead to an overestimation of the front position by up to 15% and an underestimation of the maximum thickness by up to 10% at slope θ = 22 • (e.g. Figure 15 of [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]).

Qualitatively speaking, neglecting all of these effects in our model produce the same consequences in our simulations: an underestimation of the maximum thickness of the flowing mass and an overestimation of the front position. Our simulations may therefore overestimate the front position by up to 30% and underestimate the maximum thickness at the rear of the flow by up to 35% on moderate slopes. As a result, it is difficult to quantitatively match experimental results with our simulations given that the combined effect of the gate, wall friction and processes in the direction normal to the slope are not or are poorly accounted for in our model.

The objectives here are therefore only to (i) estimate how dilatancy changes the flow profiles, velocity, forces and deposit, (ii) identify where and when dilation or compression occurs, and (iii) possibly estimate how these effects may improve the match with laboratory measurements even though precise quantitative agreement is not attainable.

Effect of the initial packing 4.3.1 Mass profiles and velocities

In classical simulations without dilatancy (i.e. K = 0, hereafter called the incompressible case), shallow depth-averaged models overestimate the flow velocity especially at the initial instants and on small slopes as shown in Figures 5 and6. The front position is strongly overestimated for θ = 0 • and by up to 40% for θ = 22 • at the initial instants while the maximum thickness is underestimated by up to 25%. The deposit is however quite well reproduced.

Let us investigate how dilatancy changes these results for both initially dense (ϕ 0-dense = 0.62) and loose (ϕ 0-loose = 0.55) columns. The most significant effect of dilatancy occurs at the rear of the mass, near the back wall, changing significantly the maximum flow thickness (Figures 5 and6). For all times and slopes, the maximum thickness is higher (smaller) for initially dense (loose) packing than without dilatancy (K = 0). These results are in very good agreement with those of discrete element modeling of dense and loose granular columns (see Figure 9 of [START_REF] Kermani | Simulation of Collapse of Granular Columns Using the Discrete Element Method[END_REF]), supporting the dilatancy law used in our continuum model. At θ = 0 • , 10 • , the rear part of the mass is better reproduced with initially dense packing, corresponding to the experimental case where ϕ 0 = 0.62. The maximum mass thickness of the deposit is however slightly overestimated for initially dense packing at θ = 22 • (Figure 6(d)). At θ = 0 • , there is almost no difference in the mass profiles near the front between the incompressible, dense and the loose cases. At θ = 10 • and up to 0.48 s at θ = 22 • , for both initially loose and dense packing, the front travels further and its thickness is slightly larger than in the incompressible case. As time increases at θ = 22 • , the incompressible mass profiles start to fall between the dense and loose cases (Figures 6(c At the beginning of the collapse, the mass velocity is similar with or without dilatancy (Figures 7(a part has already stopped and the initially dense mass in the frontal zone travels much faster than the loose one (Figure 8(c)). Later on, at t = 1.8 s when the experimental mass has already stopped, the dense and loose cases have higher velocities than the incompressible case (Figure 8(d)). The complex relative behavior of the initially dense, initially loose and incompressible mass velocity can also be observed for example at t = 0.36 s for θ = 0 • where in some positions, near the rear of the flow, the incompressible case is in between the loose and the dense cases while in other positions, near the front, the incompressible velocity is smaller than in the initially loose and dense cases (Figure 7(b)).

For all slopes, the front position of the initially dense mass is always slightly larger than in the loose case (Figure 9(a)). This is in agreement with DEM simulations of granular flows over curved topography [START_REF] Cagnoli | Grain size and flow volume effects on granular flow mobility in numerical simulations: 3-D discrete element modeling of flows of angular rock fragments[END_REF]. The simulated mass stops before the experimental mass at small inclinations θ = 0 • , 10 • . Actually, for these inclinations, simulations stopped at t = 0.48 s and t = 0.66 s, respectively while experiments stopped at t = 0.66 s and t = 1.02 s. At larger slopes (θ = 16 • , 22 • ) however, although almost all the mass is at rest at the deposit time of the experiment, there is still a very slow motion of the front in simulations (see velocity profiles for θ = 22 • in Figure 8(d)). The full stopping of the mass takes a very long time (> 7.5 s). Namely, it is achieved at t = 40 s, even if the velocity is smaller than 10 -3 m s -1 from t = 19 s. The deposit profile is however almost the same as the thickness profile at t = 1.8 s, as shown in the inset in Figure 6(d).

The front velocity is known to be strongly overestimated during the initial instants in shallow depth-averaged models (Figure 9(b)) in agreement with Figure 9 of [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF]. For θ < 22 • , there is almost no difference between the front velocity of the loose and dense cases, even though the front velocity in the dense case is slightly higher than in the loose case. However at θ = 22 • , the change of curvature of the front velocity as a function of time is very well reproduced by the initially dense simulation (corresponding to the experimental configuration) and not by the loose one. This change of regime from a decelerating flow to an almost steady flow at slow velocity is a crucial and poorly understood point in granular flow behavior, as discussed in [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF][START_REF] Farin | Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles: Insights from laboratory experiments[END_REF][START_REF] Delannay | Granular and particleladen flows: from laboratory experiments to field observations[END_REF][START_REF] Farin | Insights into the generated seismic signal and dynamics of granular flows on horizontal and sloping beds[END_REF]].

Dilatancy and friction

The different behaviors of the initially dense and loose cases can be understood by looking at the volume fraction, dilatancy angle, friction coefficient and forces involved (Figures 10 and12). The initial difference in volume fraction ϕ in the dense and loose cases tends to disappear near the front and over most of the mass as time increases (Figures 10(a), 10(b), 12(a) and 12(b)). In these regions, the mass behaves in the same way in the two cases in terms of dilatancy, as if the memory of the initial state had been mostly forgotten. Only the rear part of the flow maintains the initially dense and loose packing fractions. In this region, in the dense (loose) case, ψ > 0 (ψ < 0) so that the mass dilates (contracts) (Figures 10(c), 10(d), 12(c) and 12(d)) and the friction coefficient is large (small) at all times (Figures 10(e), 10(f), 12(e) and 12(f)). The volume fraction near the back wall stays however constant because at these instants the velocity is zero (Figure 7). At θ = 0 • , the front dilates during the initial instants (t = 0.06 s) and then contracts (t = 0.36 s) in both initially dense and loose cases (Figures 10(c angles are |ψ| 8 • initially but they can reach much higher values during the flow, for example near the front at t = 0.36 s where ψ 16 • (Figure 10(d)). Smaller dilatancy angles are observed at intermediate times for θ = 0 • and at all times for θ = 22 • (Figures 12(c) and 12(d)) with values of ψ smaller than 1 • . Indeed, at higher slopes, the volume fraction ϕ becomes rapidly close to the equilibrium volume fraction ϕ eq . Actually, as the slope increases, the system gets closer to the slopes where steady uniform flows are obtained. As observed on the horizontal slope, at θ = 22 • the front dilates at the initial instants (Figure 13), then the mass experiences dilation at the front, contraction just behind and dilation at the rear part (Figure 13(c) and 13(d)). Finally almost the whole mass contracts before stopping (Figure 13(f)).

The combined effect of µ(I) rheology and dilatancy leads to a variable effective friction coefficient (equation (2.13)). At θ = 0 • the increase of friction at the initial instant (t = 0.06 s) near the front results from the µ(I) effect, while at t = 0.36 s its decrease to much smaller values at the front is related to contraction (Figures 10(e) and 10(f), respectively). At θ = 22 • , the µ(I) effect dominates due to small values of dilatancy, leading to an effective friction coefficient that increases towards the front (Figure 12(e) and 12(f)). Interestingly, without dilatancy, µ(I) always increases towards the front, i.e. towards smaller thicknesses in our setup, while dilatancy can completely change this behavior leading to decreasing friction at the front for flows on small slopes. Increasing the slope angle leads, in the dense case, to increased mean dilation of the mass (Figure 14). Even though dilatancy angles are smaller at θ = 22 • than on a horizontal plane, the dilation of the mass lasts longer so that the overall dilation is higher reaching almost 10%, in good agreement with experimental results. The mass dilation is overestimated in the model but the general trend is quite well reproduced for all slope angles. However, in the model, the contraction of the mass after the dilation is underestimated. This could be partly due to the depth-averaging process that does not account for the difference of behavior in the direction perpendicular to the slope. Interestingly, the results obtained in the initially loose case (compression) are very different from the experiments that correspond to the dense case.

As a result, taking into account dilatancy and the initial dense volume fraction measured in the experiments improves the match between simulated and experimental mass thickness profiles and mean dilation/compression as a function of time.

Forces

At the beginning of the flow, all the forces play a significant role and dilatancy has little affect on them as shown in Figures 15(a) and 15(c). The largest effect of dilatancy is observed on the friction force and on the pressure gradient forces while acceleration, inertia and gravity forces are less sensitive to dilatancy. Later on, the flow at θ = 0 • mostly results from a balance between pressure gradient and friction forces, even though the acceleration force is still present (Figure 15(b)). At θ = 22 • , the flow is controlled essentially by the equilibrium between gravity and friction forces while the pressure gradient force is still present (Figure

15(d))

. In that case, the acceleration and inertial forces balance each other out. At these later times, dilatancy mostly affect the forces at the rear part of the mass as observed above for example on thickness profiles (Figure 15(b) and 15(d)). Note that at θ = 22 • acceleration and inertia forces are also slightly affected by dilatancy which may contribute to explain the sensitivity of the front velocity to dilatancy at this angle (see Figure 9(b)).

f acc (∂ t (ϕhv)) f inert (∂ x (ϕhv 2 )) f press (-∂ x (g c ϕh 2 /2)) f grav (-g c ϕh∂ x b) f fric (-τ b sgn(v))
dense loose

Effect of the dilation coefficient

The value of the dilation coefficient K involved in the computation of the dilatancy angle (tan ψ = K(ϕ-ϕ eq )) was determined empirically for fluid-grain mixtures and may not be the best value for dry granular flows. Figure 16(a) shows that increasing K increases the effect of dilatancy on thickness profiles up to K 4, but hardly changes the thickness profiles for K ≥ 4. Increasing K changes the velocity profiles in a more complex way (Figure 16(b)). Indeed in the dense case, at t = 0.36 s and θ = 0 • , the maximum velocity first increases from K = 0 to K = 1 and then decreases for K = 4 and K = 8. At the rear part of the mass, the behavior is the opposite.

Increasing K decreases (increases) the volume fraction ϕ near the rear part of the mass in the dense (loose) case (Figure 17(a)) and increases (decreases) significantly the dilatancy angles even between K = 4 and K = 8 (Figure 17(b)). This leads to an increase of the effective friction at the rear of the mass and a decrease of this coefficient near the front in the dense case (Figure 17(c)). As a result, stronger friction forces are observed in the dense case near the rear of the mass (Figure 17(d)). In this region the thickness gradients are higher for higher K, leading also to higher pressure gradient forces. Note that the increase of friction compensates the increase of the driving force due to surface gradients. The overall dilation of the mass increases when K increases but the qualitative behavior is very similar (Figure 18(a)). Similar behavior is observed at the rear of the mass at higher slopes such as θ = 22 • (Figure 19(a)). The front position and shape are almost insensitive to the value of K if K > 1. The velocity profiles are more affected by the value of K (especially for K > 4) in the loose case than in the dense case at θ = 22 • (Figure 19(b)). This can also be observed on the front position (Figure 20(a)) and front velocity (Figure 20(b)). Interestingly, the change in curvature of the front velocity and the development of a low velocity regime observed in the experiments is reproduced by the dense simulation and is only slightly sensitive to the value of K for K > 1. On the contrary, this behavior is observed in the loose case for K = 1 becomes very different from the experiments for K > 4 (dashed red and blue lines in Figure 20(b)). At these instants at θ = 22 • , the front seem to have a different behavior than the rest of the mass behind, with its own dynamics (see forces in Figure 21(a)). Figure 21(b) shows that in the dense case for all K and in the loose case for K = 1, where we observed a change in curvature of the front velocity similar to the experiments, the front dilates whereas in other cases, the front compresses. The ability to reach a slow, almost steady regime, could As observed on horizontal planes, the dilation of the whole mass increases for increasing K but remains stable for values of K > 4 (Figure 18(b)). Concerning the dilation of the whole mass, values of K smaller than 4 seem to better reproduce the results for all the slopes (Figure 18).

Mimicking lateral walls

The model does not take into account wall effects that significantly change the flow dynamics of granular column collapses as shown by [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF]. In the 2D case (horizontal-vertical), an approximation of wall effects can be calculated using the in-plane flow assumption from the 3-D equations [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF]. For depth-averaged models, wall friction is usually taken into account by adding a friction term in the friction coefficient [START_REF] Taberlet | Super stable granular heap in thin channel[END_REF][START_REF] Jop | Crucial role of sidewalls in dense granular flows: consequences for the rheology[END_REF][START_REF] Ionescu | Viscoplastic modeling of granular column collapse with pressure-dependent rheology[END_REF]] In (a) we show forces computed with dilation coefficient K = 4.09, the legend is the same as in Figure 15: orange, f acc (∂ t (ϕhv)); cyan, f inert (∂ x (ϕhv 2 )); green, f press (-∂ x (g c ϕh 2 /2)); purple, f grav (-g c ϕh∂ x b); strawberry, f fric (-τ b sgn(v)), full (dotted) lines are related to the initially dense (loose) case.

µ eff,w = µ eff + µ w h w , (4.1)
where w is the width between the lateral walls and µ w is the Coulomb friction coefficient between the grains and lateral walls. As mentioned in the introduction, this approximation may lead to large errors. In fact, the friction term depends on the factor µ w p(z)sign(v) (cf. [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF]), so it is bigger at the bottom and smaller at the surface due to the pressure p(z). The approximation in (4.1) ignores this z-dependence since it is written in terms of the total height (see the discussion and numerical tests in [START_REF] Fernández-Nieto | Narbona-Reina, 2D granular flows with the mu(I) rheology and side walls friction: a well balanced multilayer discretization[END_REF]). The adjustment of the parameter µ w may be used to artificially compensate for this error. A first general observation is that, contrary to dilatancy effects, lateral wall friction affects the whole mass and not only the rear part, even at θ = 0 • (Figure 22). These effects are however stronger near the rear part of the mass as h is higher there. As observed in [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF], the wall effects straighten the mass, increasing the thickness at the rear and decreasing the runout distance, thus leading to overall steeper mean slope of the deposit. The match with experimental results is better when adding wall effects at θ = 0 • but at θ = 22 • , the mass is much too steep with wall effect. This is due to the depth-averaged (i.e. monolayer) model that cannot properly account for wall effects as shown in [START_REF] Fernández-Nieto | Narbona-Reina, 2D granular flows with the mu(I) rheology and side walls friction: a well balanced multilayer discretization[END_REF] (their Figure 11). Indeed, using a multilayer model [START_REF] Fernández-Nieto | Narbona-Reina, 2D granular flows with the mu(I) rheology and side walls friction: a well balanced multilayer discretization[END_REF] or a model solving the full Navier-Stokes equations [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF] strongly improves the wall effect description as wall friction applied to the flow only involves the flowing thickness and not the whole thickness. At θ = 0 • , wall effects are of the same order of magnitude as dilatancy effects. Wall effects do not change the qualitative behavior of the flow as shown by the volume fraction, dilatancy angle or friction coefficient (Figure 23). These effects decrease the overall dilatancy of the mass by a few % (Figure 24). As for dilatancy, wall effects strongly change the front velocity during the low velocity regime after the mass deceleration (Figure 25). As a result, the front velocity in this regime provides an excellent constraint for the models.

Conclusion

We have proposed a shallow depth-averaged model for dry granular flows that takes into account volume fraction variations and dilatancy. This model was deduced from the twophase model proposed in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] by removing the fluid phase. The dilatancy law is derived from the works of [START_REF] Roux | Texture-dependent rigid plastic behavior[END_REF] and [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF]. The equations are solved numerically by a well-balanced finite volume method that uses a numerical flux resolving exactly two families of steady states at rest and satisfying a semi-discrete energy inequality. Friction is dealt with by the apparent topography method. The dilatancy law is solved by a simple explicit finite difference step. A small numerical viscosity is applied to stabilize (i) the singularity when the thickness approaches zero, which occurs at the front, and (ii) the coupling between convection of volume fraction and dilatancy. This small viscosity is shown to have a negligible effect on the calculated fields.

We performed a series of simulations of granular column collapse experiments on planes with inclination varying from horizontal to θ = 22 • [START_REF] Farin | Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles: Insights from laboratory experiments[END_REF]. By simulating the collapse of initially dense and initially loose columns, we investigated the effect of dilatancy on the flow dynamics and deposits. We compared these results with incompressible simulations where dilatancy was not taken into account. Overall, the simulations obtained for the dense case show better agreement with experimental results that indeed involve initially dense columns. We showed that the strongest effect of dilatancy in dry granular column collapses is observed at the rear part of the flow. The maximum thickness of the initially dense column is higher than that of the incompressible case which is itself higher than the initially loose column. This is in very good agreement with the discrete element modeling results of [START_REF] Kermani | Simulation of Collapse of Granular Columns Using the Discrete Element Method[END_REF], supporting the dilatancy laws used in our continuum model. At the front, the flow behaves in a similar way for both the initially dense and loose cases, with similar values of the volume fraction. The flow however keeps some memory of its initial packing fraction. At high slopes the front simulated in the initially dense case travels further than in the loose case.

During the flow, the dilation/compression pattern is complex. At the initial instants, the front dilates for all slopes and for both the initially loose and dense cases. Later on, the mass may contract at the front while it dilates behind the front and contracts in the rear part. Before stopping, most of the mass contracts. Contraction at the front at θ = 0 • is strong at the beginning and end of the collapse, with dilatancy angles ψ reaching up to 16 • . Dilatancy angles at θ = 22 • are much smaller (a few degrees) at all times but as the flow duration is longer, the overall dilatancy of the mass is greater as the slope increases. It can reach up to about 10% at θ = 22 • , in good agreement with laboratory measurements [START_REF] Martin | Continuum viscoplastic simulation of a granular column collapse on large slopes: µ(I) rheology and lateral wall effects[END_REF]. For smaller angles, dilatancy is overestimated and the compression of the whole mass in the dense case, following its first dilation, is underestimated. This could be partly due to the depth-averaging process that does not account for the difference of behavior in the direction perpendicular to the slope.

A strong effect of dilatancy is observed on the front velocity at θ = 22 • . At this angle, laboratory experiments show the emergence of a regime of slow, almost steady flow, after the deceleration phase of the mass [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF]. This regime is reproduced by our simulation when using dilatancy for the initially dense case but not for the incompressible case. For the initially loose case, this regime can be obtained for a specific value of the dilation coefficient K. In this regime, very small changes of the friction, gravity or pressure gradient forces, almost in equilibrium, can change the flow acceleration. As a result, the front velocity in this regime is very sensitive to the flow rheology. In particular we showed that the experimental results are well reproduced only for flows in which dilation occurs at the front. This suggests that dilatancy may play a key role in the development of such regimes.

The effective friction depends on both the µ(I) rheology and the dilatancy angle ψ. At the front, the effective friction increases as I increases. This is always true at θ = 22 • where dilatancy angles are small. However, at θ = 0 • , the effect of front contraction may dominate, leading to a friction coefficient decreasing towards the front (e.g. with decreasing flow thicknesses). As a result, taking into account dilatancy significantly changes the value of the friction coefficient in the model.

The coefficient K involved in the dilatancy law that best reproduces experimental results is slightly smaller than the coefficient proposed for grain/fluid mixtures in [START_REF] Pailha | A two-phase flow description of the initiation of underwater granular avalanches[END_REF]. Varying this coefficient only quantitatively change the results while the overall behavior of the mass is qualitatively the same. We found that dilatancy effects are of the same order of magnitude as wall effects at θ = 0 • , even though taking into account wall effects in shallow depth-averaged models [START_REF] Jop | Crucial role of sidewalls in dense granular flows: consequences for the rheology[END_REF][START_REF] Taberlet | Super stable granular heap in thin channel[END_REF] may not be accurate [START_REF] Fernández-Nieto | Narbona-Reina, 2D granular flows with the mu(I) rheology and side walls friction: a well balanced multilayer discretization[END_REF].

For all slopes, the front position of the initially dense mass is always slightly larger than in the loose case as observed in DEM simulations [START_REF] Cagnoli | Grain size and flow volume effects on granular flow mobility in numerical simulations: 3-D discrete element modeling of flows of angular rock fragments[END_REF], contrary to what is observed for immersed granular columns [START_REF] Rondon | Granular collapse in a fluid: role of the initial volume fraction[END_REF][START_REF] Bouchut | A twophase solid-fluid model for dense granular flows including dilatancy effects: comparison with submarine granular collapse experiments[END_REF]. Indeed, the coupling between fluid pore pressure and compression/dilation increases (decreases) significantly the friction coefficient in dilation (compression) for the initially dense (loose) case. Dilatancy effects are much greater for fluid/grain mixtures than for dry granular flows. In the literature some references dealing with dilatancy in biphasic model can be found, see for instance [START_REF] Wang | Dilatancy and compaction effects on the submerged granular column collapse[END_REF][START_REF] Meng | Modeling dynamic flows of grain-fluid mixtures by coupling the mixture theory with a dilatancy law[END_REF][START_REF] Monsorno | Time-accurate calculation of two-phase granular flows exhibiting compaction, dilatancy and nonlinear rheology[END_REF][START_REF] Lee | A two-phase flow model for submarine granular flows: With an application to collapse of deeply-submerged granular columns[END_REF]. In these works the dilatancy effects are analyzed in particular configurations such as for submerged case or uniform flows. Nevertheless these effects need to be deeply investigated even though solving these equations numerically is very challenging. Our results provide new insight into the dilatancy behavior of dry granular flows that may help better understand dilatancy in grain/fluid mixtures.

The velocity, volume fraction and dilatancy law are then supposed to not depend on z up to errors in O( 2 ), they are written in general as ξ = ξ + O( 2 ).

We calculate the integration of (A.1a) and use the boundary conditions to get:

∂ t (h φ) + ∂ x (h φv) = O( 2 ). (A.8)
From the dilatancy closure (A.2) and the mass equation (A.1a) we obtain the evolution of the solid volume fraction:

∂ t φ + v∂ x φ = -φγ tan ψ + O( 2 ). (A.9)
The pressure is obtained from the vertical momentum equation using the vertical components of the boundary condition (A.6) that leads to the hydrostatic pressure:

p = φρ s g cos θ(b + h -z) + O( 2 ).
Finally we perform the averaging of the tangential component of the momentum equation using the tangential component of the boundary conditions for the bottom (A.5) and the free surface levels (A.6). It reads

ρ s φ(∂ t v + v∂ x v) = -ρ s g cos θ φ∂ x ( b + h) + h 2 ∂ x φ -ρ s g cos θ φ tan δ eff sgn(v) + O( 2 ).
(A.10) The dry granular model studied in the paper is based on equations (A.8)-(A.10).

Remark A.1. Notice that the dilatancy law (A.2) is recovered up to order 3 . Namely, integrating equation (A.2) and using the non-penetration condition at the bottom (A.4), we get the normal velocity:

v z = v∂ x b + (z -b)( γ tan ψ -∂ x v) + O( 3 ).
This expression of the vertical velocity gives us also the asymptotic dilatancy law. Moreover, in the precedent work [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] dealing with the biphasic model, this expression is useful to calculate the excess pore pressure (see section 3.5 of the paper).

Remark A.2. The model presented here can be seen as a particular case of the model in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] when the fluid is replaced by the vacuum which means to set ρ f = 0. The absence of fluid also makes to consider β = 0 since now there is no friction between phases. Indeed, (A.8)-(A.10) is the same model as those in equations (4.10-4.13) in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] choosing ρ f = 0 and β = 0. Note that since β = 0 implies p e fm = 0, this choice can be seen as a formal way to obtain the system presented here. By taking ρ f = 0 and β = 0 in the original system proposed in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF], we have to see that some quantities like u f , 1ϕ, h f are not zero, but when looking at the equations for ϕ, v, h m , all the terms involving the fluid have ρ f or β in factor, thus it remains a closed system in terms of ϕ, v, h m and we do not have to compute h f and u f .

B Numerical fluxes for the shallow water system with volume fraction

In this appendix we define numerical fluxes for the shallow water system with volume fraction

∂ t (ϕh) + ∂ x (ϕhv) = 0, (B.1) ∂ t (ϕhv) + ∂ x ϕhv 2 + 1 2 g c ϕh 2 + g c ϕh∂ x z = 0, (B.2) ∂ t ϕ + v∂ x ϕ = 0, (B.3)
where h ≥ 0 is the thickness, ϕ > 0 the volume fraction, v the velocity, g c > 0 the component of acceleration due to gravity normal to the slope of the inclined plane, and z the (total) topography (satisfying here ∂ t z = 0). The system (B.1), (B.2), (B.3) can be interpreted as the classical full gas dynamics system with γ = 2 (by setting the density to ρ = ϕh) with force -g c ∂ x z, where the roles of energy and entropy have been reversed. Our system is very similar to the Ripa system, for which a numerical method has been proposed in [START_REF] Desveaux | Well-balanced schemes to capture non-explicit steady states: Ripa model[END_REF], differing only by the factor in front of the topography term. The consequence is that the steady states differ in the two systems. A conservative form for (B.3) is

∂ t (hϕ α ) + ∂ x (hϕ α v) = 0, (B.4)
for some α with α = 1 in order to be independent of (B.1). It is important to recall that even for weak solutions with discontinuities, the conservation laws (B.4) (when α varies) are all equivalent, because of the combination with (B.1). This is why we can formally write (B.3).

The system is completed with an entropy (energy) inequality Note that the equation (B.8) for steady states at rest is not well-defined for discontinuous unknowns. Several discontinuous solutions can be obtained depending on the way we understand the nonconservative products.

∂

Hydrostatic reconstruction scheme

Our numerical fluxes are constructed by generalizing the hydrostatic reconstruction scheme [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF][START_REF] Bouchut | A new model for shallow viscoelastic fluids[END_REF][START_REF] Bouchut | A multi well-balanced scheme for the shallow water MHD system with topography[END_REF][START_REF] Bouchut | A subsonic-well-balanced reconstruction scheme for shallow water flows[END_REF][START_REF] Bouchut | A robust well-balanced scheme for multi-layer shallow water equations[END_REF]. Let us define the conservative variable U = ϕh, ϕhv, hϕ α (B.10)

for some fixed α ≥ 1/2, α = 1. We consider a given scheme for the problem without topography F(U l , U r ) = F 0 (U l , U r ), F 1 (U l , U r ), We shall assume that the volume fraction flux is given by the classical upwind passive transport flux F 2 (U l , U r ) = F 0 (U l , U r )ϕ α-1 l if F 0 (U l , U r ) ≥ 0, F 0 (U l , U r )ϕ α-1 r if F 0 (U l , U r ) ≤ 0, (B. [START_REF] Delannay | Granular and particleladen flows: from laboratory experiments to field observations[END_REF] this assumption being satisfied by most solvers, such as for example the Suliciu relaxation solver [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF] that we use in practice. It is well-known that this property (B.13) implies that the volume fraction satisfies the minimum and maximum principles. The entropy and entropy flux of the system are where (x) + = max(0, x), ∆z = z rz l , and λ ≥ 0 is a constant parameter. The classical way of reconstructing is to take λ = 0. The value λ = 1 is intended to somehow take into account the v 2 /2 term in the steady state relation (B.9). Other values of λ appear to be irrelevant. We then have 0 ≤ ∆z * + ≤ (∆z) + , 0 ≤ ∆z * -≤ (-∆z) + , (B.17 We assume now α > 1, (B. [START_REF] Heyman | Compressibility regularizes the µ(I)rheology for dense granular flows[END_REF] and define

η(U ) = ϕh v 2 2 + g c ϕ h 2 2 , G (
h * l =        h l - α -1 α -1/2 ∆z * + if α -1 α -1/2 ∆z * + ≤ h l 1 -

Figure 1 :

 1 Figure 1: Notation and sketch of the domain. The slope angle θ is negative on this plot.

  (f0) dilat. dense (f1) dilat. dense (f2), N = 1 (c) θ = 0 • , t = 0.36 s θ = 22 • , t = 1.38 s

Figure 2 :

 2 Figure 2: Effects of the regularization of ϕ eq . (a) ϕ eq calculated with (2.9) and with regularization functions (f0), (f1) and (f2) for several values of the exponent N , shown near the front position at inclination angle θ = 22 • at time t = 1.38s; (b) Inertial number shown near the front position at inclination angle θ = 22 • . Influence of such regularization on the computed thickness for dense initial packing on a plane inclined at (c) θ = 0 • and (d) θ = 22 • with dilation coefficient K = 4.09 and no lateral wall effect.

eff , and either ϕ l = ϕ r or h l ϕ 1

 1

  visc. dilat. dense no visc. (a) θ = 0 • , t = 0.36 s θ = 22 • , t = 1.38 s

Figure 3 :

 3 Figure 3: Effects of the numerical viscosity on the computed thickness for dense initial packing on a plane inclined at θ = 0 • , 22 • with dilation coefficient K = 4.09 and no lateral wall effect.

Figure 5 :

 5 Figure 5: Thickness profiles computed for initially dense and loose packing on planes inclined at θ = 0 • (left) and θ = 10 • (right) with dilation coefficients K = 0., 4.09 and no lateral wall effect.

Figure 6 :v

 6 Figure 6: Thickness profiles computed for initially dense and loose packing on a plane inclined at θ = 22 • with dilation coefficients K = 0., 4.09 and no lateral wall effect.

Figure 7 :

 7 Figure 7: Velocity profiles computed for initially dense and loose packing on a horizontal plane (θ = 0 • ) with dilation coefficients K = 0., 4.09 and no lateral wall effect.

v

  ) and 8(a)) even though the thickness profile is already affected by dilatancy at these initial times for θ = 22 • (Figure6(a)). At t = 1.38 s for θ = 22 • , the simulated rear (cm/s) dilat. dense (K = 4.09) dilat. loose (K = 4.09) no dilat. dense (K = 0) no dilat. loose (K = 0) t = 1.8 s

Figure 8 :

 8 Figure 8: Velocity profiles computed for initially dense and loose packing on a plane inclined at θ = 22 • with dilation coefficients K = 0., 4.09 and no lateral wall effect.

Figure 9 :

 9 Figure 9: Time series of front positions (a) and front velocities (b) for various plane inclination angles and for initially dense and loose packing with dilation coefficient K = 4.09 and no lateral wall effect.

  ) and 10(d)). The spatio-temporal distribution of dilation/compression is complex (Figures 11).

IFigure 10 :

 10 Figure 10: Time series of (a)-(b) volume fraction, (c)-(d) tangent of the dilation angle, and (e)-(f) effective friction with inertial number and corresponding thicknesses in the inset, for initially dense and loose packing and for the incompressible case, at slope angle θ = 0 • . Simulations are made with dilation coefficients K = 0, 4.09 and no lateral wall effect.

  (K = 4.09) dilat. loose (K = 4.09) no dilat. dense (K = 0) no dilat. loose (K = 0)

Figure 11 :

 11 Figure 11: Zoom of dilatancy angle profiles computed for initially dense and loose packing on the horizontal plane (θ = 0 • ) with dilation coefficients K = 0., 4.09 and no lateral wall effect.

Figure 12 :

 12 Figure 12: Time series of (a)-(b) volume fraction, (c)-(d) tangent of the dilation angle, and (e)-(f) effective friction, for initially dense and loose packing on a plane inclined at θ = 22 • with dilation coefficients K = 0., 4.09 and no lateral wall effect.

21

 21

Figure 13 :

 13 Figure 13: Zoom of dilatancy angle profiles computed for initially dense and loose packing on a plane inclined at θ = 22 • with dilation coefficients K = 0., 4.09 and no lateral wall effect.

Figure 14 :

 14 Figure 14: Time series of relative volume variation for various plane inclination angles and for initially dense and loose packing with the dilation coefficient K = 4.09 and no lateral wall effect.

  θ = 22 • , t = 0.48s

Figure 15 :

 15 Figure 15: Forces computed for initially dense and loose packing on planes inclined at θ = 0 • , 22 • with dilation coefficient K = 4.09 and no lateral wall effect.

Figure 16 :

 16 Figure16: Thickness and velocity profiles computed at t = 0.36 s for various dilation coefficients K = 0., 1., 4.09, 8. and for initially dense and loose packing on a plane inclined at θ = 0 • without lateral wall effect.

Figure 17 :Figure 18 :

 1718 Figure 17: (a) Solid volume fraction, (b) dilatancy angle, (c) effective friction coefficient and (d) relevant forces profiles computed at t = 0.36s for various dilation coefficients K = 0., 1., 4.09, 8. and for initially dense and loose packing on a plane inclined at θ = 0 • without lateral wall effect.

Figure 19 :

 19 Figure19: Thickness and velocity profiles computed at t = 0.48 s for various dilation coefficients K = 0., 1., 4.09, 8. and for initially dense and loose packing on a plane inclined at θ = 22 • without lateral wall effect.

Figure 20 :

 20 Figure 20: Time series of front positions and velocities for various dilation coefficients K = 0., 4.09, 8. and for initially dense and loose packing on a plane inclined at θ = 22 • without lateral wall effect.

Figure 21 :

 21 Figure21: Zoom near the front of forces and dilatancy angle computed for initially dense and loose packing on a plane inclined at θ = 22 • at time t = 1.2 s and no lateral wall effect. In (a) we show forces computed with dilation coefficient K = 4.09, the legend is the same as in Figure15: orange, f acc (∂ t (ϕhv)); cyan, f inert (∂ x (ϕhv 2 )); green, f press (-∂ x (g c ϕh 2 /2)); purple, f grav (-g c ϕh∂ x b); strawberry, f fric (-τ b sgn(v)), full (dotted) lines are related to the initially dense (loose) case.

  θ = 0 • , t = 0.48s (deposit) θ = 22 • , t = 0.48s θ = 22 • , t = 1.8s

Figure 22 :

 22 Figure 22: Thickness profiles with lateral wall friction effect and initially dense packing on planes inclined at θ = 0 • , 22 • with dilation coefficients K = 0., 4.09.

  with dilatancy, µ(I) friction and dense initial packing -nx = 5000) (d) effective friction

Figure 23 :Figure 24 :

 2324 Figure23: Profiles computed at t = 0.36s with lateral wall friction effect and initially dense packing on a horizontal plane (inclination θ = 0 • ) with dilation coefficients K = 0., 4.09.

Figure 25 :

 25 Figure 25: Time series of front velocity with lateral wall friction effect and initially dense packing on a plane inclined at θ = 22 • with dilation coefficients K = 0., 4.09.

2 +

 2 [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF] and(B.6). Therefore the inequality (B.6) is a way to encode the minimum and maximum principle on ϕ, by taking respectively ψ(ξ) = (kξ) + and ψ(ξ) = (ξk) + for an arbitrary constant k ≥ 0.The velocity equation in nonconservative form is obtained from (B.1), (B.2),ϕ(∂ t v + v∂ x v) + g c ϕ∂ x (h + z) g c (h + z) = cst, ϕ = cst. (B.9)

F 2 (

 2 U l , U r ) , (B.11) corresponding to the conservative problem with flux F (U ) = ϕhv, ϕhv 2 + 1 2 g c ϕh 2 , hϕ α v . (B.12)

U ) = ϕh v 2 2 +

 2 g c ϕh 2 v. (B.14)The entropy and entropy flux of the system with topography areη(U, z) = η(U ) + g c ϕhz, G(U, z) = G(U ) + g c ϕhzv. (B.15)The hydrostatic reconstruction scheme is defined via reconstructed states. These are defined as follows. We define first ∆z * + = ∆z -

) and moreover 0 ≤

 0 (∆z) + -∆z * + ≤ (∆z) + -∆z --∆z) + -∆z * -≤ (-∆z) + --∆z -

21 )c ϕ l h 2 l 2 -g c ϕ * l (h * l ) 2 2 = 2 =

 21222 Note that in particular,h * l = (h l -∆z * + ) + if ϕ r ≥ ϕ l , h * r = (h r -∆z * -) + if ϕ r ≤ ϕ l . ∆z ≥ 0 and ϕ r ≤ ϕ l and α-1 α-1/2 ∆z * + ≥ h l 1 -I) the expansion (B.57) is trivial since the left and right-hand sides vanish. In domain (II) we have ϕ * l = ϕ l , h * l = (h l -∆z * + ) + . According to (B.17) we have ∆z * + = O(∆z), thus (assuming h > 0) for sufficiently small U l -U and ∆z, we haveh * l = h l -∆z * + . It follows that g g c ϕ l (h + o(1))∆z * + = g c ϕh∆z * + + o(∆z). (B.58) Next, since U * r = U r , U * l = U l + O(∆z), we have F 0 (U * l , U * r ) = ϕhv + o(1). If v = 0 then F 0 (U * l , U * r )/v l = ϕh + o(1)and it follows thatF 0 (U * l , U * r )g c (∆z) + -∆z * + v l = g c ϕh (∆z) + -∆z * + + o(∆z). (B.59)Adding this to (B.58), we get the desired second line of (B.57). In the case v = 0, (B.18) shows that (∆z)+ -∆z * + = O(v l (v lv))and thusF 0 (U * l , U * r )g c (∆z) + -∆z * + v l = o(v lv). (B.60)Adding this to (B.58) we get the desired second line of (B.57), except that we have an additional o(v lv). This is indeed the weaker form of consistency when the numerical flux is not quite smooth.In domain (III) we haveh * l = h lg c (ϕ lϕ * l ) (h * l ) 2 2 + g c ϕ l h 2 l -(h * l ) 2 2 = g c ϕh∆z * + + o(∆z). (B.63)We therefore conclude in exactly the same way as in domain (II). Finally, in domain (IV) we have h * l = h l + the domain (IV) imply that ϕ rϕ l = O(∆z * + ), and therefore we can write ϕ r ϕ l α-1 = 1 + (α -1) ϕ r ϕ l -1 + o(∆z * + ). (B.65)

  The condition α ≥ 1/2 is assumed because it characterizes the convexity of the energy ϕhv 2 /2 + g c ϕh 2 /2 with respect to the conservative variable (ϕh, ϕhv, ϕ α h).

				)
	for some fixed value α with
	α = 1	and	α ≥ 1/2.	(3.4)
	Note that equation (2.2) has no mathematical meaning when there are discontinuities, the
	last equation (3.3) makes it possible to replace (2.2) by a formally equivalent conservative
	equation that can handle discontinuities. It is important to recall that even for weak solutions
	with discontinuities, the conservation laws (3.3) (when α varies) are all equivalent, because
	of the combination with (3.1). This gives mathematical meaning to the nonconservative
	equation (2.2).

Table 1 :

 1

	48 0.73 0.18 0.279 7.10 -4 4.09 0.2 0.582	0.62	0.55

Experimental rheological parameters. µ 1 , µ 2 , I 0 and d are used in the friction law (equation (2.13)). µ w (grain/wall friction) is also involved when wall friction is considered (equation (4.1)). K, K, and ϕ stat are used in the dilation/compression law (equations (2.8) and (2.9)). Two initial packing fractions are considered corresponding to dense (ϕ 0-dense ) and loose (ϕ 0-loose ) packing.

  It is possible to check that this energy is convex with respect to the conservative variables (ϕh, ϕhv, hϕ α ) if and only if α ≥ 1/2 (and still α = 1). Other entropy inequalities are∂ t ϕhψ(ϕ α-1 ) + ∂ x ϕhψ(ϕ α-1 )v ≤ 0, for ψ convex, (B.6) indeed ϕhψ(ϕ α-1 )is convex with respect to (ϕh, ϕhv, hϕ α ) if and only if ψ is convex. Note however that (B.6) is slightly misleading since, as we said above, weak solutions to (B.1), (B.4) indeed satisfy (B.6) with equality (there is no dissipation), and for all ψ even if not convex. The interest of the formulation (B.6) with inequality is only for numerical schemes, because at the discrete level a scheme can never satisfy the two conservative equations (B.

	t ϕh	v 2 2	+ g c ϕhz + g c ϕ	h 2 2	+ ∂ x ϕh	v 2 2	v + g c ϕh(h + z)v ≤ 0.	(B.5)

Acknowledgements

This research has been partially supported by the Spanish Government under grants MTM 2015-70490-C2-2-R and RTI2018-096064-B-C22 with the participation of FEDER, by the ANR contract ANR-11-BS01-0016 LANDQUAKES, the USPC PEGES project, and the ERC contract ERC-CG-2013-PE10-617472 SLIDEQUAKES.

Appendixes A Derivation of the one-phase model

In this appendix we detail the derivation of the presented model (2.2)- (2.4). The idea of this model is to represent a particular case of the biphasic model proposed in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] when the fluid is replaced by the vacuum. Then, the governing system is given by the 2D Jackson's equations for the solid phase, that reads

where neither buoyancy nor friction with air are accounted for. We consider the dilatancy law to close the system given by (2.1) and (2.8):

The stress tensor is assumed to have the following structure

for p the pressure and T to be defined according to the rheological assumptions.

We complete the system with appropriate boundary conditions considering the one-layer domain in Figure 1.

At the bottom we assume the non penetration condition v • n = 0 at the bottom, (A. [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF] where n is the upward space unit normal to the topography. A solid Coulomb friction law is also considered,

where δ eff = tan δ + tan ψ is the effective intergranular Coulomb friction angle, sgn(v) = v/|v|, and the subscript τ denotes the tangential projection, u τ = u -(u • n)n for any vector u.

At the free surface we assume no tension T N X = 0 at the free surface, (A. [START_REF] Bouchut | A twophase solid-fluid model for dense granular flows including dilatancy effects: comparison with submarine granular collapse experiments[END_REF] together with the kinematic condition

where N = (N t , N X ) is a time-space normal to the free surface.

The development of the averaging process is analogous to that made in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF], so we just mention here the main steps. Indeed the derivation remains similar to the classical Shallow-Water equations for a variable density (here, this role is played by the solid volume fraction). The notation of the domain is given in Figure 1 where the granular layer is flowing over a fixed plane inclined at the angle θ, x being the direction along the slope, and z the direction normal to the slope. The equation of the bottom is given by z = b(x) and the free surface by z = b(x) + h(t, x). The velocity is then decomposed into its longitudinal and normal components, v = (v, v z ).

We assume the same asymptotic regime than in [START_REF] Bouchut | A two-phase twolayer model for fluidized granular flows with dilatancy effects[END_REF] in terms of the aspect ratio between the characteristic width and the length of the domain , supposed to be small under the shallow domain assumption. That is

Then we have

We also define the reconstructed volume fractions

, ϕ r if ∆z ≥ 0 and ϕ r ≤ ϕ l , (B.26)

Finally the reconstructed states are

The numerical fluxes for the nonconservative system (B.1), (B.2), (B.4) are defined by

(B.32) Note that because of (B.18), the divisions by v l or v r in (B.31), (B.32) are not singular, indeed the ratios tend to zero as v l or v r tends to zero. (viii) The volume fraction ϕ satisfies the minimum principle (but not the maximum principle), which means that for any k ≥ 0, the inequality ϕ ≥ k remains true if it holds initially.

Proof. For (i), conservativity follows from the fact that F 0 l = F 0 r and F 2 l = F 2 r . Then if ∆z = 0, according to (B.17

For the well-balanced property (ii), consider U l , U r a steady state at rest (B.8) with constant ϕ, which means that

Moreover, (B.28) means that ϕ * l = ϕ * r = ϕ l = ϕ r . Thus in any case U * l = U * r , and taking into account that v = 0,

For the well-balanced property (iii), we first observe that the relation hϕ

Consider now data U l , U r corresponding to a steady state at rest with hϕ 1-α = cst. Taking into account the previous computation, this means that

Let us assume that ∆z ≥ 0, the converse case being similar. Then h * r = h r , ϕ * r = ϕ r , ∆z * + = ∆z. From (B.37) we have h l ≥ h r , ϕ r ≤ ϕ l , and h l (1 -(ϕ r /ϕ l ) α-1 ) = h lh r . We are thus exactly in the case of equality in the condition involved in (B.20), and this formula gives h * l = h r . Then (B.26) gives ϕ * l = ϕ r . We deduce that U * l = U r = U * r , and we conclude with (B.35) as in the proof of the well-balanced property (ii).

The proof of nonnegativity (iv) is the same as in the standard shallow water system, using only (B.25) and (B.29).

For the semi-discrete energy inequality (v), the proof is more involved and is indeed the origin of the intricate formulas (B.20), (B.21). We assume that the homogeneous solver satisfies a semi-discrete energy inequality, which can be written

for any U l , U r , where G(U l , U r ) stands for the numerical entropy flux. Evaluating these inequalities on U * l , U * r thus gives

For the numerical entropy flux (for the problem with topography), we are going to take the formula

Since the proof of the right semi-discrete inequality is similar, we are going to prove only the left inequality. Comparing the left semi-discrete energy inequality that we need to satisfy to the second line of (B.39), we deduce that it is sufficient to satisfy the inequality

We can compute the derivative of the energy η with respect to the conservative variable U by

Then we compute

Also with (B.42),

Thus the left inequality to be satisfied (B.41) can be written

Now according to (B.31), F(U * l , U * r ) -F l has only the velocity component as nonzero. The velocity component of η (U l ) that it is multiplied by in the last line of (B.45) is then just v l . We then observe that the pressure difference in (B.31), multiplied by v l , will then just cancel the first line of (B.45). In the second line of (B.45), the component F 2 (U * l , U * r ) can be expressed with (B.13),

where the l/r notation means that it is left in case F 0 (U * l , U * r ) ≥ 0, and right in the case F 0 (U * l , U * r ) ≤ 0. Using z *z l = (∆z) + , the inequality to be satisfied (B.45) can be written finally by factorizing F 0 (U * l , U * r ) as

(B.47) When ∆z ≤ 0, the left factor is identically zero, thus let us assume that ∆z ≥ 0. Then h * r = h r , ϕ * r = ϕ r . Consider first the case when ϕ r ≥ ϕ l . Then ϕ * l = ϕ l , h * l = (h l -∆z * + ) + , and the inequality (B.47) becomes

For the case when

which proves (B.48). Consider now the case when ϕ r ≤ ϕ l (still with ∆z ≥ 0). Then h * l is given by (B.20) and ϕ * l is given by the second line of (B.26). In the inequality (B.47), a term can be rewritten

because of the value of ϕ * l given by the second line of (B.26). Thus the inequality (B.47) can be rewritten

We observe that with the definition of h * l ,

α-1 , the value of h * l is given by the first line of (B.20) while the value of ϕ * l is ϕ l (h * l /h l ) 1/(α-1) , and the left-hand side of (B.51), becomes

α-1 , the value of h * l is given by the second line of (B.20), while ϕ * l = ϕ r . The left-hand side of (B.51) becomes

According to the value of h * l given by the second line of (B.20), the expression between parentheses is either zero if h * l > 0, or nonnegative if h * l = 0. In the latter case, we have F 0 (U * l , U * r ) ≤ 0, thus (B.54) is nonpositive, which concludes the proof of (v). For the consistency (vi), we have to prove two properties. The first is that for any U , F l (U, U, 0) = F r (U, U, 0) = F (U ). This property is obvious given (i) and the consistency of the homogeneous flux F. The second property is that

thus it is sufficient to prove that

(B.66)

We therefore conclude in exactly the same way as in domain (II), and this ends the proof of (vi).

Next we prove property (vii) concerning constant data ϕ. If ϕ l = ϕ r , then by (B.28)

It follows that the update U = (U 0 , U 1 , U 2 ) will also satisfy the relation U 2 = ϕ α-1 l U 0 , which proves the claim. We finally prove the minimum principle (viii) via an interface, which means that for some σ l (U l , U r , ∆z) < 0 < σ r (U l , U r , ∆z) we have

where U k is the set of states U satisfying ϕ ≥ k and k ≥ 0 is a fixed constant.

To prove (B.67), we consider U l , U r ∈ U k . We now denote the components of a state U by U = (U 0 , U 1 , U 2 ). The property U ∈ U k involves only the components U 0 and U 2 of U and can therefore be written U 2 /U 0 ≥ k α-1 . Thus in the property on the right-hand side of (B.67), only the components F 0 l/r and F 2 l/r of the numerical flux are involved. According to (B.31), (B.32), these coincide with the components of F(U * l , U * r ). Thus the property (B.67) can be written where l/r stand for "l" if F(U * l , U * r ) (0) ≥ 0, and "r" in the converse case. According to (B.28), we have (ϕ * l/r ) α-1 ≥ k α-1 and since σ l < 0, we deduce that the first inequality in (B.69) holds whenever F(U * l , U * r ) (0) ≤ 0. In the case F(U * l , U * r ) (0) ≥ 0 and we can use (B.29), giving ϕ * l ≤ ϕ l , thus

With the identity (B.70) we conclude that

≥ ϕ α-1 l

. Remark: there is no discrete entropy inequality associated with (B.6) for arbitrary nonincreasing convex ψ .