Explicit solutions to a free interface model for the static/flowing transition in thin granular flows

This work is devoted to an analytical description of the dynamics of the static/ﬂowing interface in thin dry granular ﬂows. Our starting point is the asymptotic model derived by Bouchut et al. (2016) from a free surface incompressible model with viscoplastic rheology including a Drucker–Prager yield stress. This asymptotic model is based on the thin-layer approximation (the ﬂow is thin in the direction normal to the topography compared to its down-slope extension), but the equations are not depth-averaged. In addition to the velocity, the model includes a free surface at the top of the ﬂow and a free time-dependent static/ﬂowing interface at the bottom. In the present work, we simplify this asymptotic model by decoupling the space coordinates, and keeping only the dependence on time and on the normal space coordinate Z . We introduce a time-and Z -dependent source term, assumed here to be given, which represents the opposite of the net force acting on the ﬂowing material, including gravity, pressure gradient, and internal friction. We prove several properties of the resulting simpliﬁed model that has a time- and Z -dependent velocity and a time-dependent static/ﬂowing interface as unknowns. The crucial advantage of this simpliﬁed model is that it can provide explicit solutions in the inviscid case, for diﬀerent shapes of the source term. These explicit inviscid solutions exhibit a rich behaviour and qualitatively reproduce some physical features observed in granular ﬂows.


Introduction
Dense granular materials can behave like a solid or flow like a fluid. The static (solid) and flowing (fluid) zones are delimited by a static/flowing interface that changes with time, depending on the initial and boundary conditions and on the driving and resistive forces within the granular mass. This static/flowing transition is characteristic of granular flows and plays a key role in both industrial and geophysical contexts, for instance in erosion/deposition processes [20,23,39]. Dry granular flows can be described by a viscoplastic rheology with yield stress. With this rheology, the granular material starts to flow if the stress is above a threshold value and otherwise remains at rest. A Drucker-Prager yield stress proportional to the pressure has been shown to describe granular flows well [18,20,27,31,33,40]. The proportionality coefficient µ can be constant or not, as in the so-called µ(I) rheology [31]. In particular, this relationship makes it possible to obtain a flowing layer near the free surface, above a static layer near the bottom, in contrast with the Bingham flow law involving a constant yield stress and leading to a plug zone located on top of a sheared bottom layer [6]. The time evolution of the static/flowing interface (also called the yield surface) is implicitly described in these models. Up to recently, no explicit description of the movement of this interface was found, owing to the strong nonlinear coupling of the different terms involved in the equations. This is, however, a key issue for practical studies.
The aim of this work is to describe the evolution of the static/flowing interface in dry granular flows over inclined topography, described by the incompressible viscoplastic model with Drucker-Prager yield stress. To simplify the problem, we focus here on granular flows above a static layer made of the same grains such as heap flows (e.g. Figure 6a,b of [41]) or surface flows occuring during unsteady flows (e.g. Figure 18, t = 0.5s of [27]). In such cases, the velocity decreases almost exponentially down to zero from the flowing region to the static part. We make the thin-layer approximation, i.e. we assume that the flow thickness is small compared to its down-slope extension. This assumption is widely used in fluid dynamics and geophysics. However, as described in [15], a key point of our approach is that, contrary to usual thin-layer (i.e. shallow-water) models, the conservation equations are not depth-averaged. The flow here can depend on the normal coordinate to the topography and this dependence can vary with time. Indeed, because of the depth-averaging process, the static/flowing interface in the direction normal to the topography is often neglected in thin-layer models [13,14,22,26,36,38], i.e. a column of granular material is assumed to be either fully static or fully flowing. Here our starting point is the thin-layer, non-averaged, asymptotic model derived recently in [15] from a viscoplastic model with Drucker-Prager yield stress. In this model, the unknowns are a velocity (that can depend on time, down-slope and normal space variables X, Z) and an upper free surface and a static/flowing interface that both depend on time and on X. The change of momentum is driven by a source term S (opposite of the net force acting on the flowing material) that takes into account the effects of gravity, hydrostatic pressure gradient, and pressure-dependent friction. A positive S has a decelerating effect on the velocity, while a negative S has an accelerating effect. The source term S includes the effects of down-slope space inhomogeneities via the hydrostatic pressure gradient and the nonhydrostatic pressure factor in the friction term (see equation (2.14)). The down-slope space inhomogeneities are supposed to give rise to interesting dynamics, see for example [22]. However, the asymptotic model from [15] is still too elaborate to perform an analytical study directly.
In the present work, we devise a simpler form of the asymptotic model from [15] that allows us to perform an analytic study while keeping some of the essential physical processes. Our approach is based on decoupling the down-slope and normal coordinates. Our key assumption is to consider that the down-slope coordinate X is fixed and we therefore look at the model expressed in terms of the normal variable Z only. Doing this, we can no longer couple the part of the source term coming from the nonhydrostatic pressure because it involves the X derivative of the velocity. We therefore have to consider (at fixed X) a given source term S(t, Z), leading to a model in which the unknowns are a velocity u(t, Z) and a static/flowing interface position b(t) (measured in the direction normal to the topography). The main contribution of our paper is to show that the Z-dependence of this source term and its shape induce particularly rich dynamics on the velocity field and on the interface that reflects various physically relevant behaviours as further discussed below. Considering a source term S that is nonconstant in Z allows us to investigate the possible effects of space inhomogeneities and of nonhydrostatic pressure on the flow dynamics. Our simplified model takes the form of a free-interface problem for the static/flowing interface. We emphasize that the dynamical equations are set only in the flowing domain, i.e. above the position of this interface. This free-interface problem is formulated mathematically as a well-posed overdetermined boundary-value problem, with one additional condition with respect to a problem with a fixed boundary. This additional condition governs the time evolution of the static/flowing interface. We compute explicit solutions to the simplified model in the inviscid case, i.e. time-dependent velocity profiles and the time evolution of the interface position. We show that the time and space variations of the source term lead to a particularly rich structure that is able to capture several physically relevant effects from the starting to the arrest of the flow, including progressive starting, progressive stopping and a sudden start of a part of the granular mass. These are, in our opinion, important effects that come into play also when space inhomogeneities and nonhydrostatic terms are taken into account as in [15]. Further results on the present simplified model for a constant source term can be found in [35], where numerical solutions are obtained in the viscous case and where it is shown that important features observed in laboratory experiments of granular flow on erodible beds [23] can be reproduced.
Let us finally put the present work in perspective with the literature. Several attempts have been made to define an equation for the static/flowing interface in thin-layer granular models (see [12,29,48] and references therein). These works describe the interface evolution directly, without referring to a precise viscoplastic model. Most of them are based on prescribed velocity profiles derived from ad hoc considerations, laboratory measurements or solutions of models under steady uniform flow conditions (e.g. the linear velocity profile in [3,21,32] or the S-shaped velocity profile in [17]). A steady-state velocity profile is also assumed in [4] to derive the static/flowing dynamics from a phase transition fluidization model. However, laboratory measurements and numerical simulations indicate that the velocity profile changes with time and with the down-slope space variable, in particular during the starting or stopping phase of the flow of a granular mass [20,24,27,42,8]. Since the dynamics of the static/flowing transition are strongly related to this velocity profile (see [15] or equations (3.16) and (3.17) below)), imposing a velocity profile is only valid under very specific flow conditions. Other models in the literature deduce the position of the static/flowing interface from ad hoc phenomenological erosion laws (i.e. the exchange rate between the flowing layer and the static layer) [43] that generally appear only in the mass conservation equation [11,16]. However, the effect of this mass exchange should also be accounted for in the momentum conservation equation [12,25,29]. Furthermore, defining this exchange rate independently on the basis of phenomenological considerations can lead to models that do not have a correct energy balance [15]. Other methods oversimplify the problem by reducing the flow to that of a sliding block [28]. Finally, an equation for the static/flowing interface dynamics can be found in simple flows by defining jump conditions between the static and the flowing layers [29]. However, this equation involves the value of the velocity at the interface, which is not known in depth-averaged models without imposing a velocity profile or adding additional assumptions. In our approach, the velocity simply vanishes at the interface, while the dynamics of the interface and of the velocity are intricately coupled. More precisely, the evolution of the interface involves normal derivatives of the velocity up to the third order (see equation (3.16)), and the velocity obeys a momentum equation (3.1) with boundary conditions (3.3) on a moving domain. An important feature of our approach is that both the interface and the velocity are solved simultaneously. Our model is established under some assumptions concerning the flow regime, but it has no empirical input. Another fully coupled approach is proposed in [42]. In that work the µ(I)-rheology is resolved nonlinearly in the flowing region and with particular initial data. However the analysis of the stress is not done in the static region, and as a consequence the interface condition (3.3b) is not satisfied. This condition is necessary because the shear stress Σ XZ is continuous across the interface, see Section 2 and (2.8b). Numerical simulations of the full problem show however that the solution of [42] is not too far from the exact one. This solution of [42] can also be seen to be quite close to the solution to a partially regularised model of the incompressible µ(I)-rheology where a creep state at low inertial numbers is introduced and no precise interface location is needed [8].
The outline of this paper is as follows. In Section 2, we briefly describe the non-averaged, thin-layer asymptotic model of [15] and review some of its essential features. In Section 3, we derive our simplified model with a given source term depending on time t and on the normal space coordinate Z. In addition to the velocity, it includes the time-dependent position b(t) of the static/flowing interface as an unknown of the problem. Several properties of our simplified model are derived. In Section 4, we obtain explicit solutions to our simplified model in the inviscid case. Throughout the section, we assume that the source term has a (time-dependent) zero in space denoted by b ⋆ (t), i.e. S(t, b ⋆ (t)) = 0, so that the position b ⋆ (t) splits the flow in two layers where the net force is either driving or resistive. We identify four different dynamical behaviours depending on the time evolution of b ⋆ (t) and on the initial data, from the starting to the arrest of the flow, including progressive starting, progressive stopping and a sudden start of part of the granular mass. These scenarios are expected to appear in the general asymptotic model when the source term is nonlinearly coupled to the velocity. In each case, we determine an explicit solution (normal velocity profile u and interface position b as a function of time), we establish some properties of the solution and we graphically illustrate our results. Finally, in Section 5, we conclude the paper by providing an overview of all the dynamical behaviours exhibited in Section 4.

Non-averaged thin-layer asymptotic model
A non-averaged thin-layer model for viscoplastic flows over an inclined rigid bed has been derived in [15]. It makes it possible to describe the evolution of the static/flowing interface in granular flows. It is formulated in the coordinates X in the direction tangent to and Z normal to the topography (see Figure 1 showing the case of flat topography). The thin-layer model is derived from an incompressible (div u = 0) viscoplastic rheology with Drucker-Prager yield stress [31] where σ is 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 with u the velocity vector. Here the norm of a matrix The coefficients µ s > 0 and ν ≥ 0 are the internal friction and the kinematic viscosity, respectively. The viscosity ν can be time and space dependent, in order to include the case of µ(I) rheology, as described in [27]. In this case the viscosity is given by 2ν = (µ(I) − µ s )p/ Du , so that 2νDu + µ s p Du Du = µ(I)p Du Du . In the rheological model (2.1), we recall that the ratio Du/ Du is multivalued, meaning that it can take on an arbitrary (time and space dependent) trace-free symmetric matrix value of norm less than or equal to unity whenever Du = 0. It is known that the viscoplastic model (2.1) can be linearly well-posed or ill-posed [9,40], depending on the data (viscosity, strain rate and pressure values). Note that for the inviscid model ν = 0, the µ(I) model defined above, as well as for the constant viscosity model ν = cst > 0, the system is ill-posed for small strain rate Du, and in particular close to the static/flowing interface, see respectively [45,9,40]. At the time being there is no incompressible model which is always well-posed and with comparable relevance. Some compressible models have been proposed recently [10,46].
Assuming that the thickness of the layer, the curvature of the topography, the viscosity, and the velocity are small, that the internal friction angle is close to the slope angle and that the pressure is convex with respect to the normal coordinate Z, a formal expansion [15] of the governing equations leads to the following thin-layer model describing a flow including a static layer. The thickness h(t, X) evolves classically according to mass conservation as where u(t, X, Z) is the velocity in the direction tangent to the topography and b(t, X) is the position of the interface between the static part Z < b(t, X) (where we set u ≡ 0) and the flowing part Z > b(t, X). Considering a non-zero slipping velocity for the solid part (which means that u would be constant for Z < b(t, X)) would be of interest also, such a situation can be seen on [7, Fig. 4], but this is out of the scope of the present paper. The momentum balance equation writes where g > 0 is the gravitational acceleration, θ(X) < 0 the slope angle of the topography, and Σ XX , Σ XZ the components of the stress σ in (2.1) in the inclined frame, expressed as where p(t, X, Z) is the pressure, coupled to the velocity u by a nonlinear relation (see (2.13) below). The inertial terms have been neglected in (2.3) and leading-order approximations have been invoked in (2.4). In particular, Σ XX retains only the (opposite of the) hydrostatic part of the pressure p, whereas the whole pressure p is kept in Σ XZ because of the Z-derivative in (2.3) that makes this term dominant in view of the variations at the scale of the thin layer. In (2.4), the term sgn(∂ Z u) is multivalued, meaning that this term takes on an arbitrary (time and space dependent) value between −1 and 1 whenever ∂ Z u = 0. This reflects the yield stress rheology of (2.1). The system is completed by the no-stress condition at the free surface Σ XZ = 0 at Z = h.
where the source term S, which depends on t, X, Z, is the opposite of the net force (excluding the viscous force) in the flowing layer. Assuming a negative angle θ < 0 and an increasing velocity profile in the flowing layer, i.e. ∂ Z u > 0 for all Z ∈ (b, h), the source term S is expressed as (more general formulas are (4.30), (4.31) in [15]) The source term S takes into account the effects of gravity, hydrostatic pressure gradient, and pressure-dependent friction. A positive S has a decelerating effect on the velocity (resistive net force), while a negative S has an accelerating effect (driving net force). The following boundary conditions are set for all t > 0. At the static/flowing interface Z = b, the velocity vanishes and the shear stress Σ XZ is continuous across it, so that where the values at Z = b are to be understood as the limit as Z → b with Z > b. At the free surface Z = h, the viscous stress vanishes, see (2.5), hence Since the material below the static/flowing interface is at rest, the boundary conditions (2.8) mean that the velocity is continuously differentiable across the interface in the viscous case ν > 0, whereas the velocity is only continuous across the interface in the inviscid case ν = 0. An additional "static equilibrium condition" completes the model. It states that the net forces at the interface Z = b (including friction) are resistive, so as to maintain the static part of the flow at rest. The condition (2.10) together with the interface conditions (2.8) indeed replace the momentum balance (2.3), (2.4) in the static part 0 < Z < b and across the interface (which contain the yield stress formulation in the static layer). A rigorous proof of equivalence between (2.8), (2.10) and the momentum balance in the static layer and across the interface is provided in the appendix of [15]. The quantities Φ(t, X), κ(t, X, Z) in that appendix correspond here to Φ = g(sin θ + ∂ X (h cos θ)), κ = µ s p. (2.11) The definition (2.7) of S corresponds thus to S = Φ − ∂ Z κ. From that appendix, note in particular the value of Σ XZ in the static layer Σ XZ (t, X, Z) = κ(t, X, b(t, X)) − (b(t, X) − Z)Φ(t, X), for 0 < Z < b(t, X), (2.12) that makes the net force in (2.3) identically vanish in the static layer. The condition (2.10) comes from (A12) in [15], and it expresses exactly that Σ XZ in (2.12) is below the yield stress κ in the static layer Z < b. The proof of the appendix of [15] is established under the assumption that the pressure p is convex with respect to Z, or equivalently that ∂ Z S ≤ 0 according to (2.7). If this condition is not satisfied, the rheological model (2.1) could be unstable and lead to shear band instabilities such as those that can be observed in simulations, see [40]. In the latter case we do not expect our thin-layer model to be valid. An implicit assumption made here (and in the appendix of [15]) is that the velocity is continuous across the interface. This is obviously necessary in the viscous case ν > 0. Although one could think that it is not strictly necessary in the inviscid case, we nevertheless only consider continuous velocities. Laboratory experiments or simulations of the viscoplastic model always show continuous velocity profiles. Indeed, velocities which are discontinuous across a line are not allowed with the rheology (2.1), it is proved in [34] that even without viscosity, there is a solution for which the velocity has square integrable first-order derivatives in space. Furthermore, it can be observed with simulations such as those done in [28] that, as ν → 0, the (continuous) solution to the present thin-layer model with ν > 0 tends to the continuous solution of the thin-layer model with ν = 0. Thus if discontinuous solutions exist, they are unphysical and unstable by perturbation. We shall see in Section 3 that continuous solutions exist for all the initial conditions we consider here.
An important observation in our above model is that the time evolution of the interface position b(t, X) is implicitly governed by the boundary conditions (2.8), (2.9), (2.10). It is thus coupled to the evolution of the velocity u. Note also that if S(t, X, Z) is negative for all Z, the condition (2.10) cannot be satisfied, meaning that there is no solution with static/flowing interface b; instead the whole granular mass should flow down immediately (this situation is not the object of our study). On the contrary if S(t, X, Z) is positive for all Z, then the net force is resistive, leading to progressive stopping of the flow.
Because of the Z-derivative in (2.7) and of variations at the scale of the thin layer (see [15] for the orders of magnitude), the pressure p has to be expressed at higher order with nonhydrostatic corrections. Neglecting the viscosity, [15,Theorem 4.1] gives The pressure is thus coupled to the unknown velocity u and depends nonlinearly on ∂ Z u.
Plugging this expression into the expression (2.7) of S, we obtain 14) The first term S 1 is independent of Z and takes into account the gravity and the hydrostatic part g cos θ(h − Z) of p in (2.13). It is the leading term. Considering only this term in S, the static equilibrium condition (2.10) is reduced to | tan θ| ≤ µ s . The second term S 2 is also independent of Z and is a correction to S 1 due to down-slope variations of h. Finally, the third term S 3 involves velocity derivatives and results from nonhydrostatic effects and down-slope space inhomogeneities. This term is responsible for the Z-dependence of S. The aim of this paper is to investigate the possible effects on the flow of the dependence of S on Z due to the last term S 3 in (2.14). The relevance of the nonhydrostatic correction in (2.13) has been confirmed in [40,VIII(B)].
If the material is flowing in a channel of width W , friction along the lateral walls of the channel can be taken into account by adding a Coulomb friction term that is proportional to the pressure divided by the channel width (see equations (19) and (17) of [40]). This should lead to an additional term in S with Z-dependency since the pressure depends on Z. Thus in this situation, the Z-dependency of S would result both from nonhydrostatic effects and from lateral friction. Indeed friction increases with pressure, thus with depth, leading possibly to a positive S near the bottom even though a negative S could occur near the surface. These wall effects are not included in the present study because the additional momentum friction term depends on the (multivalued) sign of the velocity and therefore requires a special analysis of the static equilibrium conditions in the solid layer.
Note that the asymptotics of [15] assumes that | tan θ| and µ s are close (the difference being at most of the order of the thickness of the layer). This is necessary in order for S (indeed S 1 ) to be of moderate amplitude and it is used to prove that the static/flowing interface continues to exist and to move gradually at positive times. Otherwise the whole material layer would flow down immediately (| tan θ| ≫ µ s ) or would stop very quickly (| tan θ| ≪ µ s ). In particular, if | tan θ| = µ s and ∂ X h = 0, the source term is reduced to the third term in (2.14). A rough idea of the values of this third term can be obtained by considering the velocity profile derived in [17], where u = u(t, X) is an average velocity. However, the expression obtained is rather complicated and, in any case, this profile u( η) is anyway not a solution to (2.6). The errors produced when modifying the velocity profile are discussed after equation (3.22). Since the shape in Z of the last term in (2.14) is unknown, although it should satisfy ∂ Z S ≤ 0 (see the comments after equation (2.12)), we shall consider generic sources S satisfying this condition.

Simplified model with variable source term
In this section we introduce our simplified model and prove a few of its qualitative properties.

Derivation of the simplified model
The key idea behind the derivation of the simplified model from the asymptotic model of the previous section is the decoupling of the coordinate X in the down-slope direction from the coordinate Z in the normal direction. Considering that X is fixed, the equation (2.6) with unknowns u and b and with boundary conditions (2.8), (2.9), (2.10) is a closed system in the variables (t, Z) set in the flowing layer b < Z < h, provided that h and S are known. Assuming that h does not depend on time, we can consider the system in the (t, Z) variables in the flowing layer, written as follows. The equation for the velocity u in the down-slope direction is where h > 0 is the fixed thickness of the domain and ν ≥ 0 is the kinematic viscosity. The viscosity ν can be time and space dependent, in order to include the case of µ(I) rheology, as described in [27]. The source term S(t, Z) is assumed to be given, defined for all t ∈ [0, T ] and all Z ∈ [0, h] and satisfying The values of S that are meaningful are only those for b(t) ≤ Z ≤ h, but since b(t) is unknown, we prefer to assume that S is given for all 0 ≤ Z ≤ h. In (3.1), the position b(t) ∈ (0, h) of the static/flowing interface is an unknown of the problem, along with the velocity u(t, Z) in the flowing domain b(t) < Z < h. We complete the equation (3.1) with the following boundary conditions for all t > 0, Note that in the present formulation of the model, u is only defined for Z ≥ b and therefore the derivatives of u at b are always understood as being taken from above, i.e. from the flowing side. In the inviscid case ν = 0, only (3.3a) remains since the other two conditions hold trivially. Typical expected velocity profiles are presented in Figure 2. The initial condition is given by where b 0 ∈ [0, h] is the given initial position of the static/flowing interface (i.e. b(0)). Indeed one would like that b(0 + ) = b 0 , but we shall see that this is not always possible, see in particular Subsection 4.5. Indeed in the formulation (3.1), (3.3), there is no obvious reason for b to be continuous with respect to time. The initial velocity profile u 0 is assumed to satisfy so that the initial velocity increases over (b 0 , h). Whenever useful we extend u 0 by zero in the interval [0, b 0 ]. We always assume that b(t) ∈ (0, h) for all t ≥ 0. The special case where b(t) reaches one of the boundaries 0 or h will be discussed further on whenever it becomes relevant. The simplified model is completed with the "static equilibrium condition" (2.10), that becomes here  .6) is obtained by considering a constant source S > 0 with ν = 0, a case that has been considered in [10].
Remark 3.2. In the (t, Z) formulation (3.1), it is not possible to replace the source term S by its velocity-dependent value (2.14). This is because X is fixed and we cannot express the X-derivatives of h and u. We thus have to consider an empirical source term S that we assume to be given. This means that, in our simplified model, all the X-dependences are somehow hidden in the shape of the source term S(t, Z). The main scope of this paper is to describe the effect of the shape of S in terms of Z on the dynamics of both the velocity u and the interface b.

Properties of the simplified model
We need to set some conditions so that the simplified model (3.1), (3.3), (3.4), (3.6) with unknowns u(t, Z) and b(t) defined for Z > b(t) will be consistent with the original asymptotic model of Section 2. For that model, we expect that for all t ∈ [0, T ], the velocity is increasing in the flowing part (b(t), h). This will be the case under the assumption that the source term is nonincreasing in space. As discussed after equation (2.12), this condition is necessary for the validity of our asymptotic model with a source term. More explicitly, we require that Notice that we employ hereafter the following mathematical terminology.
In the results below, we assume implicitly that u(t, Z) is sufficiently smooth in its definition domain b(t) ≤ Z ≤ h. We will build such smooth solutions in Section 4. . Moreover, assume that the viscosity ν is nonzero, that the initial condition satisfies (3.5a), and that the source term S satisfies (3.7). Then we have that Proof. By differentiating (3.1) with respect to Z, we obtain Together with (3.7), this leads to With the boundary conditions (3.3b), (3.3c) and (3.4), (3.5a), this leads to Knowing that ν > 0, this formulation is suitable for applying the strong maximum principle to the function ∂ Z u. In order to deal with a time-independent interval, we perform the following change of coordinates As ν > 0, we deduce that V > 0, i.e. (3.8) holds true.
Our second result shows that in the viscous case, the "static equilibrium condition" (3.6) is automatically satisfied. . Moreover, assume that the viscosity ν is nonzero and that (3.2), (3.5a) and (3.7) hold. Then (3.6) holds. (3.14) Using (3.3b) and the assumption ν > 0, we get ∂ t u(t, b(t)) = 0. Then evaluating (3.1) at Our third result is an explicit differential equation on the interface position b(t) that holds under regularity assumptions on the velocity and the source term.
Proof. We again differentiate (3.3a) with respect to time, leading to (3.14).
1) If ν > 0, as in the proof of Lemma 3.2, we get (3.15) and the first identity. By differentiating the boundary condition (3.3b) with respect to time, we similarly obtain Then we differentiate (3.1) with respect to Z and evaluate it at (t, b(t)). We obtain thus with (3.18), . Thus, by continuity of ∂ t u and S in space, we obtain ∂ t u(t, b(t)) = −S(t, b(t)). With (3.14) and the assumption ∂ Z u(t, b(t)) = 0 this gives the result (3.17).
When the denominator does not vanish in (3.16), (3.17), these differential equations provide valuable information on the nature of the interface dynamics: it is driven by the normal derivative of the velocity at the interface in the inviscid case and by the third normal derivative in the viscous case. However, the formulas can only be evaluated when the Zprofile of the velocity is known, which means we must also know the solution u of our model. The viscous formula (3.16) can also be written aṡ where S tot = S − ∂ Z (ν∂ Z u) is the total (opposite) net force, including the viscous force from (3.1). Equations (3.16), (3.17) also indicate that, as discussed in the introduction, a model that does not approximate the shape of the velocity profile well (i.e. not only the mean value but also the derivatives), will not be able to accurately predict the position of the interface. This is illustrated in [2] where the approximate value of the velocity is reproduced by the model throughout the flow thickness, whereas the shape of the velocity profile is very different from the one calculated with discrete element models. In that case, Figure 10 of [2] shows that the static/flowing positions computed with the two models differ. A similar conclusion is obtained in [1] when comparing the interface dynamics in an Herschel-Bulkley fluid and in an averaged model.
Finally, note that although Lemma 3.3 provides a differential equation on b(t), that could be used instead of the overdetermined boundary conditions, it is often preferable to use the latter (see [35] where it is shown that only the overdetermined boundary conditions lead to numerical stability). Indeed, the overdetermined boundary conditions do not require strong regularity assumptions and remain valid if the denominator in (3.16), (3.17) vanishes. In particular, the formula (3.17) for the inviscid case givesḃ(t) ≥ 0 because of (3.6). However, this is wrong in general: we show in Proposition 4.7 that, in the inviscid case, it is possible to have solutions withḃ(t) < 0 for all t, together with ∂ Z u(t, b(t)) = 0 and S(t, b(t)) = 0 (thus leading to an undetermined formula "0/0" in (3.17)).

Analytical study of the inviscid case
In this section we consider the simplified model (3.1), (3.3), (3.4), (3.6) in the inviscid case, i.e. when ν = 0, and we build an explicit solution (u, b) in several configurations. We assume that the source term S(t, Z) admits a unique (time-dependent) zero b ⋆ (t) ∈ [0, h], i.e. S(t, b ⋆ (t)) = 0. We will prove below that b(t) ≤ b ⋆ (t). Thus the value Z = b ⋆ (t) divides the flowing layer (b(t), h) into an upper sublayer (b ⋆ (t), h) where the net force is driving and a lower sublayer (b(t), b ⋆ (t)) where the net force is resistive. One important outcome of the results presented in this section is that the function b ⋆ (t) provides a monitoring of the position b(t) of the static/flowing interface: as time goes on, b(t) has the tendency to approach b ⋆ (t).
We investigate several types of behaviour for b ⋆ (t) and initial conditions, which in turn allow us to model different possible dynamics. In each case, we provide an explicit solution by first building the velocity u, then the position b of the static/flowing interface. Illustrations are given for a linear initial velocity profile and a source term that is linear in space. Moreover, as in Section 3 we assume that the source term S is continuous in time and space and is nonincreasing in space, i.e. (3.2), (3.7) hold. The initial condition (3.4) writes

Simplified model in the inviscid case
where u 0 satisfies (3.5) and b 0 ∈ [0, h] is the given initial position of the interface. We seek for a solution (u, b) that satisfies the monotonicity condition and the static equilibrium condition (3.6), i.e.
The source term S(t, Z) is assumed to have a unique (time-dependent) zero b ⋆ (t) ∈ [0, h], hence satisfying S(t, b ⋆ (t)) = 0. Together with the monotonicity of S in Z resulting from (3.7), this assumption is equivalent to 6a) S(t, Z) < 0 for all Z > b ⋆ (t). (4.6b) Moreover, since S is continuous with respect to its two arguments, b ⋆ (t) is continuous likewise. Furthermore, we get the following reformulation of the static equilibrium condition (4.5).
Lemma 4.1. Assume that the source term S is continuous with respect to its two arguments (3.2), satisfies the monotonicity condition (3.7) and admits a unique zero b ⋆ (t) in space. Then, the condition (4.5) is equivalent to Proof. This is immediate with (4.6).
Let us give some comments. the inequality (4.7) provides a barrier for the static/flowing interface that prevents the flow from reaching a full stop.
With the above assumptions, we are now going to provide an explicit solution (u, b) to (4.1), (4.2), (4.3), (4.4), (4.5), in several configurations. First we investigate the case of a nondecreasing zero b ⋆ , with the initial condition b 0 < b ⋆ (0) and we seek for a continuous and nondecreasing b(t). Next, we consider a decreasing zero b ⋆ , with b 0 = b ⋆ (0) and we seek for a continuous and decreasing b(t). Then, we handle a more general case than the previous one by considering a decreasing zero b ⋆ , but with b 0 < b ⋆ (0) and we seek for a continuous b(t). Finally, we consider a nondecreasing zero b ⋆ , with b 0 > b ⋆ (0). In this last case, the position b(t) of the static/flowing interface is discontinuous at the initial time, with b(0 + ) = b ⋆ (0).

Nondecreasing zero b ⋆ (t) and b 0 < b ⋆ (0)
We here investigate the case of a nondecreasing zero b ⋆ (t) of the source term S and we look for a position b(t) of the static/flowing interface that is also nondecreasing. This configuration allows us to model a progressive stopping of the flow.
We consider a source term S with a unique zero b ⋆ (t) (i.e. (4.6)) that satisfies We seek for a continuous function b(t), with b(0) = b 0 . We first build the velocity u solution to the evolution equation (4.1), under the assumption that the position b(t) of the static/flowing interface is nondecreasing. Then, we deduce explicitly the static/flowing interface position b(t) by analyzing the boundary conditions. The following lemma provides the explicit expression of the velocity. In order for (4.1) to be valid over the interval of integration (τ, Z), τ ∈ [0, t], we need that b(τ ) ≤ Z for all τ ∈ [0, t]. Since b is nondecreasing, this condition reduces to b(t) ≤ Z. The integration path is presented in Figure 3. Integrating (4.1) yields (4.10). Next, in order to determine b(t), we consider an extended velocity defined on [b 0 , h] and denoted byũ. The extended velocity is the right-hand side of (4.10), that is (4.11) Lemma 4.2 says that u is equal toũ for all Z ∈ [b(t), h]. Since b(t) is nondecreasing, we are looking for b(t) ∈ (b 0 , h). In order to get the boundary condition (4.2), b(t) must satisfỹ u(t, b(t)) = 0. We note that the assumptions (3.7) and (3.5a) imply that ∂ Zũ > 0 for all Z ∈ (b 0 , h), hence (4.4) is automatically satisfied. The equations (4.1) and (4.3) are also consequences of the formula (4.10). Therefore, it remains only to clarify if one can find a continuous nondecreasing function b(t) such thatũ(t, b(t)) = 0, which satisfies b(0) = b 0 and (4.5). This is what gives the following lemma. Proof. The uniqueness follows from the previously stated property ∂ Zũ > 0. Then, we note thatũ is continuous in both t and Z and that by (4.9) we have b ⋆ (t) > b 0 . Hereafter, we prove the existence of b(t) ∈ [b 0 , b ⋆ (t)) such thatũ(t, b(t)) = 0. Indeed, we show thatũ(t, b 0 ) < 0 for t > 0 andũ(t, b ⋆ (t)) > 0. At first, by (3.5b), we havẽ Using again b ⋆ (t) > b 0 and (3.5), we get u 0 (b ⋆ (t)) > 0. Moreover, by (4.9a), b ⋆ (τ ) ≤ b ⋆ (t) for all τ ∈ [0, t], thus according to (4.6) we have S(τ, b ⋆ (t)) ≤ 0. Thereforeũ(t, b ⋆ (t)) > 0 as claimed. We deduce the existence of b(t) and it only remains to prove that it is continuous and nondecreasing. Sinceũ(t, b(t)) = 0 withũ continuous, we infer the continuity of b. Then, differentiating with respect to t gives ∂ tũ (t, b) + ∂ Zũ (t, b)ḃ = 0. Since ∂ tũ = −S, we have ∂ tũ (t, b(t)) = −S(t, b(t)) ≤ 0 by Lemma 4.1 and since ∂ Zũ > 0, we getḃ ≥ 0.
The synthesis of Lemmas 4.2 and 4.3 is stated in the following main result of this subsection.  and satisfies b(t) < b ⋆ (t) for all t ∈ [0, T ].

Illustration
For the purpose of illustrating our results, we investigate the most simple case of linear data. They are written in dimensionless units, meaning that appropriate rescaled variables have to be used. On the one hand we consider a linear initial velocity where b ⋆ (t) is continuous and satisfies (4.9). According to Proposition 4.4, we have a unique solution to our problem. The velocity is written, for all t ∈ [0, T ], (4.17) The expression of b(t) is determined by (4.14), leading to In accordance with our result, the function defined by (4.18) is nondecreasing, as can be checked directly with the assumption (4.9). Let us consider a constant value b ⋆ (t) = b ⋆ > b 0 for all t ∈ R + . Then (4.17) and (4.18) give and The expression (4.20) satisfies b 0 < b(t) < b ⋆ for all t > 0 and b(t) → b ⋆ as t → ∞. We remark that evaluating (4.19) meaning that the velocity profiles intersect at Z = b ⋆ for all t ∈ R + . The solution (u, b) is plotted in Figure 4. Since b(t) < b ⋆ < h, the static/flowing interface never reaches the top of the domain, modelling a flow that does never reach a complete stop.

Decreasing zero
We assume here that the zero b ⋆ (t) of the source term S is decreasing, starting from the initial position b 0 of the static/flowing interface. This configuration allows us to model a progressive start of the flow with a thickening fluid part.
, and (4.4) holds. For Z = b(t) one has t ′ = t, thus the latter equality provides the additional boundary condition (4.28) (it is also possible to use (4.23)). We can remark additionally that for t > 0, the limits of ∂ Z u on both sides of b ⋆ (0) = b 0 differ, Therefore, as long as ∂ Z u 0 (b 0 ) > 0, the velocity u is not continuously differentiable at b 0 .

Illustration
In order to illustrate the configuration of Proposition 4.7, we consider linear initial data and source term in the form (4.15) and (4.16), for some b ⋆ satisfying (4.22). We take for where 1 < b 0 < h. Then b = b ⋆ (t) and the velocity is determined for all t ∈ [0, T ], by

Decreasing zero
We now handle a more general case than previously in Subsection 4.3. We consider still a decreasing zero b ⋆ (t) of the source term, but starting from an initial value b ⋆ (0) larger than the initial position b 0 of the static/flowing interface. This configuration allows us to model a change of variation in the static/flowing interface position, from an initial thickening to a progressive thinning of the static part of the flow.
We consider a source term S with a unique zero b ⋆ (t) (i.e. (4.6)) that satisfies We seek for a continuous function b(t) with b(0) = b 0 . Let us first state a general result which follows from Lemma 4.5. Proof. According to Lemma 4.1, one has b(t) ≤ b ⋆ (t) for all t > 0. Thus the assumption b(t) < b ⋆ (t) for all t ∈ I means that b(t) = b ⋆ (t) for t ∈ I. Arguing by contradiction, the result is immediate with Lemma 4.5.
In the context of (4.35) and for a continuous b with b(0) = b 0 < b ⋆ (0), one has necessarily b(t) < b ⋆ (t) for t small enough. Then there are two possible situations: Case (i) means that there is no intersection between the graphs of b and b ⋆ . In case (ii), there is an intersection at a time t ⋆ . Then, from t ≥ t ⋆ we encounter the situation of Subsection 4.3. Herein, we focus on case (ii), that can be formalized as (4.36) Therefore, the position b of the static/flowing interface is increasing over [0, t ⋆ ] and is de- In what follows, we first build the velocity u under the assumption (4.36), then we determine the position b(t) of the static/flowing interface. We suppose that there exists t ⋆ ∈ (0, T ) such that (4.36) holds. Then, for all t ≥ t ⋆ , the velocity is given by For all t ∈ [0, t ⋆ ], the velocity is given by Proof. The geometric setting is presented in Figure 7. The domain Z ≥ b(t) where the velocity u is defined is filled with oblique lines. It is split in three sets A, B, C corresponding to (4.37a), (4.37b), (4.38) respectively. If t ≤ t ⋆ then b(t) is increasing, therefore the velocity is given by Lemma 4.2, which gives (4.38). If t ≥ t ⋆ , we are in the configuration of Subsection 4.3. Applying Lemma 4.6 we deduce the formula (4.37a) and But according to (4.38) applied for t = t ⋆ , we have Replacing the value of u(t ⋆ , Z) given by (4.40) in (4.39) gives (4.37b).
Proof. We denote by P rop(t) the following property: there exists a unique zero b(t) of Z → u(t, Z), such that b 0 < b(t) < b ⋆ (t). We consider the set I = t + > 0 such that ∀t ∈ (0, t + ), P rop(t) holds . (4.41) The set I is convex since for any t + 1 ∈ I, t + 2 ∈ (0, t + 1 ) ⇒ t + 2 ∈ I. Thus I is an interval, and one of the three situations must occur: In order to prove the first statement of the lemma, we have to establish that (a) cannot happen, i.e. that I = ∅. In other words we need to prove that for t > 0 small enough, u(t, .) admits a unique zero b(t) in (b 0 , b ⋆ (t)). The uniqueness follows from the property ∂ Zũ > 0, which is obtained with (3.5a) and (3.7). In order to get existence, we first notice that according to (4.35), for t small enough one has b ⋆ (t) > b 0 . Then, let us prove that for t small enoughũ(t, b 0 ) < 0 andũ(t, b ⋆ (t)) > 0. We have by the definition (4.11) and by (3.5b) thatũ(t, b 0 ) = − t 0 S(τ, b 0 )dτ . Thus, by (4.6) we deduce that S(τ, b 0 ) > 0 and that u(t, b 0 ) < 0. Next, since S is continuous by (3.2), there exists a constant C 1 ≥ 0 such that S(τ, Z) ≤ C 1 for all τ ∈ [0, T ] and all Z ∈ [0, h]. This leads toũ(t, b ⋆ (t)) ≥ u 0 (b ⋆ (t)) − C 1 t. This expression tends to u 0 (b ⋆ (0)) as t → 0. Since by (4.35b) and (3.5) this limit is positive, we deduce that for t small enough we haveũ(t, b ⋆ (t)) > 0, concluding the existence of a zero b(t) in the interval (b 0 , b ⋆ (t)) and the proof that I = ∅.
It remains to prove the last statement concerning the case t ⋆ < T . We are going to prove that in this case we haveũ(t ⋆ , b ⋆ (t ⋆ )) = 0, giving the claim. Indeed, the assumption t ⋆ < T yields by definition of t ⋆ that there exists a sequence t n ∈ (t ⋆ , T ) such that t n → t ⋆ and not satisfying P rop(t n ), which means thatũ(t n , .) has no zero in (b 0 , b ⋆ (t n )) (recall that uniqueness always holds because ∂ Zũ > 0).

Illustration
In order to illustrate the proposition, we consider a linear initial velocity in the form (4.15) and a linear source term in the form (4.16). We set  (4.45) and for all t ∈ [t ⋆ , T ], b(t) = b ⋆ (t) and

Increasing zero
We present here a case where the static/flowing interface position b(t) is initially discontinuous with respect to time. We assume that the zero b ⋆ (t) of the source term is increasing, starting from an initial value lower than the initial position of the static/flowing interface. This configuration allows us to model a sudden starting of a part of the initially static material, as well as a progressive stopping of the flow.
We consider a source term S with a unique zero b ⋆ (t) (i.e. (4.6)) that satisfies b ⋆ (t) is increasing, We recall that according to Lemma 4.1, we are looking for an interface position b that Together with (4.47b) and since b ⋆ is continuous, this condition implies that b(t) cannot tend to b 0 as t → 0 + , so that we have b(0 We are going to take the minimal possible jump. Thus we consider an interface position b that satisfies For consistency with this configuration, the initial condition (4.3) has to be modified to Proof. It is identical as that of Lemma 4.2. The domain of integration is presented in Figure  9.
* * * * We skip the proof since it is almost identical to that of Lemma 4.3. Using Lemmas 4.12 and 4.13 we conclude the following proposition. and satisfies b(t) < b ⋆ (t) for all t ∈ (0, T ].

Illustration
We consider linear data of the form (4.15) and (4.16), with b ⋆ continuous satisfying (4.47).
The initial velocity is extended over [0, h] by (4.51). We take with 0 < b ⋆ (0) < b 0 . Then according to Proposition 4.14, there is a unique solution (u, b) to our problem. We first determine the velocity for all t ∈ [0, T ] such that b(t) ≥ b 0 . By (4.52), in this situation we have (4.55) By writing the zero of this expression, we get This formula is valid if it satisfies b(t) ≥ b 0 , which can be rewritten .
Then, for all t ∈ (0, T ] such that (4.57) is not satisfied (which corresponds to small times), we have to look for b(t) such that b ⋆ (0) < b(t) < b 0 . By Lemma 4.12, the velocity is then given by We deduce the value of b(t) by getting the zero of (4.58b), (4.59) The velocity profiles are linear for all t ∈ [0, T ] such that b(t) ≥ b 0 . They are piecewise linear for all t such that b ⋆ (0) < b(t) < b 0 . The velocity and interface position are presented in Figure 10. The position b(t) of the static/flowing interface jumps from b 0 to b ⋆ (0) instantaneously. We also notice that b(t) increases to b ⋆ (t). This illustrates the monitoring effect of b ⋆ on b.

Overview and outlook
Starting from the non-averaged thin-layer asymptotic model of [15] describing flows of granular materials with static/flowing transition and viscoplastic rheology with Drucker-Prager yield stress, we have derived a simplified one-dimensional time-dependent model keeping the normal coordinate to the topography. The simplified model intrinsically describes the time evolution of the velocity profile with respect to the normal coordinate. It involves a source term that represents the opposite of the net force in the flowing layer, including gravity, pressure gradient, and internal friction. Given the source term, the simplified model gives at the same time the evolution of the velocity profile in the flowing layer and that of the position of the static/flowing interface. We have proved several properties of consistency of the simplified model. We have also performed an analytical study of the inviscid case for a source term variable in time and space. In this inviscid context, we have derived explicit solutions in various situations, showing the ability of the simplified model to represent different types of interface dynamics. In particular, we have shown that the shape of the source term, especially its zero b ⋆ (t) in space and its time dependence, is a key factor in determining the evolution of the interface. This zero divides the material layer into an upper sublayer (b ⋆ (t), h) where the net force is driving and a lower sublayer (0, b ⋆ (t)) where the net force is resistive. We have exhibited different dynamics of the flow, from starting to arrest, including progressive starting, progressive stopping and a sudden start. All of these behaviours represent relevant situations for real flows, although further studies are needed to compare the results with real data.
Our study indicates well-posedness of the reduced model described in Section 3, even if the initial Drucker-Prager model is ill-posed. One could hope then that when modifying the initial model to a well-posed one, the dynamics of the interface could be close to the one that we describe here.
The various behaviours can be summed up as follows. First, we assumed that the zero b ⋆ (t) of the source was nondecreasing with respect to time, and additionally that the initial data satisfies b 0 < b ⋆ (0). In this configuration, the position b(t) of the static/flowing interface increases until reaching the zero b ⋆ (t). Physically, this means that the static part of the flow thickens until reaching b ⋆ (t). If b ⋆ (t) remains smaller than the height h of the material layer, this prevents the flow from reaching a complete stop. On the contrary, if b ⋆ (t) reaches h, then the flow progressively reaches a complete stop. Such a situation could reflect the behaviour of flows over thick erodible beds confined in a channel, which are created by a constant flow rate upslope [47,30,20]. Indeed, in that case, steady uniform flows develop with a given flowing thickness above a static layer and a given free surface slope that depends on the imposed flow rate. Due to wall effects, the source term S should depend on Z and should have a zero that may change with time due to changes in the free surface gradient. The static/flowing interface position in the resulting steady uniform flow should stabilize at the position of the zero of S after a transition phase that possibly resembles the evolution described above.
In a second configuration, we assumed that b ⋆ (t) was decreasing, and additionally that the initial data satisfies b 0 = b ⋆ (0). In this case, the position of the static/flowing interface b(t) follows b ⋆ (t), i.e. b(t) = b ⋆ (t) for all times, and the velocity profiles have both a nonlinear and a linear part. As b ⋆ (t) decreases, b(t) also decreases. Therefore, the fluid part of the flow thickens. If the zero b ⋆ (t) of the source term tends to 0, then the position b(t) of the static/flowing interface tends to 0 as well, meaning that the whole material layer flows. Therefore, if b ⋆ (0) = h, this models a flow that starts progressively before fully flowing.
Next, we considered a decreasing zero b ⋆ (t) as previously, but with initial data satisfying b 0 < b ⋆ (0). In this configuration, the position b(t) of the static/flowing interface starts to increase (with b(t) < b ⋆ (t)) until reaching the zero b ⋆ (t) at some time t ⋆ . Then we are in the situation of the previous case and, for larger times, b(t) remains equal to b ⋆ (t). This models a progressive thickening of the static part, followed by a progressing thinning. Therefore, in this configuration, we obtain a change in the sign of the variation of the position of the static/flowing interface.
Finally, we assumed that the zero b ⋆ (t) of the source term was increasing, with the initial data satisfying b 0 > b ⋆ (0). In this case, b(t) jumps instantaneously to the value b ⋆ (0), i.e. b(0 + ) = b ⋆ (0), and then b(t) increases with b(t) < b ⋆ (t). This special case involves a static/flowing interface position which is initially discontinuous with respect to time. It models a sudden starting of a part of the granular mass. This behaviour may reflect what can happen when a granular avalanche in a channel flows over an initially static layer, as in the granular collapse experiments of [39,23]. In that case, due to wall effects, S possibly changes sign along the Z-axis and the zero of S could possibly change in time due to changes in the free surface gradient. We indeed observe a relatively rapid starting of part of the initially static grains of the erodible bed that may be due, under certain conditions, to the scenario described above. The starting of the granular mass initially at rest is however more progressive in the experiments, possibly because of viscous effects (e.g. see Figure 16 of [35]).
All of the above configurations lead to the following two important conclusions. First, the zero b ⋆ (t) of the source plays the role of a barrier, since as stated in Lemma 4.1, we must have b(t) ≤ b ⋆ (t). Thus it can prevent the flow from reaching a complete stop. Second, the interface position b(t) always "tries to follow" the zero b ⋆ (t). If b ⋆ (t) increases, then b(t) also increases, while remaining lower than b ⋆ (t), and b(t) can reach b ⋆ (t) in large time. If b ⋆ (t) decreases, then at some time t ⋆ , we have b(t ⋆ ) = b ⋆ (t ⋆ ) and then b(t) = b ⋆ (t) for all times t ≥ t ⋆ .
These dynamical properties of the static/flowing interface should be interpreted more generally in the context of the coupling of the source term with the velocity dependence in the down-slope coordinate, as derived in the full thin-layer asymptotic model of [15]. This coupling with (2.14) needs to be done, at least numerically. Introducing wall effects in the model would make it possible to gain deeper insight into the static/flowing interface dynamics for flows confined in a channel, as is the case in most laboratory experiments. In particular, this could help us better understand the transition to steady uniform flows on heaps or erosion processes of granular flows over initially static beds. This will be the object of forthcoming studies. Note that [35] presents a numerical study of the simplified model in the viscous case with a constant source term and a comparison with laboratory experiments.