Ioan R Ionescu
email: ioan.r.ionescu@gmail.com

Anne Mangeney
email: mangeney@ipgp.fr

François Bouchut
email: francois.bouchut@u-pem.fr

Olivier Roche
email: o.roche@opgc.univ-bpclermont.fr

Viscoplastic modelling of granular column collapse with pressure dependent rheology

Keywords: Granular collapse, Viscoplastic rheology, Drucker-Prager plasticity, yield stress, variable viscosity, augmented Lagrangian

We propose a mechanical and numerical model of dry granular flows that quantitatively reproduces laboratory experiments of granular column collapse over inclined planes, with rheological parameters directly derived from the experiments. We reformulate the so-called µ(I) rheology in the framework of Drucker-Prager plasticity with yield stress and a viscosity η( D , p) depending on both the pressure p and the norm of the strain rate tensor D . The granular domain, the velocities, the stress deviator and the pressure fields are calculated by using a finite element method based on an iterative decompositioncoordination formulation coupled with the augmented Lagrangian method.

2-D simulations well reproduce the dynamics and deposit of the collapse of granular columns. The flow is essentially located in a surface layer behind the front, while the whole granular material is flowing near the front where basal sliding occurs. The computed runout and the slope of the deposit are in a very good agreement with the experiments.

Interpretation of the results shows that the dynamic viscosity varies from very small values to 1.5 Pa.s in the context of experimental granular collapse. The respective roles of the pressure and the strain rate in the viscosity (i.e. in the µ(I) rheology) are quantified. Essentially, the viscosity is smaller in the shallow flowing region that is located near the free surface and near the front. Using the order of magnitude of the mean viscosity during the flow (η = 1 Pa.s here) that can be obtained by simple calculation, we show that a Drucker-Prager rheology with a constant viscosity gives results very similar to the µ(I) rheology and agrees as well with experimental height profiles, while significantly reducing the computational cost. Within the range of viscosities 0.1 < η < 1 Pa.s, the dynamics and deposits are very similar, so that 1 the observed slumping behavior appears to be mainly due to the flow/noflow criterion and to the associated strain-independent part of the "flowing constitutive relation" (i. e. related to plastic effects). However, to take an unrealistic large value of the viscosity (10 Pa.s) changes strongly the results. The good results obtained with a constant viscosity have important implication for application to real geophysical flows, which could therefore possibly be modeled using a limited number of parameters.

The way the gate that maintains the initial granular column is removed has a significant impact on the flow dynamics but the deposits obtained with and without the gate are almost the same. Using adherence (zero-slip velocity) instead of friction on the base the simulations better fits the experiments near the left wall but underestimates the runout distance.

Simulations provide insight into the flow dynamics showing heterogeneities within the flow with low and high pressure zones. The yield limit has a gravitational distribution while the stress deviator reaches its maximum value at mid flow height. The regions with a high strain rate, located behind the front tip of the granular mass, near the bottom surface, are not related to high stress deviator zones.

Introduction

The mechanical behavior of dense dry granular flows is of paramount importance in very different domains such as geophysics, physics or industry. An increasing amount of theoretical and experimental studies have been devoted to this question, in particular to infer mechanical properties of geophysical granular flows from field observation. Geophysical flows (rock-falls, rock or debris avalanches, etc.) have been mostly simulated using thin layer depth-averaged models to reduce the high computational cost related to the necessary description of the real topography. Beside the prohibitive computational cost that would be required by using the equations of mass and momentum conservation without depth-averaging and/or thin layer approximation, the lack of well established constitutive relation for these complex natural materials prevented developing such models for application in real situations. Furthermore, the rheological parameters associated with complex natural materials are generally very hard to measure and the properties of the material may significantly change during the flow, in particular due to fragmentation, segregation or entrainment processes. In this context, a key point for real application is to use constitutive relations with the minimum number of parameters, that would describe the mean behavior of the flow. Despite all these difficulties, developing 3D models of granular flows is crucial because some key processes such as the flow/no flow transition present in particular in erosion and deposition mechanisms are hardly handled in thin-layer depth-averaged models. Furthermore 3D structures such as convection cells may play a role in the flows dynamics [START_REF] Brodu | Shallow granular flows down flat frictional channels: Steady flows and longitudinal vortices[END_REF], especially for natural landslides that are intrinsically transient and significantly affected by the topography fluctuations. We send the reader to the reviews [START_REF] Bingham | Fluidity and plasticity[END_REF][START_REF] Gdr Midi Group | On dense granular flows[END_REF][START_REF] Ancey | Plasticity and geophysical flows: a review[END_REF][START_REF] Andreotti | Les milieux granulaires, entre fluide et solide[END_REF] for a general introduction on granular flows.

Even though granular flows at the laboratory scale may not involve the same physical processes than those acting at the natural scale, they provide a very useful way to investigate and quantify possible mechanisms and scaling laws as well as to test constitutive relations. In particular, the transient flows obtained by the release of granular columns at the laboratory scale have been largely used in the recent years to investigate granular flow dynamics and deposits. We will simulate here the collapse of granular columns over horizontal and inclined channels performed experimentally by Mangeney et al. [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF]. These experiments are in the line of a series of experiments of granular collapse, mainly performed on horizontal planes in channels (e.g. dam-break) or over unconfined beds (Lajeunesse et al. [START_REF] Lajeunesse | Spreading of a granular mass on a horizontal plane[END_REF][START_REF] Lajeunesse | Spreading of a granular mass on a horizontal plane[END_REF], Lube et al. [START_REF] Lube | Axisymmetric collapse of granular columns[END_REF][START_REF] Lube | Static and flowing regions in granular collapses down channels[END_REF][START_REF] Lube | Granular column collapses down rough, inclined channels[END_REF], Balmforth and Kerswell [START_REF] Balmforth | Granular collapse in two dimensions[END_REF], Siavoshi and Kudrolli [START_REF] Siavoshi | Failure of a granular step[END_REF], Lacaze et al. [START_REF] Lacaze | Planar collapse of a granular column: experiments and discrete-element simulations[END_REF]). Some of them extended to granular collapse on rigid inclined beds (Hogg et al. [START_REF] Hogg | Two dimensional granular slumps down slopes[END_REF], Mangeney et al. [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF], Lube et al. [START_REF] Lube | Granular column collapses down rough, inclined channels[END_REF], Farin et al. [START_REF] Farin | Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles: Insights from laboratory experiments[END_REF]) and on erodible granular beds (Crosta et al. [START_REF] Crosta | Numerical modeling of 2-D granular step collapse on erodible and non-erodible surface[END_REF], Mangeney et al. [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF], Farin et al. [START_REF] Farin | Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles: Insights from laboratory experiments[END_REF]). Efforts to explain the scaling laws obtained in these experiments have either concentrated on thin layer modeling (Mangeney-Castelnau et al. [START_REF] Mangeney-Castelnau | On the use of Saint-Venant equations to simulate the spreading of a granular mass[END_REF], Kerswell [START_REF] Kerswell | Dam break with coulomb friction: a model for granular slumping[END_REF], Larieu et al. [START_REF] Larrieu | Raining into shallow water as a description of the collapse of a column of grains[END_REF], Doyle et al. [START_REF] Doyle | Static and flowing regions in granular collapses down channels: Insights from a sedimenting shallow water model[END_REF]), Discrete Element Methods (DEM) (Staron and Hinch [57], Zenit [START_REF] Zenit | Computer simulations of the collapse of a granular column[END_REF], Lacaze et al. [START_REF] Lacaze | Planar collapse of a granular column: experiments and discrete-element simulations[END_REF], Girolami et al. [START_REF] Girolami | A three-dimensional discrete-grain model for the simulation of dam-break rectangular collapses: Comparison between numerical results and experiments[END_REF]) and more recently continuum viscoplastic models (Crosta et al. [START_REF] Crosta | Numerical modeling of 2-D granular step collapse on erodible and non-erodible surface[END_REF], Lacaze and Kerswell [START_REF] Lacaze | Axisymmetric granular collapse: a transient 3D flow test of viscoplasticity[END_REF], Meruane et al. [START_REF] Meruane | On the role of the ambient fluid on gravitational granular flow dynamics[END_REF] and Lagrée et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF]).

A lot of numerical studies were focused on the reproduction of the experimental scaling laws for granular collapse over horizontal beds, but quantitative comparison between the simulations and the experiments during the flow and for the deposits are still poorly investigated. Indeed, while all studies reproduced the scaling laws, they often did not match the runout distance and duration of the flow. While thin layer depth-averaged models based on a Coulomb-type friction law are able to well reproduce quantitatively the final deposit for aspect ratios lower than 1 using an empirical friction coefficient slightly higher than the tangent of the friction angle of the involved particles (µ = 0.6), the spreading velocity during the first instants of the collapse is significantly overestimated (Mangeney-Castelnau et al. [START_REF] Mangeney-Castelnau | On the use of Saint-Venant equations to simulate the spreading of a granular mass[END_REF], Kerwell [START_REF] Kerswell | Dam break with coulomb friction: a model for granular slumping[END_REF], Ionescu [START_REF] Ionescu | Viscoplastic shallow flow equations with topography[END_REF]). On the other hand, 2D DEM (contact dynamics) performed by Staron and Hinch [START_REF] Staron | Study of the collapse of granular columns using 2D discrete-grains simulation[END_REF], strongly overestimates the maximum extent of the deposit (up to 40%) of dam-break granular collapse even with a very high grain/grain friction coefficient µ m = 1. Using 3D DEM (soft particle dynamics), Lacaze et al. [START_REF] Larrieu | Raining into shallow water as a description of the collapse of a column of grains[END_REF] reproduced quantitatively dam-break granular collapse of small lateral width (1 or 2 particles) using a grain/grain friction coefficient µ m = 0.35 and an empirical friction coefficient along the lateral wall, which surprisingly varies significantly depending on the size of the granular particles (µ w = 0.15 or 0.35). Similarly, 3D DEM (soft particle dynamics) of laterally wider dam-break granular collapse (Girolami [19]) reproduced quantitatively the dynamics and the deposit observed experimentally with a grain/grain friction coefficient µ m = 0.5, a particlewall friction µ w = 0.5, and an empirical additional dissipation said to be related to the resistance to rolling of the grains along the wall. The runout of unconfined granular collapse was quantitatively reproduced by Lacaze and Kerswell [START_REF] Lacaze | Axisymmetric granular collapse: a transient 3D flow test of viscoplasticity[END_REF] using 3D DEM (soft particle dynamics) with a grain/grain friction coefficient µ m = 0.5. These results suggest that DEM simulations have to be performed in 3D to reproduce quantitatively the experimental collapse of granular columns, even for dam-break configuration.

We propose here alternatively to use a continuum approach to simulate dambreak granular collapse based on viscoplastic constitutive relations in the line of recent studies, and to perform quantitative comparison with laboratory experiments. Crosta et al. [START_REF] Crosta | Numerical modeling of 2-D granular step collapse on erodible and non-erodible surface[END_REF] first simulated horizontal granular collapse using an elastoplastic constitutive relation with a Mohr-Coulomb yield rule involving a constant friction coefficient, nonassociate flow rule, and a small cohesion to avoid numerical prob-lems. They found good agreement with the experimental scaling laws on the deposit, with internal friction angles within a few degrees of the characteristic repose and avalanche angles of the material involved, even though the final maximum thickness of the deposit was slightly underestimated. However, they did not provide detailed comparison between the experiments and the simulation during the flow. On the other hand, Lacaze and Kerswell [START_REF] Lacaze | Axisymmetric granular collapse: a transient 3D flow test of viscoplasticity[END_REF], using 3D DEM, showed that the socalled viscoplastic µ(I) rheology seemed to hold during the 3D collapse of granular columns, where I is denominated as the inertial number. Indeed, a large number of studies show that the µ(I) flow law makes it possible to describe a wide range of experimental observations (e.g. GDR MiDi group [START_REF] Gdr Midi Group | On dense granular flows[END_REF], Silbert et al. [START_REF] Silbert | Granular flow down an inclined plane: Bagnold scaling and rheology[END_REF], Jop et al. [START_REF] Jop | Crucial role of sidewalls in dense granular flows: consequences for the rheology[END_REF][START_REF] Jop | A constitutive law for dense granular flows[END_REF]). Later on, Lagrée et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF] compared 2D continuum and DEM (contact dynamics) simulations of granular collapse. In their continuum model, Lagrée et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF] prescribed no-slip boundary conditions at the bottom, neglecting possible basal sliding. They showed that with a given parameter set, the continuum and discrete approaches gave very similar results. As discussed above, the DEM used in Lagrée et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF], and Staron and Hinch [START_REF] Staron | Study of the collapse of granular columns using 2D discrete-grains simulation[END_REF], strongly overestimated the runout extent of the deposit observed experimentally. With the selected parameter set, the continuum approach predicted smaller runout distance than the DEM, despite they revealed similar scaling laws, but still much larger than those obtained in the experiments (by more than 20%). They stated that the µ(I) rheology gave better results than a friction law with a constant friction coefficient µ even though, at small aspect ratio (a=1.42), their results only showed very little differences between the two laws. These differences are similar to those obtained using µ(I) with different rheological parameters (see their Figure 16). These differences are located near the front of the flow, that traveled faster for µ = cst.

For practical application to natural flows, the possible relevance of a simplified rheology (constant friction and/or constant viscosity) for granular flows over inclined slopes is a critical issue because the three parameters involved in the µ(I) rheology are very difficult to calibrate in nature and because the small viscosities involved in this rheology may induce prohibitive computational time.

We will focus here on the detailed quantitative comparison of Drucker-Prager plastic models with constant and variable viscosity with experimental results on the dynamics and deposit of dam-break granular flows. While former similar simulations were restricted to horizontal planes we will also investigate granular flows over inclined topography. First, we will present the 2D viscoplastic continuum model developed here (section 2.1) where the µ(I) rheology is reformulated in a viscoplastic context with Drucker-Prager plasticity (pressure-dependent yield stress) and a variable viscosity. Then, a quantitative comparison between the simulations and laboratory experiments of granular dam-break with small aspect ratios (typical of natural landslides) over horizontal and inclined planes will be presented to investigate the capability of the viscoplastic model to reproduce the observations (section 3). In the next section we discuss the influence of the gate and of the frictional boundary conditions. New insight into the flow properties (velocity profile, basal sliding, yield limit distribution, yielding surface, stagnant zones, evolutions of the horizontal and vertical layers) resulting from the numerical simulations (section 5), shows strong heterogeneities within the flow. Finally, in section 6, the results and their implication in terms of rheology and practical application will be discussed. In the appendix, one can find a short description of the numerical scheme used in the computations.

Modeling granular collapse 2.1 Pressure dependent viscoplastic fluids

In contrast with a newtonian fluid, which cannot sustain a shear stress, for a rigid viscoplastic fluid at rest (i.e. no strain rate), the Cauchy symmetric stress σ, must belong to a set of admissible rigid stresses. This set is usually defined through a continuous scalar function F = F ( σ ′ 0 , p) that describes the flow/no flow condition, i.e. F ( σ ′ 0 , p) ≤ 0 if and only if the fluid is at rest. We have denoted by p = - 1 3 trace(σ) the pressure, by σ ′ = pI + σ the deviatoric stress and all over the paper A = A : A/2 denotes the second invariant of a deviator A. If flowing, then the stress does not belong to the set of admissible rigid stresses and we have to define a "flowing constitutive equation". We shall neglect the second order effects so that our discussion will be limited to a subclass of (incompressible) viscoplastic fluids characterized by trace(D) = 0,

σ ′ = β( D , p)D if D = 0, F ( σ ′ , p) ≤ 0 if D = 0. ( 1
)
where D = D(u) ≡ (∇u + ∇ T u)/2 is the strain rate tensor and u is the velocity field. It should be noted that in contrast with a classical fluid constitutive equation, for a rigid viscoplastic fluid/solid, the function β is not defined and cannot be continuously extended at D = 0. We will describe here a procedure to define expressions of the constitutive functions F ( σ ′ , p) and β( D , p) such that the model (1) is consistent (i.e. the compatibility conditions of Cazacu and Ionescu [START_REF] Cazacu | Compressible rigid visco-plastic fluids[END_REF] ) are satisfied). One technique is the so-called superposition method (see Cazacu and Ionescu [START_REF] Cazacu | Compressible rigid visco-plastic fluids[END_REF]). The main assumption is that the state of stress in the material, σ ′ , can be represented as the sum of a viscous (rate dependent) contribution σ V and a (rate independent, i.e. depending only on D/ D , and not on D ) contribution S related to plastic effects,

σ ′ = σ V + S. (2)
The viscous part of the stress is taken as for a classical viscous fluid,

σ V = 2η( D , p)D, (3)
where η is a viscosity coefficient which may depend on D and p. The viscous contribution σ V has to be continuous in D and vanishes for D = 0, i.e.

η( D , p) D → 0, as D → 0.

We assume that there is flow only if the yield condition, expressed in terms of the plastic stress deviator S and p, is satisfied, i.e. F ( S , p) = 0. Since by [START_REF] Ancey | Plasticity and geophysical flows: a review[END_REF] the strain rate and the plastic deviator are colinear, the strain rate D is given though a "plastic flow rule"

D = D S S , with F ( S , p) ≤ 0, D F ( S , p) = 0. ( 5
)
Note that the above equation is rate independent and for all admissible plastic stress deviator S we must have F ( S , p) ≤ 0. In order to get S as a function of the strain rate D we have to invert the flow rule [START_REF] Bingham | Fluidity and plasticity[END_REF]. To do that we suppose that there exits a pressure dependent yield limit κ = κ(p) such that the flow/no flow condition F ( S , p) = 0 can be represented as S = κ(p), i.e. F ( σ ′ , p) = σ ′ -κ(p). We can invert now the flow rule to get

S = κ(p) D D • (6)
Note that in contrast with the viscous contribution σ V , the above stress-strain rate relation is not continuous in D, and S does not vanish necessarily for D = 0.

   σ ′ = 2η( D , p)D + κ(p) D D if D = 0, σ ′ ≤ κ(p) if D = 0. (7)
This formula clearly reflects the decomposition of the deviatoric stress as a sum of a rate-dependent viscous term and a rate-independent plastic term with yield stress. We can recast this relation in a different form by writing the strain rate D as a function of the stress deviator (see also [START_REF] Ionescu | Onset and dynamic shallow flow of a viscoplastic fluid on a plane slope[END_REF][START_REF] Ionescu | Viscoplastic shallow flow equations with topography[END_REF]):

D = 1 2η( D , p) 1 - κ(p) σ ′ + σ ′ , (8)
where [x] + = max(0, x) is the positive part. The above expression was used by Perzyna [START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF] and Duvaut Lions [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF] in extending inviscid plastic models to account for rate effects (viscoplastic regularization method). For κ(p) ≡ 0 the plastic effects are vanishing and (7) reduces to a viscous fluid model. As an example, if η is independent of D and p, (7) reduces to the incompressible Navier-Stokes model. If the plastic effects are present then different choices of yield limit κ(p) can be considered. For constant κ (i.e. κ(p) ≡ κ 0 pressureindependent plasticity) we recover the Von-Mises plasticity criterion σ ′ ≤ κ 0 , introduced to describe the plasticity of metals. If η is constant (independent of D and p), the constitutive equation [START_REF] Bouchut | A shallow model including staticflowing transition for viscoplastic Drucker-Prager materials[END_REF] recovers the classical Bingham model (see [START_REF] Bingham | Fluidity and plasticity[END_REF]).

In this paper we consider a linear dependence of the yield limit with the pressure κ, called the Drucker-Prager plasticity (flow/no flow) criterion (see [START_REF] Drucker | Soil mechanics nd plastic analysis of limit design[END_REF]):

κ(p) = κ 0 + µ s p, (9)
where κ 0 is the cohesion and µ s = tan(δ s ), with δ s the internal frictional angle. This yield criterion was constructed as a simplification of the Mohr-Coulomb plasticity criterion.

Here we will use two different models, one with a constant viscosity η (i.e. independent of D and p) that we will call the "constant viscosity model" (or the Druker-Prager fluid), and one with a variable viscosity η( D , p), chosen so that we recover the model proposed by Jop, Forterre and Pouliquen [START_REF] Jop | A constitutive law for dense granular flows[END_REF] through the inertial number I. The inertial number I, which is the square root of the Savage number or of the Coulomb number introduced by Savage [START_REF] Savage | The mechanics of rapid granular flows[END_REF] and Ancey et al. [START_REF] Ancey | A theoretical framework for very concentrated granular suspensions in steady simple shear flow[END_REF], respectively, can be interpreted as the ratio between two timescales: the inertial microscopic timescale of particle rearrangement d/ p/ρ s (where d is the grain diameter and ρ s is the grain density) and a macroscopic strain rate time scale (1/ D ), i.e.

I = 2 D d p/ρ s . (10)
Note that we have to assume that p > 0 for this to make sense. Following [START_REF] Jop | A constitutive law for dense granular flows[END_REF] we introduce the variable friction

µ(I) = µ s + µ 2 -µ s 1 + I 0 /I , (11)
where I 0 is a non-dimensional constant and µ 2 ≥ µ s is the limiting value of the friction coefficient for large I. Then the µ(I) rheology proposed in [START_REF] Jop | A constitutive law for dense granular flows[END_REF] can be written trace(D) = 0,

   σ ′ = µ(I)p D D if D = 0, σ ′ ≤ µ(0)p if D = 0. (12)
It can be identified with the constitutive equation ( 7) by considering the Drucker-Prager plasticity [START_REF] Cazacu | Compressible rigid visco-plastic fluids[END_REF] with vanishing cohesion (κ 0 = 0) and

β( D , p) = µ(I)p D , (13)
or equivalently by taking the specific dependence of the viscosity η on D and p given by 2η( D , p) D = (µ(I) -µ s )p, i.e. with [START_REF] Crosta | Numerical modeling of 2-D granular step collapse on erodible and non-erodible surface[END_REF],

η( D , p) = (µ 2 -µ s )p 2 D + I 0 k √ p , (14)
where k = d √ ρ s . Note that this viscosity η, corresponding to the decomposition [START_REF] Bouchut | A shallow model including staticflowing transition for viscoplastic Drucker-Prager materials[END_REF], differs from the "effective viscosity" of [START_REF] Jop | A constitutive law for dense granular flows[END_REF] which is β/2, in relation with the decomposition (1).

Experimental setup

Let us describe briefly the granular collapse experiments of Mangeney et al. [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF].

The experimental setup consists of a narrow channel between plexiglas walls, spaced by W = 10 cm (around 140 particles). The planar channel is 3 m long with possible inclination angles θ varying from horizontal up to 30 • (Figure 1). A rectangular granular mass of thickness h 0 = 14 cm (around 200 particle diameters) and of downslope length r 0 = 20 cm (around 286 particles), i.e. an aspect ratio a = h 0 /r 0 = 0.7, is released from a reservoir at time t = 0 s by opening a gate. The glass beads are subspherical, cohesionless, and highly rigid with a diameter d = 0.7 ± 0.1 mm. They are flowing down an inclined channel, roughened by gluing a layer of the same beads on its surface. The particle density ρ s = 2500 kg.m -3 and volume fraction ν= 0.62 of the mass were estimated, which gives an apparent flow density of ρ = νρ s = 1550 kg.m -3 , that is used in [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF]. The parameter k in the variable viscosity η( D , p) rheology ( 19) is k = 0.035 kg 1/2 .m -1/2 . The length of the deposit r f measured from the front of the initial mass located at x = 0, i.e., the runout distance, and the final thickness of the deposit at the upper wall h f have been systematically recorded as well as the time at which the front stops t f . The profiles of the granular mass as a function of time have been measured using a high-speed camera.

In the experiments, a gate is removed at the initial instant to release the granular mass. The removing of the gate is simulated here by considering a lifting velocity V b = h 0 /t b = 2.3 m.s -1 (here t b = 0.06 s is the lifting time) as measured in the experiments of Mangeney et al. [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF] and Farin et al. [START_REF] Farin | Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles: Insights from laboratory experiments[END_REF]. The barrier is considered in the numerical scheme as any other rigid boundary for which the no penetration condition (u • n = 0) is assumed but the position of the gate changes at each time step. For sake of simplicity, we neglect here the friction between the gate and the granular material, even though friction is expected to occur as shown for example in the Figures 6 and 14 of [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF] (see section 4.1 for a detailed analysis of the gate effect).

Problem statement

To model the granular collapse experiment described above we consider here the evolution equations in the time interval (0, T ), T > 0 describing the in-plane flow of a rigid visco-plastic fluid in a domain D(t) ⊂ R 2 with a smooth boundary ∂D(t).

The notation u stands for the 2-D velocity field (u = (u x , u y ), u z = 0) and σ for the 3-D stress tensor field (σ xz = σ yz = 0), p = -trace(σ)/3 is the pressure and σ ′ = σ + pI the in-plane stress deviator tensor (σ Mass and momentum conservation. The momentum balance law (in the Eulerian coordinates) reads

′ xz = σ ′ yz = σ ′ zz = 0).
ρ ∂u ∂t + (u • ∇)u -div σ ′ + ∇p = ρf in D(t), (15)
where ρ > 0 is the mass density distribution and f denotes the body forces (gravity for our purpose). Since we deal with an incompressible fluid, we have

div u = 0 in D(t). (16)
The viscoplastic fluid domain D(t) is transported with the fluid as

∂1 D(t) ∂t + u • ∇1 D(t) = 0, (17)
where 1 D(t) is the characteristic function of the domain.

Constitutive laws.

In the next we consider the constitutive equation ( 7) with D = D(u) ≡ (∇u + ∇ T u)/2 and with the Drucker-Prager plasticity [START_REF] Cazacu | Compressible rigid visco-plastic fluids[END_REF] with vanishing cohesion (κ 0 = 0). Two different choices of the viscosity η will be analyzed: the constant viscosity model (η = cst)

trace(D) = 0,    σ ′ = 2ηD + µ s p D D if D = 0, σ ′ ≤ µ s p if D = 0, (18)
and the variable viscosity model with η given by ( 14),

trace(D) = 0,      σ ′ = 2 (µ 2 -µ s )p 2 D + I 0 k √ p D + µ s p D D if D = 0, σ ′ ≤ µ s p if D = 0, (19)
with k = d √ ρ s , which is equivalent to the µ(I) rheology [START_REF] Doyle | Static and flowing regions in granular collapses down channels: Insights from a sedimenting shallow water model[END_REF]. Note that the dividing σ ′ and p by ρ s involves the dynamical viscosity and pressure η/ρ s , p/ρ s . Indeed any of the previous rheologies can be written in terms of the dynamical variables σ ′ /ρ s , p/ρ s , η/ρ s . In this way the grain density ρ s disappears, and there remains only the rheological parameters µ s , µ 2 , I 0 /d .

Boundary and initial conditions. In order to complete the above equations with the boundary conditions we assume that ∂D(t) is divided into two disjoint parts ∂D(t) = Γ b (t) ∪ Γ s (t). On the boundary Γ b (t), the fluid is in contact with a rigid body through the Coulomb friction law

u • n = 0,    σ T = -µ C [-σ n ] + u T |u T | if u T = 0, |σ T | ≤ µ C [-σ n ] + if u T = 0, (20)
where µ C is the Coulomb friction coefficient at the boundary, n stands for the outward unit normal on ∂D(t) and we have adopted the following notation for the tangential and normal decomposition of the velocity field u and surface stresses σn

u = u n n + u T , with u n = u • n, σn = σ n n + σ T with σ n = σn • n,
where σ = σ ′ -pI. Note that our boundary conditions at the bottom are different from previous studies of continuum viscoplastic simulations of granular collapse that imposed a no-slip condition (Crosta et al. [START_REF] Crosta | Numerical modeling of 2-D granular step collapse on erodible and non-erodible surface[END_REF] and Lagrée et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF]), corresponding to the limit µ C → ∞.

On the (unknown) free surface boundary Γ s (t), which is computed from the advection equation ( 17), we impose a stress free condition:

σn = 0 on Γ s (t). ( 21
)
The initial conditions are given by

u| t=0 = 0, D(0) = D 0 . ( 22
)
Finally the problem is to find the domain D(t), the velocity field u(t), the pressure p(t), and the deviatoric stress tensor σ ′ (t) satisfying equations ( 15), ( 16), [START_REF] Fortin | Méthodes de Lagrangien augmenté, application à la résolution de problèmes aux limites[END_REF] and [START_REF] Gdr Midi Group | On dense granular flows[END_REF] (or [START_REF] Girolami | A three-dimensional discrete-grain model for the simulation of dam-break rectangular collapses: Comparison between numerical results and experiments[END_REF]) with the boundary and initial conditions (20)- [START_REF] Hogg | Two dimensional granular slumps down slopes[END_REF].

Rheological parameters

The parameters used in the model, which have to be derived from the experiments are: µ s , µ 2 and I 0 for the η( D , p) rheology [START_REF] Girolami | A three-dimensional discrete-grain model for the simulation of dam-break rectangular collapses: Comparison between numerical results and experiments[END_REF], µ s and η for the constant viscosity rheology [START_REF] Gdr Midi Group | On dense granular flows[END_REF]. For the two rheologies, we have to impose the friction at the base µ b and the friction at the upper left plexigalss wall of the reservoir µ w . Note that as the simulations are 2D (downslope/vertical directions), the model does not simulate the friction at the two lateral walls bordering the channel. This effect is only accounted for empirically in the effective friction coefficients as described bellow. In Mangeney et al. [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF], the repose angle θ r = 23.5 o ± 0.5 o (µ r = tan θ r = 0.43 ± 0.01) and the avalanche angle θ a = 25.5 o ± 0.5 o (µ a = tan θ a = 0.48 ± 0.01) of the material have been measured by adding material on top of a pile at a small rate and measuring the angle of the pile with the horizontal after and before an avalanche, respectively.

The glass beads used by [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF] are very similar to those used by Pouliquen and Forterre [START_REF] Pouliquen | Friction law for dense granular flows: applications to the motion of a mass down a rough inclined plane[END_REF] and Jop et al. [START_REF] Jop | Crucial role of sidewalls in dense granular flows: consequences for the rheology[END_REF] and are expected to be characterized by the same rheological parameters. As explained in Pouliquen and Forterre [START_REF] Pouliquen | Friction law for dense granular flows: applications to the motion of a mass down a rough inclined plane[END_REF], the friction coefficients in the µ(I) rheology can be obtained by fitting the curve h stop (θ), where h stop is the maximum thickness of deposit resulting from steady uniform flows. The resulting minimum and maximum friction angles fitted by Pouliquen and Forterre [START_REF] Pouliquen | Friction law for dense granular flows: applications to the motion of a mass down a rough inclined plane[END_REF] and Jop et al. [START_REF] Jop | Crucial role of sidewalls in dense granular flows: consequences for the rheology[END_REF] are about θ s = 21 o (µ s = 0.38) and θ 2 = 32.8 o (µ 2 =0.64), respectively. Slightly higher friction angles are observed in the experiments of [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF] (compare Figure 5 of [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF] and Figure 3 of [START_REF] Pouliquen | Friction law for dense granular flows: applications to the motion of a mass down a rough inclined plane[END_REF]). This may result from the width of the channel that is much smaller (W = 10 cm) in [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF] than in [START_REF] Pouliquen | Friction law for dense granular flows: applications to the motion of a mass down a rough inclined plane[END_REF] (W = 70 cm). Indeed, the effect of the lateral walls can be modeled as an additional friction coefficient with a maximum value of the order of µ w h/W where µ w is the grain/wall friction and h is the thickness of the flowing layer (see equation (1) of Taberlet et al. [START_REF] Taberlet | Super stable granular heap in thin channel[END_REF] or equation (4.5) of Jop et al., [START_REF] Jop | Crucial role of sidewalls in dense granular flows: consequences for the rheology[END_REF]).

By assuming that the beads slip against the lateral side walls and that the induced stress is a pure solid friction with a constant coefficient of friction µ w , Jop et al. [START_REF] Jop | Crucial role of sidewalls in dense granular flows: consequences for the rheology[END_REF] measured µ w = tan(10.5 o ) = 0.18. We have therefore chosen µ w = tan(10.5 o ) = 0.18 for the friction on the plexiglass wall at the upper end left side of the reservoir. Note that this value is also consistent with the small grain/wall friction in the DEM of Lacaze et al. [START_REF] Lacaze | Planar collapse of a granular column: experiments and discrete-element simulations[END_REF] but much smaller than in Girolami et al. [START_REF] Girolami | A three-dimensional discrete-grain model for the simulation of dam-break rectangular collapses: Comparison between numerical results and experiments[END_REF]. In [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF], as the maximum flowing thickness is about 0.05 m, the additional friction contribution related to the walls should be at maximum 0.18×0.05/0.1 = 0.09, while in Pouliquen and Forterre [START_REF] Pouliquen | Friction law for dense granular flows: applications to the motion of a mass down a rough inclined plane[END_REF] this additional friction would be 0.18 × 0.01/0.7 = 0.0026 (see their Figure 10 for an approximation of the flowing depth), so that we expect a difference in the friction coefficients of about 0.1 between Mangeney et al. [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF] and Pouliquen and Forterre [START_REF] Pouliquen | Friction law for dense granular flows: applications to the motion of a mass down a rough inclined plane[END_REF]. As a result, friction coefficients corresponding to the experiments in [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF] The value of I 0 = 0.279 is taken the same as in Pouliquen and Forterre [START_REF] Pouliquen | Friction law for dense granular flows: applications to the motion of a mass down a rough inclined plane[END_REF]. Our values of the η( D , p) parameters (i.e. µ(I)) are summarized in Table 1 and are consistent with those deduced from DEM by Lacaze and Kerswell [START_REF] Lacaze | Axisymmetric granular collapse: a transient 3D flow test of viscoplasticity[END_REF] that reproduced quantitatively the scaling laws for axisymetric glass beads collapse (µ s ≃ 0.4, µ 2 ≃ 0.7 and I ≃ 0.3). Lagrée et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF] used µ s = 0.32, µ 2 = 0.6 and I = 0.4, overestimating the runout of experimental granular collapse.

µ s (δ s ) µ b (δ b ) µ w (δ w ) µ 2 (δ 2 ) I 0 k 0.
Another parameter which has to be chosen is the frictional coefficient µ C involved in the boundary condition [START_REF] Glowinski | Augmented Lagrangian and Operator Splitting method in Non-Linear Mechanics[END_REF], modeling the contact with the surrounding rigid walls. At the base of the flow µ C = µ b is the basal friction coefficient, along the upper left wall of the reservoir µ C = µ w is the wall friction coefficient. As the basal roughness was provided by gluing a single layer of the same particles to the channel base, we chose the same friction coefficient at the basal surface as for the internal friction µ b = µ s = tan(25.5 • ).

To speed up the calculation, we have to impose a minimum value of the viscosity η = 0.1 Pa.s (i.e. a kinematic viscosity ν of 4 to 6 ×10 -5 m 2 .s -1 ) for both rheologies. This is necessary with the augmented Lagrangian method, even if the problem can be well-posed (at least in the Bingham case) without viscosity [START_REF] Bouchut | Convergence of conforming approximations for inviscid incompressible Bingham fluid flows and related problems[END_REF][START_REF] Lusso | Modélisation numérique des écoulements gravitaires viscoplastiques avec transition fluide/solide[END_REF]. The values of this threshold viscosity and of the constant viscosity η are discussed in section 3.

Variable viscosity versus constant viscosity

Let us describe here some numerical settings we have used in the numerical simulations (see the Appendix for a brief description of the numerical scheme). The time step was chosen to be between ∆t = 0.0005 s and ∆t = 0.001 s, while the edge size of the mesh is between h min = 0.0033 m and h max = 0.01 m. Concerning the re-meshing process we have used a variable metric/Delaunay automatic meshing algorithm (see [START_REF] Hecht | New development in freefem++[END_REF]) such that the finest mesh is located in the zones of small pressure and of high velocity (see the blue mesh on Figure 6). We have also checked that the mass loss during the simulations is very small (smaller than 0.8%). The remove of the barrier (gate) was simulated by a moving friction less boundary condition (see section 4.1 for details and the gate influence).

Collapse over an horizontal plane θ = 0 o

Figure 2 shows that the model with variable viscosity η( D , p) well reproduces the shape of the granular mass during its spreading over an horizontal plane, with rheological parameters directly derived from the experiments. The colors, representing the norm of the velocity, show that the flow is concentrated only in a shallow sub-domain near the free surface as observed in the experiments and as obtained in the numerical simulations of Crosta et al. [START_REF] Crosta | Numerical modeling of 2-D granular step collapse on erodible and non-erodible surface[END_REF], Lacaze et al. [START_REF] Lacaze | Planar collapse of a granular column: experiments and discrete-element simulations[END_REF][START_REF] Lacaze | Axisymmetric granular collapse: a transient 3D flow test of viscoplasticity[END_REF] and Lagrée et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF]. The computed runout distance r f and the slope of the deposit are in very good agreement with the experiments while the spreading of the mass is slightly faster during the flow. This was also the case in the discrete element simulation of Lacaze et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF] (their Figure 4) when compared to their experiments. The front is slightly more rounded in our simulation than in the experiments, an effect even more pronounced in the simulation of Crosta et al. [START_REF] Crosta | Numerical modeling of 2-D granular step collapse on erodible and non-erodible surface[END_REF] (see their Figure 1a). The thickness of the mass on the left side is decreasing slightly in the simulations while it stays constant in the experiments. Figure 2 shows that already at t = 0.06 s, the velocity in the top left side of the domain is non zero. As a result, the final maximum thickness of the deposit h f is about 10% smaller than in the experiments. Similar results were obtained in the continuum and DEM simulations of Crosta et al. [START_REF] Crosta | Numerical modeling of 2-D granular step collapse on erodible and non-erodible surface[END_REF] (their Figure 6) and Lacaze et al. [START_REF] Lacaze | Planar collapse of a granular column: experiments and discrete-element simulations[END_REF] (their Figures 4, 6 and7), respectively. In our simulation, at time t = 0.76 s, the maximum velocity norm is lower than 0.03 m/s and at t f = 1.02 s an equilibrium configuration is reached. During this long period of time t ∈ [0.76, 1.02]s, a lot of slow flow events occur, changing very slightly the shape of the granular mass before the granular material reached the final equilibrium state. These processes are also observed in the experiments. Figure 3 shows that the viscosity varies from 0.1 Pa.s (minimum threshold viscosity in the model) at the free surface to about 1.5 Pa.s at the base of the flow (see equation ( 14)). As an example, the viscosity within the flowing layer at t = 0.3 s is 0.1 < η < 0.7 Pa.s. The small viscosity at the free surface is directly related to the small pressure in this zone (Figure 4 right column). While the pressure distribution is quite regular within the flowing region, the strain rate is concentrated near the front and at the base of the flowing layer (Figure 4 left column). Its maximum value is around 70 s -1 , so that the first term in the denominator of equation ( 14) is lower than 140 s -1 . The pressure is already about 200 Pa near the front, which gives a second term in the denominator of ( 14) equal to about 112 s -1 near the front. As a result, in the region of high strain rate as near the front, the viscosity arises from both the strain rate and pressure values. At the other places, the first term is much smaller than the second term and the strain rate can be neglected in the viscosity formula [START_REF] Duarte | Arbitrary LagrangianEulerian method for NavierStokes equations with moving boundaries[END_REF]. Interestingly, while the strain rate is quite localized, the viscosity varies smoothly in the flowing region. The smooth increase of the pressure with depth seems to dominate the viscosity variation. This suggests to neglect the strain rate in [START_REF] Duarte | Arbitrary LagrangianEulerian method for NavierStokes equations with moving boundaries[END_REF], giving

η ≃ (µ 2 -µ s ) k I 0 √ p. ( 23
)
This approximation is equivalent to linearizing the µ(I) law [START_REF] Crosta | Numerical modeling of 2-D granular step collapse on erodible and non-erodible surface[END_REF] for small I,

µ(I) ≃ µ s + (µ 2 -µ s ) I I 0 . ( 24
)
The order of magnitude of the viscosity can be simply obtained from equation ( 14) without performing the simulation, only from the knowledge of the initial thickness of the granular column h 0 , of the grain diameter and density, and of the parameters of the µ(I) rheology. Indeed, if we take a pressure p = ρ s gh 0 /2 = 2500×9.81×0.07 = 1716 Pa, assuming that the flowing layer has a thickness of about the half of the initial thickness, a velocity u = √ gh 0 ≃ 1.17 m.s -1 (see Figure 9 of [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF]) and a strain rate D = u/(h 0 /2) ≃ 17 s -1 , equation ( 14) gives a viscosity η ≃ 1.2 Pa.s. Note that, for this calculation, the strain rate in equation ( 14) is negligible. When looking at I and µ(I), the results are more difficult to interpret in terms of dissipation than by looking at the viscosity (Figure 5). The inertial number I only varies slightly between zero in the static region and less than 0.1 in most of the flowing region, except near the front and close to the free surface where it reaches values of about 0.8. Next to the free surface, I is not well defined because the pressure is equal to zero (see equation [START_REF] Cazacu | Steady-state flow of compressible rigidviscoplastic media[END_REF]). These values of I lead to very high values of µ(I) at the free surface. However as µ is multiplied by p, these high values are compensated by an almost zero pressure near the surface. This is the same in the very vicinity of the front (Figure 5). Within most of the flowing region, the friction coefficient is 0.5< µ(I) <0.54, corresponding to a friction angle between 26.5 o and 28 o . µ(I) increases significantly in the front zone with values higher than 0.6 in the very front and behind. The friction varies essentially from higher values near to the bottom up to smaller friction near the surface, except at the very top where the meaning of µ(I) is questionable as discussed above. This is different from the viscosity that increases regularly when going from the left of the mass to the front within the flowing region.

For the constant viscosity rheology, we tested η = 1 Pa.s, which is the order of magnitude of the 'real' viscosity derived from equation ( 14), a 10 times smaller viscosity η = 0.1 Pa.s which is the minimum viscosity that we can handle with a reasonable computational time, and a 10 times larger viscosity η = 1 Pa.s.

Figure 6 shows that very similar results are obtained when using a constant viscosity η = 1 Pa.s and the variable viscosity η( D , p). With this set of parameters, the mass is flowing slightly less rapidly than with the variable viscosity and the runout is slightly smaller (Figure 6). While at t = 0.18 s there is almost no difference between the simulations with the two rheologies, the final runout is 5% smaller with η = 1 Pa.s. This is because in the flowing layer and near the front the variable viscosity is smaller than η = 1 Pa.s (Figure 3). With the variable viscosity, a thinner front zone seems to develop at the very extremity of the mass (t = 0.3 s and t = 1.02 s in Figure 6). Because the deposit is very thin near the front and several saltating beads are observed in the experiments, it is difficult to discriminate which law provides the best fit to the experimental results at the front (Mangeney et al. [START_REF] Mangeney | Erosion and mobility in granular collapse over sloping beds[END_REF], Farin et al. [START_REF] Farin | Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles: Insights from laboratory experiments[END_REF]). The velocities distribution and the shape of the mass are very similar when using both the constant (η = 1 Pa.s) and the variable viscosity. When using η = 0.1 Pa.s, there is no strong differences at times t < 0.3 s, but at later time, while the shape of the main mass is very similar (see Figure 7), the front is thinner and is running faster. The final runout is about 5% larger than with the variable viscosity, in agreement with Lagrée et al. (2011). Still the rest of the deposit behind the front has a very similar shape with η = 0.1 Pa.s, η = 1 Pa.s, and with the variable viscosity. As a result, within this range of viscosities 0.1 < η < 1 Pa.s, the dissipation due to viscous effects is much smaller than that due to plastic deformation (i. e. strain-independent part of the "flowing constitutive relation"). For a higher though unrealistic viscosity (η = 10 Pa.s), however, viscous effects strongly change the flow dynamics and deposit (see Figure 8). Note that the computation is much faster using the constant viscosity η = 1 Pa.s than using the variable viscosity which itself is faster than using η = 0.1 Pa.s.

Collapse over an inclined plane θ = 16 o

Up to now, simulations of granular collapse using a viscoplastic rheology have been performed only on an horizontal plane. Here we simulate the collapse over a moderate slope θ = 16 o using the same set of parameters. The results are qualitatively similar to those obtained on an horizontal plane. The simulation with the variable viscosity reproduces quite well the mass spreading (Figure 9). The collapse of the mass at the left upper channel end is however too fast in the simulations, and the front is also advancing more rapidly than in the experiments. The shape of the final deposit is very well reproduced, except next to the left wall where the computed maximum thickness is smaller, and in the front that is slightly overestimated in the simulations. One possible explanation is that, for θ = 16 o , the flowing thickness (corresponding to the warmer colors in Figure 9) is thicker than for θ = 0 o (Figure 2). The maximum flowing thickness is about 0.1 m for θ = 16 o while it is about 0.05 m for θ = 0 o . As a result, according to Taberlet et al. [START_REF] Taberlet | Super stable granular heap in thin channel[END_REF] and Jop et al. [START_REF] Jop | Crucial role of sidewalls in dense granular flows: consequences for the rheology[END_REF], the effect of the lateral wall friction should be two times bigger in the experiments at θ = 16 o (see section 2.2 for more details). This is not taken into account in our 2D simulations where the lateral wall friction and the other friction coefficients are taken constant (i.e. not dependent on the flowing depth). Taking into account the real effect of the wall friction at θ = 16 o (3D simulations) would lead to smaller runout distances and to a smaller deposit thickness near the left wall, that would better match the experimental observations (see the next section).

The variable viscosity is slightly smaller near the static/flowing transition than at θ = 0 o but is almost the same within the flowing region where 0.1 < η < 0.7 Pa.s. The viscosity decreases towards the front where it is smaller than 0.5 Pa.s. As an example, at t = 0.56 s, the viscosity is lower than 0.3 Pa.s over the 10 cm behind the front. In the flowing region, the viscosity gradually increases with the distance perpendicular to the free surface, seemingly mainly influenced by the increase of the pressure. Indeed, the viscosity mainly shows patterns similar to the pressure field (see section 5), masking the influence of the strain rate. The strain rate distribution shows interesting features, with high strain rates concentrated initially at the base near the front, and further along most of the flowing region. The strain rate is distributed within almost the whole thickness of the flow near the front and over a significant thickness near the base beyond the front (Figure 11 left column). The constant viscosity model leads to smaller front propagation velocities and therefore to smaller runout distances, in better agreement with experiments. At t = 0.23 s, periodic patterns of the strain rate distribution develop, possibly related to flow instabilities (see section 5). Similar patterns are observed for I and µ(I) (Figure 12).

Figures 13 shows again that the results are very similar for the constant (η = 1 Pa.s) and variable viscosity rheologies. For the constant viscosity the mass spreads slower and the runout distance is smaller, in better agreement with experimental results.

Influence of the boundary conditions

Let us recall that, because of the very similar results obtained with the variable η( D , p) and the constant viscosity (η = 1 Pa.s), we will use a constant viscosity with η = 1 Pa.s in all the following simulations because it required less computational time.

Influence of the gate

The removing of the gate at the initial instant has always been thought to have negligible influence on the granular collapse dynamics and deposit. To investigate this effect, we considered two scenarios: in the first one the mass is released instantly and in the second one the granular mass is in frictionless contact with a gate, removed at the speed V b = 2.3 m.s -1 in the direction perpendicular to the inclined bed.

In the scenario with the gate the flow initiates near the bottom of the granular column (dark blue on Figure 14), while in the scenario with no gate (light blue mesh), the mass spreads from the top of the granular mass. As expected, at the beginning of the flow the scenario with the gate is much closer to the experimental data represented in pink lines.

The initial lifting of the top right edge of the mass in the experiments (t = 0.06 s in Figure 14 and15) is not reproduced in the simulations. This is because in the simulations, there is a frictionless contact between the mass and the gate, while in reality friction occurs at the contact.

During all the spreading, in the scenario with the gate, the front is advancing less rapidly and the thickness of the mass is generally smaller near the front than when the gate is not taken into account. At later time (e. g. t = 0.3 s in Figure 14), the thickness of the mass on the left part of the domain is the same, but the We show here that the scenario with the gate better fits the dynamics of the experimental collapse. Interestingly, the final deposit is however the same in the two scenarios: the gate has an influence on the flow dynamics that happens to disappear on the deposit. The same qualitative behavior is observed for granular collapse over inclined beds (Figure 15). The influence of the gate is even higher at θ = 16 o than at θ = 0 o . These results are in good agreement with the simulations of Carbopol flows over a horizontal plane performed by Martin Rentschler [51]. The influence of the gate may explain why Lacaze et al. [START_REF] Lacaze | Planar collapse of a granular column: experiments and discrete-element simulations[END_REF] obtained a faster front propagation in their discrete element simulations than in their experiments while reproducing the experimental runout distance, even though their way of removing the gate was different (see their Figure 4). However, for a very slow gate removal the runout distance is significantly shorter (see Mériaux [START_REF] Mriaux | Two dimensional fall of granular columns controlled by slow horizontal withdrawal of a retaining wall[END_REF]).

Influence of the friction on the left wall and at the base

Increasing the friction on the left wall from µ w = tan 10.5 o to µ w = tan 25.5 o in the simulations has only very little influence on the flow. In particular, it does not improve the overestimation of the collapse on the left part of the granular mass. Indeed as discussed before, the pressure near the half upper part of the left wall is small, so that the friction force is still always small whatever the friction coefficient at the wall. Near the bottom of the left wall, the pressure is higher but their is no motion. In the experiments, the presence of lateral walls may significantly change this picture. Indeed the presence of the lateral walls may increase the pressure in their vicinity, which may increase the friction force near the left wall and thus reduce the vertical collapse of the mass in its vicinity.

Our simulations indicate that part of the mass is sliding at the base near the front for a basal friction coefficient µ b = 0.48 (see Figures 2 and8), contrary to what was assumed in other simulations where no-slip (i.e. adherence) was prescribed at the bottom (e.g. Lagree et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF], Crosta et al. [START_REF] Crosta | Numerical modeling of 2-D granular step collapse on erodible and non-erodible surface[END_REF]). To investigate the influence of the basal boundary condition we have prescribed at the base either real adherence (zero-velocity) or a very high friction coefficient. Imposing adherence reduces the t = 0.3 s t = 1.02 s collapse of the mass near the left wall leading to a maximum thickness of the deposit h f about 4% higher than with the a basal friction coefficient µ b = 0.48 (see Figure 16). It also reduces the runout distance by 10%. As a result, the simulation with adherence better fits the results near the left wall but underestimates the runout distance. If, instead of prescribing adherence, we impose a very high friction at the bottom (µ b = 4.8) we obtain a maximum thickness near the left wall and a runout distance that are in between the simulation with adherence and that with a basal friction µ b (Figure 17). In fact, near the front the normal stress is very small and even with a very high friction coefficient the friction force is still small enough for the balance between this dissipative force and the driving forces to make it possible for the material to slip on the basal surface. As a result, sliding at the base is allowed near the front, contrary to the case of adherence. Note that in Lagrée et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF], the mass spreads less rapidly in their continuum simulation than in their DEM simulations. One explanation could be that in their continuum model, they imposed a no-slip boundary conditions while there is a significant slip velocity near the front. 5 Insight into the flow dynamics

Velocity distribution

Figure 18 shows the distribution of the horizontal velocity at θ = 0 • as well as its vertical profiles at some fixed locations along the plane (at x = 10, 15, 30 and 40 cm) (pink lines). In the first snapshot, at t = 0.06 s, just after the gate is lifted, the maximum horizontal velocity is reached at depth, close to the bottom. This is the effect of the gate removal which induces a flow starting next to the bottom. At later time, the flow is concentrated only in a shallow sub-domain near the free surface as observed with the norm of the velocity field. The maximum horizontal velocity is reached near the free surface. In the flowing region, the horizontal velocity has a parabolic profile with an exponential decay near the flow/no-flow transition. Near the front, the granular mass slips at the contact with the bottom surface, i.e. the basal velocity is not equal to zero (for example at t = 0.22 s and x = 30 cm or at t = 0.35 s and x = 40 cm). Figure 18 shows that the horizontal velocity could change direction (i.e. negative horizontal velocity) in the quasi-static region (e.g. at t = 0.18, 0.22, 0.3 s). However these negative velocities are small and close to the expected numerical error.

The maximum downward and upward vertical velocities are about 0.6 m.s -1 and 0.06 m.s -1 , respectively (Figure 19). As a result, the downward vertical velocity can be as high as half of the horizontal velocity. Figure 19 clearly shows, at t = 0.18 s, an interface between a zone with no vertical velocity and a zone with a small downward velocity. Under this rough interface, a localized zone of upward velocity develops (dark red color and negative vertical velocity in Figure 19 at t = 0.18, 0.3 s). Within the small vertical velocity regions underlying the main flowing layer (red colors), the vertical velocity distribution is not regular and exhibit spatio-temporal fluctuations at a scale of a few centimeters or more.

Yield limit distribution and stagnant zones

The yield limit κ = µ s p, which has the same distribution as the pressure p (plotted in Figure 4 (right)), has an almost gravitational distribution. The stress deviator distribution σ ′ , plotted in Figure 20 (right), is rather different. The maximum value is reached in the middle part, while near the free surface and in the left corner the stress deviator is decreasing rapidly. Except at the left bottom corner, we remark that the yield stress and the stress deviator have the same order of magnitude (i.e. the same color on the color scale) which means that the yielding surface separating flow and no-flow zones could change a lot during the experiment. For t = 0.3 s and t = 0.42 s the yield limit (pressure) distribution presents an unexpected zone of weakness (low pressure). This (non-static) sucking zone, which is localized on a slope line, appears before the granular mass reached a configuration close to equilibrium and it is associated to a slightly positive vertical velocity. The static/flowing transition can be seen in the left side of Figure 20, where we have plotted the distribution of the yield function F (σ) := σ ′ -µ s p. The zone corresponding to a positive yield function (dark red) represents the flowing region of the granular mass, while the strain rate evolution is plotted in Figure 4 (left). We remark that for t = 0.18 s and t = 0.30 s the regions with a high strain rate are located at the left of the tip of the granular mass, near the bottom surface. In these regions the inertial and gravitational forces have to be accommodated with the frictional forces giving an important shear of the granular mass. As we can see by comparison with Figure 20 (right) these high deformation regions are not related to high stress deviator zones. This confirms the nonlinear character of the constitutive law [START_REF] Bouchut | A shallow model including staticflowing transition for viscoplastic Drucker-Prager materials[END_REF].

For t = 0.3 s or for t = 0.42 s the yielding surface is not sharply determined. Indeed, we remark that at the left of the main flowing region (on the right of the granular mass near the free surface) there are some zones corresponding to a nonvanishing strain rate in Figure 4 (left). These small regions are related to the "weakness zones" of low yield limit (or pressure) described above. Note that in all the computations presented above we have not used any special techniques for tracking the yield surface. For a sharp investigation of the static/flowing transition one should use some of the re-meshing techniques developed for stationary flows (see for instance [START_REF] Saramito | An adaptive finite element method for viscoplastic fluid flows in pipes[END_REF][START_REF] Cazacu | Steady-state flow of compressible rigidviscoplastic media[END_REF]) at each time step. This will imply however an important increase of the computational effort. deposit behind the front is made of the layers initially located in the half deeper part. This may be different for columns with high aspect ratios (see Thomson and Huppert [START_REF] Thompson | Granular column collapses: further experimental results[END_REF]). This kind of observation may be important to help interpretation of real deposits. The deformation of the artificial initially vertical layers is very similar to that obtained by Lagree et al. [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF] that used the µ(I) rheology (see their figure 9), except near the base where the boundary conditions are different (sliding friction here and no-slip in their simulations). As a result, the deformations of the inner vertical layers are not specifically related to the µ(I) rheology.

Deformation of horizontal and vertical layers

Discussion and conclusion

We have proposed a mechanical and numerical model for granular flows based on the Drucker-Prager plasticity (pressure dependent yield stress) with either a constant viscosity or a variable viscosity η( D , p), depending on both the pressure p and the norm of the strain rate tensor D , that comes from the reformulation of the so-called µ(I) rheology.

We showed that the Drucker-Prager plasticity criterion and a variable or constant viscosity well reproduce quantitatively the dynamics and deposit of the granular collapse of columns with small aspect ratio (here a = 0.7) over horizontal and inclined planes, using rheological parameters derived from the experiments, without any fitting procedure. Note that this is obtained while taking into account the initial removal of the gate that releases the granular column. The Drucker-Prager rheology with a constant viscosity (η = 1 Pa.s) and the variable viscosity η( D , p) (i.e. the µ(I) rheology) provide very similar results. The numerical simulation shows that the main differences are located near the front. Indeed, the parameter I of the µ(I) rheology mainly varies in this region. Interestingly, the region near the front is characterized by thin flows over a rigid bed, a situation close to that where the µ(I) rheology was established. While the variable viscosity η( D , p) rheology may be crucial to simulate different regimes such as the collapse of high aspect ratio columns (Lagree et al., [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF]), the strong similarity of the two rheologies in the case studied here is an important result in terms of application to natural landslide, that have typical aspect ratios lower than 1 (see Lucas et al., [START_REF] Lucas | Frictional weakening in landslides on Earth and on other planetary bodies[END_REF] for a compilation of real data). An intermediate linearized approach [START_REF] Hughes | Lagrangian-Eulerian finite element formulation for incompressible viscous flows[END_REF], [START_REF] Ionescu | Onset and dynamic shallow flow of a viscoplastic fluid on a plane slope[END_REF] can also be used. Indeed, for the Drucker-Prager rheology with constant viscosity, only 2 rheological parameters are needed, compared to 3 for the variable viscosity (and 2 for the linearized law). Furthermore, calculations with a constant viscosity η = 1 Pa.s requires much less computational time (at least for the numerical approach presented here). A crucial contribution of the expression of the viscosity η( D , p) (equation ( 14)) derived from the µ(I) rheology is to provide the order of magnitude of the viscosity in a specific situation, without having to perform the simulation. As an example, equation [START_REF] Duarte | Arbitrary LagrangianEulerian method for NavierStokes equations with moving boundaries[END_REF] gives the value of the viscosity (η = 1 Pa.s here) that makes it possible to reproduce the experimental results using a constant viscosity, while numerical simulations show that the variable viscosity ranges from values up to η = 1.5 Pa.s in the deep quasistatic regions down to zero near the free surface. Discriminating the constant and variable viscosity rheologies in the regimes of granular collapse of small aspect ratio on inclined beds would necessitate more detailed comparison with measured velocity profiles.

The gate has a significant impact on the flow dynamics but the deposits obtained with and without the gate are almost the same. The deposit seems thus to depend only on the initial energy of the column and on the mean dissipative forces at work during the flow, provided gate removal is not too slow (Mériaux [START_REF] Mriaux | Two dimensional fall of granular columns controlled by slow horizontal withdrawal of a retaining wall[END_REF]).

As the slope of the bed increases, these models slightly overestimate the flow velocity and runout distance. This may be due to lateral wall friction that increases with increasing thickness of the flowing layer at higher slopes. This effect is not accounted for in the 2D simulations. As in discrete element modeling, 3D viscoplastic modeling of the wall effect will certainly increase the dissipation due to the wall as the slope increases but will also increase the pressure, in particular in the vicinity of the upper left wall of the reservoir, thus increasing the friction in this region. This will certainly reduce the overestimated collapse of the columns near the left wall. On the other hand, simulations show that there is a significant sliding of the mass at the bottom near the front. As a result, prescribing adherence near the base as was done in previous studies seems not appropriate. To show that we have investigated the influence of the friction on the left wall and at the base. To do that we have prescribed at the base either real adherence (zero-velocity) or a very high friction coefficient. Imposing adherence at the base reduces the collapse of the mass near the left wall but also reduces the runout distance by 10%. As a result, the simulation with adherence better fits the results near the left wall but underestimates the runout distance. By choosing a very high friction coefficient at the bottom, sliding is allowed near the front and the computed runout distance is in between the simulation with adherence and that with a normal friction. Indeed, near the front the normal stress is very small and even with a very high friction coefficient the friction force is still small enough for the balance between this dissipative force and the driving forces to allow basal sliding.

Except in the very first moments during the gate removing period, the yield limit has an almost gravitational distribution while the stress deviator reaches its maximum value at mid flow height. The "constitutive" flowing region (i.e. with a positive yield function) near the free surface is globally well delimitated and corresponds to a non-vanishing velocity and to a non-vanishing strain rate zones. However, during the flow the yield limit distribution presents a zone of weakness (low pressure), localized on a slope line and associated to a slightly positive vertical velocity.

The regions with a high strain rate are located at the left of the tip front of the granular mass, near the bottom surface, where an important shear of the granular mass is expected as a the result of the action of the driving (inertial and gravitational) forces and the frictional forces. These high deformation regions are not related to high stress deviator zones, confirming the strong nonlinear character of the viscoplastic law.

The numerical simulation of the 2D collapse of granular columns shows that only a shallow layer of granular material is flowing. Despite the heterogeneities observed within the flow and especially near the static/flowing transition, essentially, the velocity profiles have a typical shape with a maximum velocity at the free surface, a Bagnold-like to linear profile and an exponential decrease of the velocity near the static/flowing transition. The pressure seems to be close to hydrostatic in this flowing region. These results provide new constrains for developing thin layer models with a flowing layer overtopping a static zone (see e.g. Bouchut et al. [START_REF] Bouchut | A shallow model including staticflowing transition for viscoplastic Drucker-Prager materials[END_REF][START_REF] Lusso | A simplified model of thin layer static/flowing dynamics for granular materials with yield[END_REF]).

Step 3. Finally, we update the stress deviator σ ′k,n and the tangential stress σ k,n T through σ ′k,n = σ ′k,n-1 + r(D(u k,n ) -γk,n ),

σ k,n T = σ k,n-1 T -r f (u k,n T -δ k,n ).
Solving the Stokes type problem at step 1 is a standard problem in fluid mechanics and there exists many technics to do it (see for instance [START_REF] Pironneau | Finite element methods for fluids[END_REF]). Here, D is discretized by using a family of triangulations (T h ) h made of finite elements (h > 0 is the discretization parameter representing the greatest diameter of a triangle in T h ). We denote by V h the FE space for the velocity field u k,n , by W h the FE space for the pressures field p k,n , by Q h the FE space for the stresses deviators σ ′k,n and for the strain rate multipliers γk,n and R h the space for the for slip rate multipliers δ k,n and for tangential stresses σ k,n T . Note that these finite element spaces cannot be chosen independently. For instance, in two dimensional computations, if V h =[continuous P2] (i.e. the space of continuous functions which are polynomial of degree 2 on each triangle) then we have to choose W h =[continuous P1], Q h =[discontinuous P1] and R h = [continuous P2]. This choice is used in all simulations presented in this paper.

Note that if the convergence is achieved then the iterative solution of the algorithm satisfies the nonlinear system ( 26)- [START_REF] Lacaze | Planar collapse of a granular column: experiments and discrete-element simulations[END_REF]. Indeed, if the convergence criterion σ ′k,n = σ ′k,n-1 , σ k,n

  Following the stress decomposition we get the constitutive scalar function β( D , p) = 2η( D , p) + κ(p)/ D . Now we can compute σ ′ = 2η( D , p) D + κ(p) and (1) reads trace(D) = 0,

Figure 1 :

 1 Figure1: Experimental setup: morphometric and control parameters measured in the experiments. The initial mass (light gray) with initial thickness h 0 = 14 cm and width r 0 = 20 cm is released on a plane with inclination θ by opening very rapidly a gate at time t = 0 s. It forms a deposit with a length r f from r 0 (runout distance), and a final maximum thickness h f . The gate is represented by a thick line perpendicular to the plane.

  should be µ s ≃ 0.38 + 0.1 = 0.48 = tan(25.6 o ) and µ 2 ≃ 0.64 + 0.1 = 0.74 = tan(36.5 o ). Finally, we chose here µ s = tan(25.5 o ) ≃ 0.48 and µ 2 = tan(36 o ) ≃ 0.73.

  48 (25.5 o ) 0.48 (25.5 o ) 0.18 (25.5 o ) 0.73 (36 o ) 0.279 0.035 kg 1/2 .m -1/2

Figure 2 :

 2 Figure 2: Simulation using the variable viscosity η( D , p) (i. e. µ(I)): comparison of the simulated granular mass and the experimental results (pink line) at different times for granular collapse over an horizontal plane (θ = 0 • ). The colors represent the distribution of the norm of the computed velocity field |u| (in m/s).

Figure 3 :

 3 Figure 3: Viscosity (in Pa.s) calculated with the µ(I) rheology at different times for granular collapse over an horizontal plane (θ = 0 • ).

Figure 4 :

 4 Figure 4: (left column) Strain rate D (in s -1 ), and (right column) pressure p (in Pa), calculated with the variable viscosity η( D , p) (i. e. µ(I)) at different times for granular collapse over an horizontal plane θ = 0 • .

Figure 5 :

 5 Figure 5: Values of I and of µ(I) at different times for granular collapse over an horizontal plane θ = 0 • .

Figure 6 :

 6 Figure 6: Comparison between the simulations using the variable viscosity η( D , p) (i.e. µ(I)) (light gray) and the constant viscosity with η = 1 Pa.s (blue mesh) for granular column collapse at θ = 0 o . The experimental results are represented in black lines. The inserts at the right side correspond to the zooms near the front of the mass.

.

Figure 7 :

 7 Figure 7: Comparison between the simulations using the variable viscosity η( D , p) (i.e. µ(I)) (light gray) and the constant viscosity with η = 0.1 Pa.s (blue mesh) for granular column collapse at θ = 0 o . The experimental results are represented in black lines. The inserts at the right side correspond to the zooms near the front of the mass.

.

Figure 8 :

 8 Figure 8: Comparison between the simulations using the variable viscosity η( D , p) (i.e. µ(I)) (light gray) and the constant viscosity with η = 10 Pa.s (blue mesh) for granular column collapse at θ = 0 o . The experimental results are represented in black lines.

Figure 9 :

 9 Figure 9: Simulation using the variable viscosity η( D , p) (i.e. µ(I)): comparison of the simulated granular mass and the experimental results (pink line) at different times for granular collapse over a plane inclined at θ = 16 • . The colors represent the distribution of the norm of the computed velocity field |u| (in m/s).

Figure 10 :

 10 Figure 10: Viscosity (in Pa.s) calculated with the µ(I) rheology at different times for granular collapse over a plane inclined at θ = 16 • .

Figure 11 :Figure 12 :

 1112 Figure 11: (Left column) Strain rate (in s -1 ), and (Right column) Pressure (in Pa), calculated with the variable viscosity η( D , p) (i.e. µ(I)) at t = 0.15 s, t = 0.23 s, and t = 0.56 s for granular collapse over a plane inclined at θ = 16 • .

Figure 13 :

 13 Figure 13: Comparison between the simulations using the variable viscosity η( D , p) (i.e. µ(I)) (light gray) and the constant viscosity with η = 1 Pa.s (blue mesh) at different times for granular collapse over a plane inclined at θ = 16 • . The experimental results are represented in black line. The inserts at the right side correspond to the zooms near the front of the mass.

Figure 14 :

 14 Figure 14: Influence of the gate on granular collapse over a horizontal plane θ = 0 • . Comparison of the thickness of the granular mass simulated with the gate (light blue) and without the gate (blue mesh), as well as the experimental results (pink line) at different times.

Figure 15 :

 15 Figure 15: Influence of the gate on granular collapse over an inclined bed θ = 16 • . Comparison of the thickness of the granular mass simulated with the gate (light blue) and without the gate (blue mesh), as well as the experimental results (pink line) at different times.

Figure 16 :

 16 Figure 16: Influence of the boundary conditions on the base on granular collapse over a horizontal plane θ = 0 • . Comparison of the thickness of the granular mass simulated in the "no gate" scenario with friction (dark blue) and with adherence (light blue mesh), as well as the experimental results (black line) at different times.

t

  = 0.3 s t = 1.02 s

Figure 17 :

 17 Figure 17: Influence of the boundary conditions on the base on granular collapse over a horizontal plane θ = 0 • . Comparison of the thickness of the granular mass simulated in the "gate" scenario with friction (light blue) and with high friction (blue mesh), as well as the experimental results (black line) at different times.

Figure 18 :

 18 Figure 18: The evolution of the distribution of the horizontal velocity u x (color scale in m/s) and vertical profiles (pink lines) of the horizontal velocity u x (z) at four fixed locations x = 10, 15, 30 and 40 cm along the plane, for simulation with the Drucker-Prager model with constant viscosity η = 1 Pa.s over horizontal plane θ = 0 • . Note that the vertical profiles are represented in the plane transverse to the simulated flow.

Figure 19 :

 19 Figure 19: The evolution of the distribution of the vertical velocity u y (color scale in m/s) and vertical profiles (blue lines) of the vertical velocity u y (z) at three fixed locations x = 10, 15 and 30 cm along the plane, for simulation with the Drucker-Prager model with constant viscosity η = 1 Pa.s over horizontal plane θ = 0 • . Note that the vertical profiles are represented in the plane transverse to the simulated flow.

Figure 21 showsFigure 20 :

 2120 Figure21shows how artificial initially horizontal and vertical layers deform during the flow. At time t = 0.18 s, the left part (next to the left wall) of the top horizontal layers are collapsing vertically more rapidly than between x = 0.05 and 0.1 m. Going further right, the top layers get thinner and kind of converge toward a stamp point at the surface, corresponding to the advancing upper right corner of the initial column. Under and next to this point, the horizontal layers are bending, first going down while thinning and then going up while thickening towards the front. Interestingly, the left 2/3 of the deposit is made of the same succession of layers than the initial mass, even though these layers are deformed. On the contrary, the right 1/3 of the

Figure 21 :

 21 Figure 21: Deformation of the initially horizontal layers (left) and initially vertical layers (right) calculated for granular collapse on horizontal plane θ = 0 • (color scale in m).

Table 1 :

 1 Rheological parameters.

Acknowledgements

The authors are thankful to Yoël Forterre, Christophe Ancey, Guillaume Ovarlez, John de Bruyn, and Maxime Farin for fruitful discussions. This work has been partially funded by the ANR contract ANR-11-BS01-0016 LANDQUAKES, the USPC project PAGES, and the European Research Council ERC contract ERC-CG-2013-PE10-617472 SLIDEQUAKES.

Appendix-Numerical Method

We propose here a numerical algorithm to solve the dynamic flow problem described above. We use a time implicit (backward) Euler scheme for the time discretization of the field equations, which gives a set of nonlinear equations for the velocity u, the deviatoric stress tensor σ ′ and the pressure p. At each iteration in time, an iterative algorithm is used to solve these nonlinear equations. Note that in the case of the proposed viscoplastic model, numerical difficulties arise from the nondifferentiability of the viscoplastic and friction terms. That means that one cannot make use of the numerical techniques developed for Navier-Stokes fluids. To overcome these difficulties the iterative decomposition-coordination formulation coupled with the augmented Lagrangian method of [START_REF] Glowinski | Augmented Lagrangian and Operator Splitting method in Non-Linear Mechanics[END_REF][START_REF] Fortin | Méthodes de Lagrangien augmenté, application à la résolution de problèmes aux limites[END_REF] is adapted here.

To include frictional effects in the algorithm we have first regularized the friction law [START_REF] Glowinski | Augmented Lagrangian and Operator Splitting method in Non-Linear Mechanics[END_REF]. This can be done (see also [START_REF] Ionescu | Onset and dynamic shallow flow of a viscoplastic fluid on a plane slope[END_REF][START_REF] Ionescu | Augmented Lagrangian for shallow viscoplastic flow with topography[END_REF]) by introducing a small frictional viscosity η f << 1 in the Coulomb friction law and writing [START_REF] Glowinski | Augmented Lagrangian and Operator Splitting method in Non-Linear Mechanics[END_REF] as

where, as before, [ ] + is the positive part. Note that using this regularization, the friction law has the same mathematical structure as the visoplastic constitutive equation ( 8) and we can use the same iterative decomposition-coordination formulation. For all the computations of the present paper the regularization frictional viscosity coefficient was chosen to be η f = 0.1 Pa.s.m -1 . An ALE (Arbitrary Lagrangian-Eulerian) description is used to compute the evolution of the fluid domain D. For the sake of simplicity we will not present here the ALE formulation of the numerical scheme. As a mater of fact, the numerical algorithm proposed here makes use only of a Stokes-type problem at each time step, and the implementation of the Navier-Stokes equations in an ALE formulation is rather standard (see for instance [START_REF] Hughes | Lagrangian-Eulerian finite element formulation for incompressible viscous flows[END_REF][START_REF] Maury | Characteristics ALE method for the unsteady 3D Navier-Stokes equations with a free surface[END_REF][START_REF] Maronnier | Numerical simulation of three dimensional free surface flows[END_REF][START_REF] Duarte | Arbitrary LagrangianEulerian method for NavierStokes equations with moving boundaries[END_REF]).

Time discretization

Let ∆t be the time step and let us denote by u k , σ ′k and p k the values of the unknowns u(k∆t), σ ′ (k∆t) and p(k∆t). Suppose that we have computed all these variables at time t = (k -1)∆t.

The time implicit (backward) Euler scheme for the field equations ( 15), ( 16) and ( 8) gives the following nonlinear equations for the unknowns u k , σ ′k and p k ρ(

while the boundary conditions read

The algorithm at each time step

Let us fix the iteration in time, k. In order to describe the algorithm let r, r f > 0 be the augmented Lagrangian coefficients. We start with

T and let u k,n-1 , σ ′k,n-1 and σ k,n-1 T be known. We suppose also that the strain rate multipliers γk,n-1 : D → R 3×3 S , and the slip rate multipliers δ k,n-1 : Γ b → R 2 introduced below, are also known. The convergence is achieved when the difference between u k,n , σ ′k,n , p k,n , σ k,n T and u k,n-1 , σ ′k,n-1 , p k,n-1 , σ k,n-1 T is small enough.

Step 1. The first step consists in solving the following linear equation of Stokes type for the velocity field u k,n and the pressure p k,n :

with the boundary conditions

Step 2. First we update the viscosity coefficient η = η( D(u k,n ) , p k,n ) and the yield limit κ = κ(p k,n ). Then, we compute the strain rate multipliers γk,n and the slip rate multipliers δ k,n γk,n = 1 2η + r 1 -κ σ ′k,n-1 + rD(u k,n ) + (σ ′k,n-1 + rD(u k,n )), ( 33)

according to the decomposition-coordination formulation coupled with the augmented Lagrangian method.