A free interface model for static/flowing dynamics in thin-layer flows of granular materials with yield: simple shear simulations and comparison with experiments

Flows of dense granular materials comprise regions where the material is ﬂowing, and regions where it is static. Describing the dynamics of the interface between these two regions is a key issue to understand the erosion and deposition processes in nat-ural environments. A free interface simpliﬁed model for non-averaged thin-layer ﬂows of granular materials has been previously proposed by the authors. It is a coordinate-decoupled (separated variables) version of a model derived by asymptotic expansion from an incompressible viscoplastic model with Drucker-Prager yield stress. The free interface model describes the evolution of the velocity proﬁle as well as the position of the transition between static and ﬂowing material. It is formulated using the coordinate Z in the direction normal to the topography and contains a source term that represents the opposite of the net force acting on the ﬂow, including gravity, pressure gradient, and internal friction. In this paper we introduce two numerical methods to deal with the particular formulation of this model with a free interface. They are used to evaluate the respective role of yield and viscosity for the case of a constant source term, which corresponds to simple shear viscoplastic ﬂows. Both the analytical solution of the inviscid model and the numerical solution of the viscous model (with a constant viscosity or the variable viscosity of the µ ( I ) rheology) are compared with experimental data. Although the model does not describe variations in the ﬂow direction, it reproduces the essential features of granular ﬂow experiments over an inclined static layer of grains, including the stopping time and the erosion of the initial static bed, which is shown to be closely related to the viscosity for the simple shear case.


Introduction
Understanding and theoretically describing the static/flowing transition in dense granular flows is a central issue in research on granular materials, with strong implications for industry and geophysics, in particular in the study of natural gravity-driven flows.Such flows (e.g., landslides or debris avalanches) play a key role in erosion processes on the Earth's surface and represent major natural hazards.In recent years, significant progress in the mathematical, physical, and numerical modeling of gravity-driven flows has made it possible to develop and use numerical models for the investigation of geomorphological processes and assess risks related to such natural hazards.However, severe limitations prevent us from fully understanding the processes acting in natural flows and predicting landslide dynamics and deposition, see e.g., [1].In particular, a major challenge is to accurately describe complex natural phenomena such as the static/flowing transition.
Geophysical, geotechnical, and physical studies have shown that the static/flowing transition related to the existence of no-flow and flow zones within the mass plays a crucial role in most granular flows and provides a key to understanding their dynamics in a natural context, see [1] for a review.This transition occurs in erosion-deposition processes when a layer of particles flows over a static layer or near the destabilization and stopping phases.Note that natural flows often travel over deposits of past events, which may or may not be made of the same grains, and entrain material from the initially static bed.Even though erosion processes are very difficult to measure in the field, e.g., [2][3][4], entrainment of underlying material is known to significantly change the flow dynamics and deposition, e.g., [5][6][7][8][9][10].
Experimental studies have provided some information on the static/flowing dynamics in granular flows, showing for example that the presence of a very thin layer of erodible material lying on an inclined bed may increase the maximum runout distance of a granular avalanche flowing down the slope by up to 40% and change the flow regimes [9,11,12].In these experiments that mimic natural flows over initially static beds (Figure 1), quasi-uniform flows develop when the slope angle θ is a few degrees lower than the typical friction angle δ of the involved material [12].Figures 2 and 3 show new data extracted from the experiments performed by [12] on the change with time of the static/flowing interface position b and the velocity profiles U(Z) within the granular mass, where Z is the direction normal to the bed.At a given position X along the plane, the flow is shown to excavate the initially static layer immediately when the front reaches this position.The static/flowing interface rapidly penetrates into the static layer, reaching a lowest position (that depends on the slope angle) for significantly high slopes and then rises almost linearly or exponentially toward the free surface until the whole mass of material stops (Figure 2).A theoretical description of these observations is still lacking.In particular, questions remain as to what controls the change with time of the static/flowing interface position and the velocity profiles and how these characteristics are related to the rheology of the granular material and the initial and boundary conditions.[9,12] to study granular column collapse over inclined planes covered by an initially static layer made of the same grains as those released in the column.For slope angles θ a few degrees smaller than the typical friction angle δ of the involved material, a quasi-uniform flow develops behind the front., and θ = 24 • (black vertical crosses).Time t = 0 s corresponds to the time when the front of the flowing layer reaches the position X = 90 cm.The approximate upward velocity ḃ of the static/flowing interface is indicated in m/s for each slope angle.These new results have been extracted from the experiments performed by [12,13] for granular columns of initial down-slope length r 0 = 20 cm, initial thickness h 0 = 14 cm, and width W = 20 cm (i.e., volume V = 5600 cm 3 ).Note that the position of the free surface when the mass stops, i.e., when b = h, represented by the upper point for each angle, slightly depends on the slope angle and decreases as the slope angle increases, from about 0.025 m at θ = 19 • to about 0.018 m at θ = 24 • . .Time t = 0 s corresponds to the time when the front of the flowing layer reaches the position X = 90 cm.These new results have been extracted from the experiments performed by [12,13] for granular columns of initial down-slope length r 0 = 20 cm, initial thickness h 0 = 14 cm, and width W = 20 cm (i.e., volume V = 5600 cm 3 ).
In order to alleviate the high computational costs required to describe the real topography and the rheology, which both play a key role in natural flow dynamics, thin-layer (i.e., the thickness of the flow is assumed to be small compared to its down-slope extension) depth-averaged models are generally used to simulate landslides [14].Such models have been rigorously derived for arbitrary topography, but the static/flowing transition has generally been neglected, e.g., [15][16][17][18].Several attempts have been made to describe this transition in thin-layer (i.e., shallow) models, in particular by deriving an equation for the static/flowing interface position or by establishing erosion/deposition rates, e.g., [19,20].However, these approaches are generally based on debatable phenomenological laws and/or are too schematic to be extended to natural flows [21].The simplifications used in some previous models lead to an inconsistent energy Equation [22].Indeed, some of these models prescribe a given velocity profile to deduce the entrainment rate (see e.g., [23] or references in [22]).However, in transient flows, the velocity profile changes with time as shown for example in the multilayer shallow model of [24].A better understanding of the non-averaged case is necessary before defining more physically relevant depth-averaged models including the static/flowing transition.This is why we focus here on non-averaged thin-layer models.
Recent work has shown that viscoplastic flow laws with yield stress describe well granular flows and deposits in different regimes, from steady uniform flows [25][26][27][28] to transient granular collapse over rigid or erodible horizontal beds [29][30][31] and inclined beds [24,32,33] or accelerating/decelerating channel flows [34].Very good agreement with experiments is obtained even though these laws, and in particular the so-called µ(I) rheology, where I is the inertial number [28], have been shown to be ill-posed in the quasi-static regime, i.e., for small I (and also for large I) [35].These promising results may be related to the use of a coarse mesh, whereby simulations avoid the ill-posedness by damping the fast-growing high wavenumbers [33].In any case, the quasi-static regime near the static/flowing interface is known to be very complex, involving strong and weak force chains and local rearrangement of particles, e.g., [36,37] and has been widely studied, in particular using discrete element methods, e.g., [1,38,39].This regime is not accurately described by the proposed viscoplastic laws involving a simple yield stress.Questions however remain as to whether such 'simple' viscoplastic laws are able to describe quantitatively the change with time of the static/flowing interface position and the velocity profiles observed experimentally for flows over an initially static bed and how the viscosity affects these processes.
Based on such a viscoplastic model with yield stress, an analytic expansion in the non-averaged thin-layer regime is provided in [40], giving a theoretical basis for equations describing the static/flowing interface dynamics in a dry granular material.A key point in this approach that fundamentally differs from previous thin-layer models (see [21] for a review) is that even if the flow is assumed to be thin, the normal coordinate Z is still present (there is no depth-averaging).The static/flowing interface is free in this approach, in the sense that it has no explicit equation defining its evolution.Its time variation is instead implicitly determined by an extra boundary condition on the velocity at this interface.The model from [40] being still too complicated, however a formal decoupling of the coordinates X in the down-slope direction and Z normal to the topography is proposed in [41], leading to a model with only the coordinate Z, but including a source term S that represents the opposite of the net force acting on the flow, including gravity, pressure gradient, and internal friction.We propose here to evaluate the dynamics produced by the free interface model of [41], for the most simple case involving a constant source S.This corresponds to assuming no dependency on the down-slope coordinate X in the model of [40].It is also equivalent to considering simple shear solutions to the original viscoplastic model, which shows that in this case the model is valid without a smallness assumption.Although the down-slope coordinate X plays a key role in real flows, taking into account the fluctuations of the free surface, topography, and inflow information, we must first understand the dynamics when X is not involved.We show that the analysis of this simple shear system provides new insight into the change of the static/flowing interface position and the velocity profiles with time.
Because the boundary conditions and interface evolution are formulated in an uncommon way, specific methods must be used to deal numerically with the free interface model.We introduce two numerical methods that give similar grid-independent results, thereby showing that our model is well-posed for the case considered here involving only Z dependency.
The paper is organized as follows.The model of [40,41] for the case of a constant source (corresponding to simple shear flows of the viscoplastic model) is recalled in Section 2.Then, in Section 3, we derive an analytical solution for the inviscid model with a constant source term, which partly reproduces the experimental observations and shows explicitly how the static/flowing interface position and the velocity profiles are related to the flow characteristics and to the initial and boundary conditions.In Section 4, we introduce two new numerical methods for the simulation of the viscous model, with a constant viscosity or the variable viscosity associated with the µ(I) rheology.The numerical results show that, as opposed to the inviscid model, viscosity makes it possible to reproduce the initial penetration of the static/flowing interface within the static bed and the exponential shape of the velocity profiles near this interface.Finally, in Section 5, these results are discussed based on comparison with former experimental and numerical studies on granular flows, showing the key features and the limits of the non-averaged thin-layer model for providing a better understanding and modeling of laboratory and natural flows.

Origin of the Model: Viscoplastic Rheology with Yield Stress
The initial governing equations express the dynamics of an incompressible (∇ • U = 0) viscoplastic material with Drucker-Prager yield stress, with the rheological law where U is the velocity, −g gravitational acceleration, σ the stress tensor (normalized by the density), p the scalar pressure (also normalized by the density), and DU the strain rate tensor DU = (∇U + (∇U) t )/2.Here the norm of a matrix . The coefficients µ s > 0 and ν ≥ 0 are the internal friction and the kinematic viscosity, respectively.In general ν and µ s could depend on DU and p, but we shall consider µ s constant here since the so-called µ(I) law can be written in this form for some viscosity ν, see Section 4.5.

Free Interface Model with Source Term
A non-averaged thin-layer model for flows described by Equations ( 1) and (2) with a static/flowing interface has been derived in [40].It is formulated for the coordinates X in the direction tangent to the topography and Z normal to the topography.The topography is described by its angle θ(X) with respect to the horizontal (see Figure 4 for a flat topography).For the case of simple shear (no dependency on X), no smallness assumption is necessary ([40] Section 2.3) and the model can be written as the following free interface problem with a source term.The unknowns are U(t, Z) (velocity in the direction tangent to the topography) defined for b(t) < Z < h, where h > 0 is the given constant thickness of the layer and b(t) the unknown position of the static/flowing interface, 0 < b(t) < h.They obey the following equations, with the following boundary conditions for all t > 0, The source S(t, Z) is assumed to be given and represents the opposite of the net force, excluding the viscous force (see [41] and (10) below).The viscosity ν can a priori depend on t and Z in the model.The velocity is extended by defining U(t, Z) = 0 for 0 ≤ Z ≤ b(t), corresponding to the static domain.The initial condition is formulated as where U(0, Z) is the limit of U(t, Z) as t ↓ 0. The initial velocity U 0 (Z) is given for Z ∈ [0, h], such that for some b 0 ∈ (0, h) (initial position of the interface), the function and ∂ Z U 0 > 0 for Z ∈ (b 0 , h).The system is completed by the static equilibrium condition that states that the net force must be resistive in the static layer.The velocity U obtained when solving (3)-( 6) is expected to satisfy ∂ Z U > 0 for b(t) < Z < h.As proved in [41], this is the case as soon as ∂ Z S ≤ 0 and ν > 0.
The unusual nature of the system (3)-( 6) lies in the fact that the interface position b(t) is free, meaning that it has no given equation relating ḃ(t) ≡ db/dt to other quantities.Its evolution is instead implicitly governed by the boundary conditions (4a) and (4b).An interesting property of the model ( 3)-( 6) is that a differential equation can nevertheless be derived for the time evolution of the position of the interface, valid under some conditions, as follows (see [41]).

•
If ν > 0 then and whenever S(t, b(t)) = 0, one has These formulas show the strong interrelationship between the velocity profile in the direction perpendicular to the inclined plane and the evolution of the static/flowing interface position.Erosion, i.e., penetration of the static/flowing interface into the erodible bed, occurs when ḃ(t) < 0, and deposition when ḃ(t) > 0.
The full model of [40] is indeed more general than (3)-( 6), since it includes a possible dependency on X.It is derived under the assumptions that the thickness of the layer, the curvature of the topography, and the viscosity are small, the internal friction angle is close to the slope angle, the velocity is small, and the pressure is convex with respect to the normal coordinate Z. Then an intermediate model with formal decoupling of the X and Z coordinates has been proposed in [41], corresponding to a non-constant source term S in (3)-( 6).An important motivation for our simple shear study is the possibility to later deal with the full X-dependent case.
Under the present assumption of simple shear, the source term S is indeed constant (see below Section 2.3).Interestingly, Equation (9) shows that ḃ(t) ≥ 0 (i.e., the static/flowing interface never drops) when ν = 0 because S(t, b(t)) ≥ 0 and ∂ Z U > 0. This means that for the inviscid model, a layer flowing on top of an initially static layer will not erode, i.e., the static material will not be put into motion (this property is however not true when S depends on Z, in which case we can have ḃ(t) < 0 together with ∂ Z U(t, b(t)) = 0 and S(t, b(t)) = 0, see [41]).For the viscous model (Equation ( 8)), and still with constant source term S, erosion or deposition can occur, depending on whether ∂ 2 ZZ (ν∂ Z U)(t, b(t)) is positive or negative, respectively.

Free Interface Model with Constant Source
The most simple case in the free interface model ( 3)-( 6) is when we take the source term S to be constant in time and space.Indeed, according to ([40] Section 2.3), when all quantities are independent of X, the pressure becomes hydrostatic p = g cos θ(h − Z), and the source takes on a constant value where g > 0 is gravitational acceleration, θ > 0 the constant slope angle (note that this sign convention differs from [40,41]), and µ s = tan δ > 0 the friction coefficient with δ the friction angle related to the material.The source term is thus the result of the balance between the driving force due to gravity (−g sin θ < 0) and the friction force (µ s g cos θ > 0).Indeed according to ([40], Proposition 2.1), for the source term (10), the solution to the system (3)-( 6) is an exact solution to the two-dimensional original viscoplastic model ( 1) and ( 2) without dependency on X.This means that the free interface model ( 3)-( 6) with a constant source term (10) exactly describes the simple shear flows of ( 1) and ( 2), corresponding to the configuration shown in Figure 4. We observe that with (10) and µ s = tan δ, the static equilibrium condition ( 6) is simplified and becomes tan θ ≤ µ s , i.e., θ ≤ δ, (11) meaning that the internal friction must at least neutralize the gravity force due to the slope.This is a necessary condition for the existence of a solution to the free interface model ( 3)- (6).Indeed, if this condition is not satisfied, we expect that "b = 0", meaning that the entire layer immediately flows down the inclined plane.The remainder of the paper is devoted to the evaluation of this constant source term model.

Analytical Solution
For the inviscid model where ν = 0, the free interface model ( 3)-( 6) can be written with the following boundary condition for all t > 0, and with the initial condition ( 5) and the static equilibrium condition (6).Moreover, if the source term is chosen constant and uniform as in Equation ( 10), i.e., with µ s = tan δ, θ ≤ δ, then we can infer an analytical solution.Specifically, the solution to Equations ( 5), ( 12) and ( 13) is given by Equation (15) shows that the velocity profile (at all times when a flowing layer exists) has the same shape as the initial velocity profile: it is just shifted towards decreasing velocities at constant speed S and clipped below the value 0. Note that the velocity profile only depends on g, θ, δ through S in (14).If S = 0, the solution is steady, whereas if S > 0, the velocity decreases with time until the flow stops.Furthermore, the interface position b(t) results from the following implicit equation: This equation has a unique solution in [b 0 , h] for all times t ≤ t stop , where t stop is the time when the whole mass stops, defined by The complete stopping of the flow occurs at t = t stop .For all times t > t stop , the velocity U can be extended by setting U(t, Z) = 0 for all Z ∈ [0, h], and b(t) = h.

Choice of the Parameters and Initial Conditions
In order to compare the analytical solution to the results presented in Figures 2 and 3 extracted from the experiments performed by [12], we have to define the friction angle δ, the slope angle θ, the thickness of the granular layer h, the thickness of the initially static layer b 0 , and the initial velocity profile U 0 (Z).Glass beads of diameter d = 0.7 mm were used, with repose and avalanche angles of about 23 • and 25 • , respectively.Because wall effects are known to increase the effective friction for granular flows in channels such as those of [12], we use here a friction angle of δ = 26 • [27,32,42].We perform different tests by varying the slope angle θ from 19 • to 24 • .Indeed, for slope angles θ much lower than δ, the interface reaches the top of the layer very quickly, preventing us from observing a possible erosion of the static layer.We set b 0 = 5 mm, which is the thickness of the initial static layer in the experiments, and h = 0.02 m, corresponding to the mean thickness at the position X = 90 cm where the measurements were performed, see Figure 3. Thus unless specified, we always take The objective here is only to compare the order of magnitude and the general trend of the analytical solution to the experimental results, since the experiments are more complicated than the uniform granular layer and the initial conditions defined in the model.In particular, the thickness of the granular layer in the experiments may vary by up to 20% during the flow and slightly depends on the inclination angle (see Figure 3).Furthermore, the initial velocity profiles and the maximum velocity also depend on the inclination angle, whereas we impose here the same velocity profile for all the tests.
Velocity profiles in experimental granular flows have been extensively measured in different regimes, see e.g., [25].For free surface flows over rigid inclined beds, the velocity profiles vary with inclination, thickness of the flow, and time.Essentially, the velocity profiles may vary from a linear profile for thin layers over small slope angles to Bagnold-like profiles for higher inclinations.The same general trend is observed for thicker flows (see Figure 5 of [25]).For surface flows over a pile of static grains, the velocity profiles roughly exhibit an upper linear part in the flowing layer and a lower exponential tail near the static/flowing interface (see Figure 6 of [25]).This is consistent with the measurements shown in Figure 3. Furthermore, experimental results suggest that the shear rate ∂ Z U is almost constant and equal to 0.5 g d (see e.g., Equation ( 11) of [25]).As a result, for a linear profile of slope α 1 (see case (a) below), we choose α 1 = 70 s −1 , which is consistent with the velocity profiles measured at t = 0 s (see Figure 3).
In order to investigate the different possible profiles of the velocity, we choose three initial velocity profiles defined, for all Z ∈ [b 0 , h], as: In each case, the maximum velocity is U 0 (h) 1 ms −1 .For each profile, Equation ( 9) provides explicitly the time evolution of the static/flowing interface position as follows: These formulae are valid as long as t ≤ t stop = U 0 (h)/S, i.e., b(t) ≤ h.

Results and Comparison with Experiments
We compare here our analytical results with the experimental data shown in Figures 2 and 3.The velocity profiles at the position X = 90 cm from the gate at different times for slope angles θ = 19 • , 22 • , 23 • , and 24 • were taken from the experiments of [9,12,13].The position of the static/flowing interface has been calculated using a velocity threshold of 1cm/s.The material has been assumed to be static when the particle velocity is lower than this threshold.In order to observe velocity profiles on one side of the flow (through the transparent channel wall) and monitor the evolution of the interface separating flowing and static grains, black beads have been used as tracers at a volume fraction of about 50%.Measurements have been made from successive frames during short time intervals of 0.01-0.04s ( [13]).
The b(t) curves are plotted in Figure 5 and the U(Z) profiles in Figures 6 and 7.The evolution of the static/flowing interface position b(t) obtained from the analytical solution (Figure 5) reproduces to a certain extent the experimental observations.The shape of b(t) is directly related to the velocity profile, as demonstrated by Equation (9).For the analytical solution, depending on the initial velocity profile, the stopping time is in the range 2.75-3 s for θ = 24 • , 1.35-1.5 s for θ = 22 • , and 0.75-0.85s for θ = 19 • , whereas t stop 3.4 s, t stop 1.4 s, and t stop 0.9 s in the experiments, respectively.
As a result, the stopping time is well reproduced by the analytical solution, even though its strong increase for θ = 24 • is underestimated in the analytical solution.On the other hand, the penetration of the static/flowing interface within the initially static bed is not reproduced by the analytical solution that instead predicts a static/flowing interface rising towards the free surface at all times.The decrease of the velocity with time is relatively well reproduced up to about 0.8 s.With the analytical solution, the velocity profiles maintain the same shape whereas the maximum velocity decreases, as roughly observed in the experiments.The decrease of the maximum velocity in the experiments and with the analytical solution are very similar (Figure 7).At later times (t ≥ 1 s) and at θ = 22 • and θ = 24 • , the experiments show a clear change in the velocity profile (see Figure 3c at t = 1 s and t = 2 s) that is not reproduced by the constant shape of the velocity profiles predicted by the analytical solution.Furthermore, the maximum velocity decreases much faster in the experiments.In the experiments, the velocity profiles seem to be closer to linear for smaller slopes (θ = 19 • and θ = 22 • ) and more exponential for θ = 24 • .Referring to Equation (16) and Figure 5, this may explain why b(t) measured experimentally has an exponential shape at θ = 24 • , whereas it is closer to linear for smaller slope angles.

Numerical Solution to the Viscous Free Interface Model
In this section, we consider a numerical solution to the free interface model ( 3)-( 6) with viscosity ν > 0 and constant and uniform source term S given by (14).We first present two robust numerical methods to obtain the numerical solution.Then, we describe the main features of the transient regime by a scale analysis.Finally, we compare our numerical results with experiments, first with constant viscosity and then with variable viscosity deduced from the µ(I) viscoplastic rheology.
Given the unusual form of our model with a free moving interface, some comments are in order.From the mathematical standpoint, it is not obvious that the boundary formulation (4) gives a well-posed problem.Actually, it is known [35] that the µ(I) rheology can lead to an ill-posed problem.We claim that for the simple shear configuration considered here, the model is well-posed (although this would not be the case if X dependency were involved, see [40]).To support this claim, we observe that we consider two different numerical methods which both deliver stable, grid-independent solutions.Moreover, the second numerical method is related to an optimal problem under a very classical form (see (26) below), for which well-posedness is well-known.From a numerical standpoint, we point out that both numerical methods are computationally effective.Let us mention in passing that none of these methods uses the differential Equation ( 8) which would require an intricate treatment of the third-order derivative, but use instead the boundary conditions (4).

Discretization by Moving the Interface
The first numerical method involves rewriting (3), ( 4), ( 5) for a normalized coordinate 0 ≤ Y ≤ 1.We perform the change of coordinates that leads to the differential relations Here, ∂ t and ∂ τ denote the differentiation with respect to time at constant Z and Y, respectively.The change of coordinates (19), and hence the discretization method presented in this subsection, is appropriate as long as there is a flowing layer so that h − b(τ) > 0. Another discretization method dealing with the stopping phase when b(t) reaches the total height h is presented in Section 4.2.The discretization method considered in this subsection has the advantage of tracking explicitly the position of the static/flowing interface, whereas the method of Section 4.2 requires post-processing to evaluate the interface position.
Using the change of coordinates (19), Equation (3) becomes and the boundary conditions (4) become We split the space domain (0, 1) into n Y cells of length ∆Y with n Y ∆Y = 1 and denote the center of the cells by Y j = (j − 1/2)∆Y, for all j = 1 . . .n Y .For n ≥ 0, the discrete times t n are related by t n+1 = t n + ∆t n , where ∆t n is the time step (chosen according to the CFL condition (25) below) and t 0 = 0. We write a finite difference scheme for the discrete unknowns U n j U(t n , Y j ) and b n b(t n ), for all j = 1 . . .n Y and all n ≥ 1, using the initial conditions U 0 (and b 0 ) to initialize the scheme.For all n ≥ 0, given (U n j ) 1≤j≤n Y and b n , the equations to compute (U n+1 j ) 1≤j≤n Y and b n+1 are for all j = 1 . . .n Y , with together with the boundary conditions In Equation ( 22), the values ν n j+1/2 are assumed to be known, and they are computed as interface values from a given variable viscosity ν(t, Z).Equations (24a) and (24c) are used to provide the ghost values U n+1 0 and U n+1 n Y +1 needed in Equation ( 22) for j = 1 and j = n Y , respectively, whereas Equation (24b) is used to determine b n+1 as described below.We observe that in Equation ( 22), the diffusive term is treated implicitly in time, and the first-order derivative of U is treated explicitly using upwinding.As a result, we impose the CFL condition This CFL condition is evaluated approximately using the value ḃn− 1 2 from the previous time step, since ḃn+ 1  2 is unknown at the beginning of the time step; for n = 0, the value 0 is used (hence, no CFL condition is initially enforced, but the time step is taken small enough).Typical values are ∆t 0 = 10 −4 s for the initial time step and ∆Y = 10 −4 .
The solution to Equations ( 22)-( 24) is obtained as follows.We can solve the system ( 22) together with the boundary conditions (24a) and (24c) for any value of b n+1 since it is a linear system in U n+1 with right-hand side containing U n , the source term and the advection term (remember that the advective derivative is treated explicitly).This leads to a tridiagonal linear system with a matrix that results only from the time derivative and the diffusive terms.This matrix, which is diagonally dominant, has an inverse with nonnegative entries.Thus, the solution (U n+1 j ) 1≤j≤n Y can be expressed linearly with nonnegative coefficients in terms of the coefficients (a n+1/2 j ) 1≤j≤n Y which appear on the right-hand side and which depend on the still unknown interface position b n+1 .According to (23) and since U n j is nondecreasing with respect to j (it should remain nondecreasing because we assume ∂ Z S ≤ 0), U n+1 j is thus, for all j = 1 . . .n Y , a nondecreasing continuous and piecewise linear function of b n+1 , with two different formulae corresponding to whether b n+1 is greater or smaller than b n .In particular, the remaining boundary condition (24b), which is equivalent to U n+1 1 = 0 owing to (24a), determines a unique solution b n+1 , see Figure 8.The value of b n+1 can be computed explicitly by solving the linear system twice with two different right-hand sides and using linear interpolation.The first solve uses the right-hand side evaluated with the temporary value b n for b n+1 , yielding a temporary value for U n+1 1 .If the obtained value is negative (left panel of Figure 8), the second solve is performed using the value h for b n+1 , otherwise the value 0 is used (right panel of Figure 8).Once b n+1 is known, we determine the entire profile U n+1 j for all j = 1 . . .n Y by linear interpolation.

Discretization by an Optimality Condition
The second numerical method for solving conditions (3), ( 4) and ( 5) with positive viscosity ν > 0 uses a formulation as an optimal problem set on the whole interval (0, h), valid under the assumption where we consider a no-slip boundary condition at the bottom.Note that this condition becomes relevant whenever the static/flowing interface reaches the bottom, a situation encountered in our simulations as shown in Figure 10c.In (26), the static/flowing interface position b(t) no longer appears explicitly, but has to be deduced from the velocity profile as We split the space domain (0, h) into n Z cells of length ∆Z with n Z ∆Z = h, and denote the center of the cells by Z j = (j − 1/2)∆Z, for all j = 1 . . .n Z .The discrete times t n , for all n ≥ 0, are related by t n+1 = t n + ∆t n , where ∆t n is the time step (chosen according to the CFL condition (30) below) and t 0 = 0. We discretize the problem (26) using a finite difference scheme by writing for the discrete unknowns for all j = 1 . . .n Z .The boundary conditions, which are discretized as U n n Z +1 = U n n Z at the free surface (Z = h) and as U n 0 = −U n 1 at the bottom (Z = 0), are used to provide the ghost values involved in the discretization of the diffusive term.The problem (28) is solved in two steps as Owing to the explicit discretization of the diffusive term, we use the CFL condition We could also consider an implicit discretization to avoid any CFL condition, but each time step would be more computationally demanding.Finally, the thickness of the static layer can be evaluated as However, to prevent b n from being influenced by small values of U and for better accuracy, we prefer to use where C 0 is an appropriate constant of the order of S/ν, in relation to Equations (4a), (4b) and (7).

Scales in the Transient Regime
We assume that the source term S is constant, and that the viscosity ν is also constant.Then, the solution to the viscous model ( 3)-( 5) depends on the constant S, the total thickness h, the viscosity ν, the initial thickness of the static layer b 0 , and the initial velocity profile U 0 (Z).We introduce dimensionless quantities, denoted by carets (hats), where τ is a time scale, l a space scale, and u = lα 1 with α 1 being the order of magnitude of the initial shear rate.In order to write Equation (3) in dimensionless form, we take The dimensionless equation is then According to the above scaling analysis, if the velocity profile is initially linear with shear α 1 , and if then the dimensionless problem has no parameter left and the velocity U is therefore a fixed profile (in terms of Z − b 0 ), as well as b − b 0 .We conclude with (33) that the quantities t c , b 0 − b min , ḃ∞ are proportional to τ, l, l/τ, respectively.We obtain the following proportionality factors from the numerical simulation Note the dependency on the ratio α 1 /S, which is due to the homogeneity of the problem ( 3) and ( 4) with respect to (U, S) (a new solution can be obtained by multiplying (U, S) by a positive constant, with b unmodified).

Results and Comparison with Experiments for Constant Viscosity
We consider a constant and uniform source term of the form (14).The case of a linear initial velocity profile is simulated using different values of the viscosity ν (taken constant) and slope angles θ.The results obtained with the two methods of Sections 4.1 and 4.2 are always identical, thus we shall not specify which one is used.As discussed in Section 4.3, the static/flowing interface always penetrates the initially static layer, in contrast to what has been observed for the inviscid model, within a length b 0 − b min proportional to ν, according to (36).In other words, the flow excavates the static bed.However, for ν < 10 −5 m 2 s −1 , the penetrating length is too small to be observed.For ν = 10 −5 m 2 s −1 (Figure 10a), b 0 − b min = 8 × 10 −4 m for θ = 24 • (with t c = 5 × 10 −2 s) and b 0 − b min = 2 × 10 −4 m for θ = 19 • (with t c = 4 × 10 −3 s).Thus, the static/flowing interface penetrates only slightly within the initially static layer.As the viscosity increases (Figure 10b), the static/flowing interface penetrates deeper into the initially static layer and even reaches the bottom for ν = 10 −4 m 2 s −1 at θ = 24 • (Figure 10c).The results that best reproduce the experimental observation for the penetration of b(t) within the initially static layer are obtained with ν 5 × 10 −5 m 2 s −1 .In good agreement with the experiments, the simulations with viscosity predict that the depth and duration t c of penetration of the static/flowing interface within the initially static layer increase with the slope angle.Furthermore, the values of b(t) and t c are in reasonable agreement with those observed experimentally (Figure 10b and Figure 17 of [12]).
Qualitatively similar results are obtained using different initial velocity profiles (Figure 12).However, the shape of b(t) is affected by the choice of the initial velocity profile.For example, for an exponential initial velocity profile with θ = 24 • and ν = 5 × 10 −5 m 2 s −1 , the static/flowing interface position b(t) stagnates at an almost constant position for the first 0.5s, contrary to the case of a linear initial velocity profile (Figures 10b and 12b).
The convex shape of b(t) during the migration of the static/flowing interface up to the free surface obtained with the viscous model is very different from that observed and from that obtained with the inviscid model, which predicted a linear shape related to the linear initial velocity profile.With the viscous model, the time evolution of b(t) and the velocity profile U(t, Z) are not so clearly related to the shape of the initial velocity profile, as shown for example in Figure 12 for θ = 24 • .The velocity profiles for the viscous model are also very different from those for the inviscid model.Whatever the shape of the initial velocity profile, the velocity profiles later exhibit an exponential-like tail near the static/flowing interface, similar to that observed experimentally.They also exhibit a convex shape near the free surface which can be observed, although not as marked, in some but not all experimental velocity profiles (e.g., Figure 3c).While the maximum velocity decreases too fast at t = 0.5 s and θ = 24 • compared to the experiments and to the inviscid model, the velocity is much closer to the experiments at t = 1 s than with the inviscid model (Figure 11b).For θ = 19 • and θ = 22 • , the decrease in velocity at time t = 0.7 s is overestimated in the viscous model.
Finally, we briefly look at the stopping time t stop of the whole granular layer.The stopping time for the inviscid model and for the experiments is given in Table 1, for different slope angles.The stopping time is smaller for viscous than for inviscid flow, whatever the shape of the initial velocity profile.Given that the stopping time for the inviscid model is expressed as t stop ν=0 = U 0 (h)/S, we can evaluate numerically the difference (t stop ν=0 − t stop ν ), where t stop ν denotes the stopping time for the viscous model.Our numerical simulations show that this difference depends to a moderate extent on the viscosity (Figure 13 where the slope of the curves suggests a behavior of this difference close to ν 1/4 for the present parameters).In any case, the presence of viscosity diminishes the stopping time.Consequently, according to Table 1, the stopping time is significantly smaller in the viscous model than in the experiments.This can also be observed in Figures 10 and 12     with respect to viscosity in log scale.

Variable Viscosity
The viscoplastic description of granular materials involves fundamentally a Drucker-Prager yield stress proportional to the pressure, the coefficient being the static friction parameter µ s = tan δ.A full rheological law including this yield stress (i.e., defined for all values of the strain rate and not only close to zero) has been proposed in [28], the so called µ(I) rheology.As described in [32], this rheology can be interpreted as a decomposition (i.e., (2)) of the deviatoric stress tensor in a rate-independent pure plastic part proportional to µ s and a viscous part with pressure-and rate-dependent dynamic viscosity η given by where p dyn is the dynamic pressure, D is the strain rate tensor, , and I is the inertial number defined by with d the grain diameter as before and ρ s the grain density.The kinematic pressure p is related to p dyn by p = p dyn /ρ, with ρ = φρ s the density of the granular material, φ being the volume fraction.
If we consider a slope aligned velocity field depending only on the normal coordinate Z (simple shear flow), the arguments of Section 2.3 show that the solution to the viscoplastic model ( 1) and ( 2) is described by the system (3)-( 5) with the constant and uniform source term S given by ( 14) and ν = η/ρ.Such a flow has hydrostatic pressure p = g cos θ(h − Z) and shear rate D = ∂ Z U/2.Thus with (37), the kinematic viscosity becomes with Then the term appearing in (3) is and the part of this term proportional to µ s indeed balances in (3) the corresponding part in S from ( 14).
The nonlinearity is given according to [28] as with µ 2 > µ s the friction at large strain rate and I 0 is a constant that depends on the material used and other device specificities, see ( [27], Appendix A).Here we take the value I 0 = 0.279 proposed by [27].Note that the boundary condition (4c) is automatically satisfied and can therefore be skipped.
For the discretization, we use the method of Section 4.2 with the discrete unknowns U n j U(t n , Z j ), Z j = (j − 1/2)∆Z for j = 1 . . .n Z , n Z ∆Z = h.The Equations ( 3) and ( 41) with ( 42) and ( 40) are discretized with finite differences under the conservative form for all j = 1 . . .n Z , with for j = 0 . . .n Z , with Z j+1/2 = Z j + ∆Z/2 = j∆Z, and once U n+1/2 j has been computed, we set Since Φ n n Z +1/2 = 0, there is no need to define U n n Z +1 , in accordance with the loss of the free surface boundary condition (4c).On the left boundary (bottom), we use as before the no-slip condition, We use the CFL condition 2ν max ∆t n ≤ ∆Z 2 , with ν max defined as the maximum value of ν computed with (39).The interface position b n is computed according to (32).
As previously, we take µ s = tan(26 • ), h = 0.02 m, b 0 = 0.005 m, and d = 7 × 10 −4 m.As in [12], we take φ = 0.62.We choose µ 2 = tan(28 • ), which leads to ν(Z = b 0 ) 5 × 10 −5 m 2 s −1 , corresponding to the order of magnitude of ν taken in Section 4.4.We consider the case of a linear initial velocity profile, with slope angles θ = 19 • , 22 • , 24 • .The numerical results for the static/flowing interface position b(t) and velocity profiles U(Z) are plotted in Figures 14 and 15.In comparison with numerical results of Figures 10b and 11a,b obtained with the constant viscosity model, the shapes are similar, in particular for the initial erosion of the static bed.The behavior of the interface position b(t) reaching the free surface close to the stopping time is here closer to the prediction of the inviscid model (see Figure 5a).However, this behavior does not recover (since ν vanishes at the free surface) the expected behavior, which is roughly linear, of the experimental results of Figure 2. The comparison with the experimental results is better with the µ(I) rheology than with a constant viscosity.For clarity, the static/flowing interfaces and the velocity profiles corresponding to the different cases are plotted in Figure 16  It is important to explain the behavior of the solution close to the free surface.As said above, there is no need to impose any boundary condition at the free surface, because the pressure vanishes there.Nevertheless, the behavior close to the free surface of the solution to the parabolic problem (3) is deduced from the property that ∂ Z (ν∂ Z U) is bounded, which gives ν∂ Z U ∼ h − Z.For the model with constant viscosity, this leads to ∂ Z U ∼ h − Z and U has a (rather) parabolic profile close to the free surface.For the model with µ(I) rheology, Equation (41) implies that I tends to a finite value as Z → h, and thus with Equation (40), we obtain ∂ Z U ∼ √ h − Z, and U has a (rather) Bagnold profile close to the free surface.We conclude that in any case, the Neumann condition ∂ Z U = 0 is recovered at the free surface, even if not enforced explicitly.For the model with the µ(I) rheology, ∂ Z U tends to zero more slowly than for the model with constant viscosity, leading to a behavior closer to the prediction of the inviscid model.Nevertheless, the property ∂ Z U = 0 at the free surface, which seems to be a consequence of the incompressible viscoplastic model, makes it impossible to obtain a good representation of the experimental stopping phase.Note that the µ(I) rheology was introduced precisely to represent the Bagnold profile (without a static phase) U(Z) = c(h 3/2 − (h − Z) 3/2 ), where c is determined by the relation µ(I) = tan θ, which gives This Bagnold profile without a static phase is a steady solution to the system (3) and ( 41) with ( 42) and ( 40), but only when µ s < tan θ < µ 2 .In our framework, we have a static phase, the flow is not steady, and tan θ < µ s .The formula (46) is therefore not applicable and would give a negative c.Our choice of α 3 is nevertheless of the same order of magnitude as formula (46).

Discussion and Conclusions
We have considered the 1D (in the direction normal to the flow) non-averaged model with static/flowing dynamics that was initially proposed in [41].This model is based on the description of a granular material by a yield stress viscoplastic rheology.Here we have assumed that the source term in this model is constant, which corresponds to simple shear solutions to the viscoplastic model.Thus, all the quantities are assumed to be independent of the down-slope variable and in particular the flow thickness is constant.We have shown that this free interface model can be solved with reliable and relatively simple numerical methods.We have compared model solutions for both the inviscid and viscous models to observations from experiments on granular flow over an inclined static layer of grains.The analytical solution for the inviscid model and the numerical results for the viscous model reproduce quantitatively some essential features of the change with time of the velocity, of the static/flowing interface position, and of the stopping time of the granular mass, even though the flow thickness in the experiments is not perfectly uniform and the initial velocity profile changes with slope angle and flow thickness as discussed below.
The analytical solution for the inviscid model shows that the evolution of the static/flowing interface position is proportional to the source term and inversely proportional to the shear rate (Equation ( 9)).For the viscous model (with constant viscosity), the analysis shows that the evolution of the interface is related to the viscosity, the source term, and the first-and third-order derivative of the source term and the velocity, respectively (Equation ( 8)).Owing to the appearance of this third-order derivative, the dynamics of the static/flowing interface cannot be reduced to a simple differential equation in terms of depth-averaged quantities.While the shape of the initial velocity profile is preserved at all times for the inviscid model according to Equation ( 15), the viscous model predicts that an exponential-like tail near the static/flowing transition and a convex shape near the free surface develop.The viscous contribution enables the static/flowing interface to initially penetrate within the static layer (which is eroded), as observed in the experiments, and contrary to the inviscid model.
The viscoplastic model used here has the great advantage of involving only two parameters, i.e., the friction coefficient µ s and the viscosity ν (or the coefficient µ 2 for the µ(I) rheology), whereas the so-called partial fluidization theory, involving an order parameter to describe the transition between static and flowing material, also reproduces the erosion of the static bed and the velocity profiles predicted here with the viscous model [7], but at the cost of additional empirical equations for the time-change of a state parameter.
One of the important results of the analysis lies in the explicit expressions obtained from the analytical solution, especially for the time evolution of the static/flowing interface.The dynamics are controlled by the source term S that is constant when the pressure is assumed to be hydrostatic and when the slope and thickness are assumed to be constant.In such a case, S(t, Z) = S = g cos θ(tan δ − tan θ), with µ s = tan δ.Comparable results have been obtained from the analytical solution of thin-layer depth-averaged equations for a granular dam-break (i.e., with non-constant h) [43,44].This analytical solution predicts a granular mass front velocity that decreases linearly with time where k 0.5 (see [9]).Furthermore, the stopping time of the analytical front for the granular dam-break is Equations ( 47) and (48) for the front velocity and the stopping time of the front for a depth-averaged model of a granular dam-break are very similar to Equations ( 15) and ( 17) respectively, found here for the velocity of the flow and for the stopping time of the granular layer in the non-averaged case and without viscosity.
Although the initial penetration of the static/flowing interface into the static layer (erosion) can be reproduced by taking into account the viscosity, this leads at the same time to an underestimation of the stopping time.The viscous model better reproduces the exponential-like tail of the velocity profile near the static/flowing interface than the inviscid model, but overestimates the convexity of the velocity profile near the free surface.Furthermore, the change in shape of the velocity profile observed in the experiments during the stopping phase is reproduced by the viscous model, as opposed to the inviscid model.However, the decrease in maximum velocity near the surface is too fast in the viscous model.All these results suggest that (i) viscosity plays an important role near the static/flowing interface at depth and in this region a reasonable estimate for the viscosity is ν 5 × 10 −5 m 2 s −1 and (ii) viscous effects in the experiments seem to be much smaller near the free surface.These observations suggest a non-constant viscosity, as proposed in the so-called µ(I) flow law, e.g., [25][26][27][28].As in [32], we have used the µ(I) flow law to derive the value of the viscosity, which becomes variable and nonlinear with respect to the shear rate.This flow law reproduces approximately the value of the viscosity at the static/flowing interface (i.e., ν 5 × 10 −5 m 2 s −1 ) and at the free surface (vanishing viscosity).However, this model does not allow us to significantly improve the accuracy of the behavior of the static/flowing interface close to stopping.This is because it imposes a Bagnold behavior close to the free surface, which is in contradiction with experiments that predict a linear behavior.A more detailed numerical simulation of the granular collapse over sloping beds with the complete viscoplastic model and the µ(I) rheology gives a dynamic viscosity of around 0.2Pa• s near the front, increasing up to about 0.4Pa•s behind the front (see Figure 11 of [33] and Figure 10 of [32]), leading to a kinematic viscosity of ν = 1.3 × 10 −4 m 2 s −1 and ν = 2.6 × 10 −4 m 2 s −1 , respectively, for a density of ρ = 1550kg/m 3 .Note also that our treatment of side wall effects deserves further improvements, since we simply increased the friction coefficient as done in [27,32,42].Taking into account wall effects changes the static/flowing position, as shown in particular in [33] for the granular collapse considered here.More precisely, the static/flowing interface is found to be less deep within the flow when wall effects are accounted for (see Figures 3 and 4 of [33]).This effect could possibly increase the value of the viscosity needed in the model to reproduce the penetration of the static/flowing interface observed in the experiments.This would give values of the viscosity closer to that given by the µ(I) rheology near the front (Figure 11 of [33]).
It would be of interest to take into account the X-variations of the source S from ( [41], Equation (2.12)), thereby accounting for topography, propagation and non-hydrostatic effects.This would make it possible to study erosion in flows over complex topography in the laboratory (e.g., [45]) or at the natural scale (e.g., [46]).For the very simple case of flows over a constant slope and if the pressure is assumed to be hydrostatic (p = g cos θ(h − Z)), taking into account the X-variations is the same as (according to [41]) considering once again (3)-( 6) with S including an additional term, and as considering a thickness h = h(t, X) that evolves according to (2.2) in [41].Then U = U(t, X, Z) and b = b(t, X).However S = S(t, X) in (49) is still independent of Z. Then Equations ( 8) and ( 9) read as follows with S(t, X) defined by (49): • If ν > 0 and S(t, X) = 0, then ḃ(t, X) = −∂ 2 ZZ (ν∂ Z U)(t, X, b(t)) S(t, X) ν. (50) • If ν = 0 and ∂ Z U(t, X, b(t)) = 0, then ḃ(t, X) = S(t, X) The model with 0 < b < h is valid as long as S(t, X) ≥ 0, so as to satisfy (6).Indeed, violating this condition leads to instantaneous flow of the whole layer i.e., b = 0, and full erosion of the bed.In this case, (3) still has to be solved (with b = 0 and S < 0), but with boundary conditions (4a) and (4c) only.If at some later time S becomes positive again, b can rise above the value b = 0. Thus when S(t, X) can change sign, there is in general a succession of progressive depositions and sudden full erosions of the bed.As discussed previously, the experimental granular flow is not perfectly uniform, as illustrated for example for the granular collapse at θ = 22 • in Figure 17.In this case, tan δ − tan θ = 0.0837, so that near the front, where the h-gradient is strongest (∂ X h(t = 0.66 s, X = 90 cm) 0.125 as shown in Figure 17), S is negative and thus even the inviscid model predicts that erosion occurs.Later, when the front flows beyond the position X = 90 cm, ∂ X h decreases with values around 5 × 10 −3 .Therefore S becomes positive, and with the inviscid model, the static/flowing interface rises (i.e., deposition occurs).For the viscous model, for the case S ≥ 0 studied in this paper, occurrence of erosion ( ḃ < 0) or deposition ( ḃ > 0) is related according to (50) to the change of sign of ∂ 2 ZZ (ν∂ Z U)(t, X, b(t)), i.e., positive and negative, respectively.Therefore, when the X variations are considered and under the hydrostatic assumption, the occurrence of erosion or deposition depends on the signs of both S and ∂ 2 ZZ (ν∂ Z U)(t, X, b(t)).The influence of a Z-dependency of S on the static/flowing interface dynamics and on the erosion process in the free interface model ( 3)-( 6) is studied in [41].It would be interesting to extend the approach proposed here to 2D and possibly 3D, so as to capture the static/flowing interface in thin-layer models, as proposed in [40].The orders of magnitude assumed in [40] are indeed satisfied in the experiments discussed here, because the typical length is L = 1 m, the typical time is τ = 0.33 s, (so that L/τ 2 = g), h = 0.02 m, and ν = 5 × 10 −5 m 2 s −1 , leading to ε ≡ h/L = 0.02, tan δ − tan θ = O(ε), and the normalized viscosity ντ/L 2 10 −5 is of the order of ε 2 or ε 3 .The primary models of [40] with dependency on X and without hydrostatic assumption (giving S depending also on Z) lead however to severe nonlinearities and ill-posedness, so that the numerical treatment of these extensions with all flow-aligned variations represents a major challenge.
An important issue is also to summarize the dynamics of the normal velocity profile using a finite number of parameters (for example the interface position b, the thickness h, and the shear rate), in order to keep computational costs low enough to simulate natural situations.This could lead to a depth-averaged model, which has until now seemed inaccessible because of the dependency on the third normal derivative in the differential Equation (8).[16,17], at times t = 0.66 s, 0.78 s, 1.02 s.The value X = 0 corresponds to the position of the gate.The slope angle is θ = 22 • .Note that the vertical scale differs from the horizontal one, since the thickness of the flow represents only 7% of its horizontal extension.

Figure 1 .
Figure 1.Experimental setup from[9,12] to study granular column collapse over inclined planes covered by an initially static layer made of the same grains as those released in the column.For slope angles θ a few degrees smaller than the typical friction angle δ of the involved material, a quasi-uniform flow develops behind the front.

Figure 2 .
Figure 2. Position of the static/flowing interface b as a function of time t until the granular mass stops, measured at X = 90 cm from the gate, from experiments of granular collapse over an initially static granular layer of thickness b 0 = 5 mm on an inclined channel of slope angle θ = 19 • (green squares), θ = 22 • (blue stars), θ = 23 • (red crosses), and θ = 24• (black vertical crosses).Time t = 0 s corresponds to the time when the front of the flowing layer reaches the position X = 90 cm.The approximate upward velocity ḃ of the static/flowing interface is indicated in m/s for each slope angle.These new results have been extracted from the experiments performed by[12,13] for granular columns of initial down-slope length r 0 = 20 cm, initial thickness h 0 = 14 cm, and width W = 20 cm (i.e., volume V = 5600 cm3 ).Note that the position of the free surface when the mass stops, i.e., when b = h, represented by the upper point for each angle, slightly depends on the slope angle and decreases as the slope angle increases, from about 0.025 m at θ = 19 • to about 0.018 m at θ = 24 • .

Figure 3 .
Figure 3. Velocity profiles U(Z) at different times until the granular mass stops, measured at X = 90 cm from the gate, in experiments of granular collapse over an initially static granular layer of thickness b 0 = 5 mm on an inclined channel of slope angle (a) θ = 19 • ; (b) θ = 22 • ; and (c) θ = 24 •.Time t = 0 s corresponds to the time when the front of the flowing layer reaches the position X = 90 cm.These new results have been extracted from the experiments performed by[12,13] for granular columns of initial down-slope length r 0 = 20 cm, initial thickness h 0 = 14 cm, and width W = 20 cm (i.e., volume V = 5600 cm 3 ).

Figure 4 .
Figure 4. Simplified flow configuration consisting of a uniform flowing layer over a uniform static layer, both parallel to the rigid bed of slope angle θ.The static/flowing interface position is b(t) and the total thickness of the mass h is constant.

Figure 5 .Figure 6 .
Figure 5. Static/flowing interface position b as a function of time t for the inviscid model, different slope angles and an initially static granular layer of thickness b 0 = 5 mm, using (a) a linear, (b) an exponential, and (c) a Bagnold initial velocity profile.The experimental results from Figure 2 are shown for comparison.

Figure 7 .
Figure 7. Velocity profiles U(Z) at different times for the inviscid model, with an initially static granular layer of thickness b 0 = 5 mm over an inclined plane of slope (a,c) θ = 19 • and (b,d) θ = 24 • , using (a,b) an exponential and (c,d) a Bagnold initial velocity profile.The experimental results from Figure 3 are shown for comparison.

Figure 8 .
Figure 8. Velocity U n+1 1 versus b n+1 .The chosen value for b n+1 is determined by the intersection of the curve with the horizontal axis.

Figure 9 .
Figure 9. Schematic evolution of the thickness of the static/flowing interface as a function of time. (right).

Figure 14 .Figure 15 .Figure 16 .
Figure 14.Static/flowing interface position b as a function of time t, for variable viscosity associated with the µ(I) law, with linear initial velocity profile and slope angles θ = 19 • , θ = 22 • , θ = 24 • .The experimental results from Figure 2 are shown for comparison.

Figure 17 .
Figure 17.Thickness profile h as a function of down-slope position X in the experiments performed by[16,17], at times t = 0.66 s, 0.78 s, 1.02 s.The value X = 0 corresponds to the position of the gate.The slope angle is θ = 22 • .Note that the vertical scale differs from the horizontal one, since the thickness of the flow represents only 7% of its horizontal extension.