# Thick surface flows of granular materials: The effect of the velocity profile on the avalanche amplitude

## I General principles

### A Onset of avalanches

It is a dailylife experience that the top surface of a mass of granular matter need not be horizontal unlike that of a stagnant liquid. However, there exists an upper limit to the slope of the top surface, and the angle between this maximum slope and the horizontal is known, for non-cohesive material, as the Coulomb critical angle $`\theta _{max}`$. Above this angle, the material becomes unstable, and an avalanche at the surface might occur. The Coulomb angle is related to the friction properties through $`\mathrm{tan}\theta _{max}=\mu _i`$ where $`\mu _i`$ is an internal friction coefficient \[\].

As of today, the physical picture associated with the onset of the avalanche is still obscure. One could imagine a local scenario in which the dislodgement of some unstable grains leads by amplification to a global avalanche (see for instance \[\]). Alternately, one can think of a delocalized mechanism \[\], in which a thin slice of material is destabilized and starts to slide as a whole. In the present paper, we will focus on the latter point of view.

It has been recently suggested \[\] that the thickness of the initial gliding layer should be of the order of $`\xi `$, the mesh size of the contact force network \[\]. For simple grain shapes (spheroidal), one expects $`\xi `$ 5–10 grain diameter $`d`$. The angle at which the avalanche process actually starts is of the order of $`\theta _{max}+\xi /L`$, where $`L`$ is the size of the free surface. At the moment of onset, our picture is that this initial layer starts to slip, and is rapidly fluidized by the collisions with the underlying heap, therefore generating a layer of rolling grains on the whole surface.

Now that we have proposed a description of the initial situation, we may turn to the model scheme accounting for the further evolution of the avalanche.

### B Saturation effects for thick avalanches

Some years ago, Bouchaud, Cates, Ravi Prakash and Edwards introduced a model to describe surface flows of granular materials \[\]. The model assumes a rather sharp distinction between immobile particles and rolling particles and, accordingly, introduces the following two important physical quantities (see Fig. 1): the local height of immobile particles $`h(x,t)`$ (where $`x`$ denotes the horizontal coordinate \[\] and $`t`$ the time), and the local amount of rolling particles $`R(x,t)`$.

The time evolution of $`h(x,t)`$ is written in the form

$$\frac{h}{t}=\gamma R(\theta _n\theta )$$
(1)   
where $`\theta \mathrm{tan}(\theta )=h/x`$ is the local slope, $`\gamma `$ a characteristic frequency and $`\theta _n`$ the neutral angle of grains at which erosion of the immobile grains balances accretion of the rolling grains. For the rolling particles, Bouchaud and co-workers wrote a convection-diffusion equation \[\] that was later simplified by de Gennes as \[\]

$$\frac{R}{t}=v\frac{R}{x}\frac{h}{t}$$
(2)   
where $`v`$ is the downhill typical velocity of the flow, and is assumed to be constant.

According to the Bouchaud-Cates-Ravi Prakash-Edwards (BCRE) model, $`h/t`$ is linear in $`R`$ \[see Eq. (1)\]. This is natural at small $`R`$, when all the rolling grains interact with the immobile particles. But as explained in Refs. \[\], this cannot hold when $`R`$ becomes larger than a given saturation length $`\xi ^{}`$, since the grains in the upper part of the rolling phase are no longer in contact with the immobile grains. The length $`\xi ^{}`$ is expected to be of the order of a few grain diameters $`d`$ \[\]. This led Boutreux, Raphaël, and de Gennes to propose \[\] a modified saturated version of the BCRE Eq. (1), valid for thick surface flows and of the form

$$\frac{h}{t}=v_{uh}(\theta _n\theta )(R\xi ^{})$$
(3)   
where $`v_{uh}`$ is defined by $`v_{uh}\gamma \xi ^{}`$. The constant $`v_{uh}`$ has the dimensions of a velocity.

The description of thick avalanches modelized by Eq. (3) was discussed in Ref. \[\]. However, one might encounter situations where the local amount $`R`$ of rolling particles is rather large except in some regions of space where it takes values smaller than $`\xi ^{}`$. For such cases, various ‘generalized’ forms of the BCRE equations valid both in the large and small $`R`$ limit, and able to handle intermediate values have been proposed \[\]. As we will be concerned only with thick flows, we will henceforth use the saturated form (3).

### C Velocity profiles in thick flows

We now consider the hypothesis made in Eq. (2) that the downhill typical convection velocity of the rolling grains $`v`$ is constant. As a matter of fact, $`v`$ might vary for two reasons.

First, $`v`$ depends on the local slope $`h/x`$ of the static bed, reflecting that the mean convection velocity should increase as the sandpile is further tilted. However, in the situation we are going to consider, the slope should never depart from $`\theta _n`$ by more than a few degrees, so that the variations of $`v`$ originating in this may reasonably be taken to be negligible.

Second, $`v`$ might as well depend on the local amount of rolling particles $`R`$. This dependence is quite natural, since as soon as the thickness of the flow exceeds a few grain diameters, one would expect a velocity gradient perpendicular to the flow to establish. Such a possibility was already considered by Bouchaud et al. \[\], but, to our best knowledge, not further studied. We think that taking this velocity gradient into account does lead to an improvement of the model description of avalanches. In the forthcoming sections we will analyse the physical consequences of this modification.

If analyticity is assumed, we can expand $`v(R)`$ in powers of $`R`$, and considering only the first two contributions to be significant, we write:

$$v(R)=v_0+\mathrm{\Gamma }R.$$
(4)   
with $`\mathrm{\Gamma }`$ a constant, homogeneous to a shear rate, and $`v_0`$ a constant velocity.

When $`R`$ becomes small, Eq. (4) tells us that $`v(R)`$ becomes constant \[$`v(R)v_0`$\]. Physically this velocity should correspond to the typical convection velocity of a single grain on a bed of immobile grains. For simple grain shapes (spheroidal) and average levels of inelastic collisions, one expects this velocity $`v_0`$ to scale as $`(gd)^{1/2}`$ (where $`g`$ is the gravity) \[\]. Similarly, the shear rate $`\mathrm{\Gamma }`$ is expected to scale as $`(gd)^{1/2}/d(g/d)^{1/2}`$ \[\]. We can therefore rewrite Eq. (4):

$$v(R)=\mathrm{\Gamma }(R+d).$$
(5)   
We note that $`v_0`$ becomes negligible compared to $`\mathrm{\Gamma }R`$ as soon as $`R`$ exceeds a few grain diameters.

In our approach, the typical velocity $`v(R)`$ depends linearly on the local rolling height $`R`$ \[Eq. (5)\]. Such a form is in part motivated by the recent work of Douady et al. \[\] (see also Section IV C). It is also supported by the experimental results of Rajchenbach et al. who carried measurements in a rotating drum \[\]. These authors have found linear velocity profiles in the surface flow, with a shear rate $`\mathrm{\Gamma }`$ independent of the thickness of the flow. However, in other experiments of chute flows carried on rough inclined planes, Anzaza et al. \[\] and Pouliquen \[\] observe that the mean velocity (averaged on cross-sections) scales as a power-law of the thickness with an exponent about 3/2. In the following we will mainly focus on the linear form (4), since it allows us to give explicit analytical solutions, and shall discuss the changes that are to be brought in the case of a power-law velocity in Section IV A.

In the next section, we will derive the governing equations from the saturated BCRE equations and the above considerations on the velocity profile inside the flow.

### D Governing equations

We may define a reduced profile $`\stackrel{~}{h}`$, deduced from $`h`$ by substracting the ‘neutral’ profile $`\theta _nx`$:

$$\stackrel{~}{h}(x,t)h(x,t)\theta _nx.$$
(6)   
Using Eqs (2)(3)(5) and (6), we easily obtain the following system:

$`{\displaystyle \frac{\stackrel{~}{h}}{t}}`$ $`=`$ $`v_{uh}{\displaystyle \frac{\stackrel{~}{h}}{x}}`$ (7)   
$`{\displaystyle \frac{R}{t}}`$ $`=`$ $`\mathrm{\Gamma }(R+d){\displaystyle \frac{R}{x}}+v_{uh}{\displaystyle \frac{\stackrel{~}{h}}{x}}.`$ (8)   
In our approach, Eqs. (7) and (8) are the governing equations for surface avalanches displaying linear velocity profiles.

An important point is that we must have $`R>0`$ for Eqs. (7) and (8) to be valid. If we reach $`R=0`$ in a certain spatial domain, then Eq. (7) must be replaced in that domain by $`\stackrel{~}{h}/t=0`$.

## II Application to the simple case of an open system

### A Physical situation

We will now solve Eqs. (7) and (8) in the following simple situation: we consider a cell, of dimension $`L`$, partially filled with monodisperse grains of diameter $`d`$, as shown on Fig. 2. The heap has an initial uniform slope $`\theta _{max}`$, the Coulomb angle of the material. The origin of the $`x`$ axis is taken at the bottom of the cell, and the orientation of the axis is such that the slope of the heap is positive.

We now consider that an avalanche has just started in the cell (see Section I A), so that we have at time $`t=0`$ a layer of rolling grains in the whole cell, of thickness $`\xi `$ greater than the saturation length $`\xi ^{}`$. We may thus use the saturated equations (7) and (8) from the beginning of the avalanche.

As the rolling population will rapidly grow and become independent of the initial thickness $`\xi `$ (for $`\xi `$ small), we can as well consider the initial condition on $`R`$ to be:

$$R(x,t=0)=0.$$
(9)   
We also know the initial value of $`\stackrel{~}{h}`$:

$$\stackrel{~}{h}(x,t=0)=(\theta _{max}\theta _n)x\eta x,$$
(10)   
where $`\eta `$ is defined as the (positive) difference between the Coulomb angle and the neutral angle.

We have additional conditions in our system, due to the boundaries. At the top of the cell, there is no feeding in rolling species, so that we impose:

$$R(x=L,t)=0\text{at any time }t0.$$
(11)   
Another condition arises from the fact that grains fall off the cell at the bottom and cannot accumulate there:

$$\stackrel{~}{h}(x=0,t)=0\text{at any time }t0.$$
(12)   
### B Uphill wave in the static phase

Equation (7) can be readily solved along with conditions (10) and (12) to give:

$$\stackrel{~}{h}(x,t)=\eta v_{uh}H(xv_{uh}t)x\text{ for }0xL,$$
(13)   
where $`H`$ denotes the Heavyside unit step function \[$`H(u)=1`$ if $`u>0`$, $`H(u)=0`$ otherwise\]. This result corresponds to the uphill propagation (at constant speed $`v_{uh}`$) of a surface wave on the static phase. Let us call $`x_{uh}(t)`$ the time-dependent position of the wavefront, given by:

$$x_{uh}(t)=v_{uh}t.$$
(14)   
(where the subscript $`uh`$ stands for ‘uphill’).

The wave starts from the bottom of the cell at time $`t=0`$ and reaches the upper end at time $`t_2`$ defined by:

$$t_2L/v_{uh}.$$
(15)   
At a given time $`t`$ (smaller than $`t_2`$), the profile of the static phase can be described as follows: ahead of the wavefront \[$`x_{uh}(t)xL`$\], the profile is linear and the slope is the initial angle $`\theta _{max}`$ (since $`\stackrel{~}{h}=\eta v_{uh}`$). Behind the wavefront \[$`0xx_{uh}(t)`$\], the slope has decreased and reached the neutral angle $`\theta _n`$ ($`\stackrel{~}{h}=0`$) (see Fig. 3). For times $`tt_2`$ the slope of the static phase inside the cell is uniformly equal to the final value $`\theta _n`$, which is thus the angle of repose of our specific open cell system \[\].

### C Downhill convection of rolling grains

Substituting Eq. (13) into the evolution equation (8) for $`R`$ gives:

$$\frac{R}{t}\mathrm{\Gamma }(R+d)\frac{R}{x}=\eta v_{uh}H(xv_{uh}t)$$
(16)   
Eq. (16) is a non-linear convection equation. The rolling species are thus convected downhill, with a convection velocity dependent on the local rolling thickness $`R`$. In the spatial region $`x>v_{uh}t`$, the right-hand side (which couples the evolution of $`R`$ to that of $`\stackrel{~}{h}`$) plays the role of a source term, leading to an amplification of the avalanche. On the contrary, for $`xv_{uh}t`$, the right-hand side goes to zero, so that the material flowing through the surface $`x=v_{uh}t`$ from uphill is simply convected, without amplification nor damping.

Equation (16) can be solved analytically by using the method of characteristics \[\], which utilizes the property that certain types of partial differential equations reduce to a set of ordinary differential equations along particular lines, known as the characteristic curves. For more details on this method and on its application in the case of Eq. (16), see Appendix.

### D Propagation of boundary effects in the cell

Before we go to the precise solutions, we can try to get some physical insight of the way the avalanche is going to develop. The global shape of the rolling phase at different moments during the avalanche is of course very dependent on the boundary condition (11) for $`R`$ in the cell, but also on the condition (12) for $`\stackrel{~}{h}`$, since the evolution of $`\stackrel{~}{h}`$ and $`R`$ are coupled.

However, the effects of these boundary conditions cannot spread over the entire cell instantly after the beginning of the avalanche, and shall propagate with finite velocities. We then expect the progression of these boundary effects (one could say, the propagation of the ‘information’ on the boundaries) to control the evolution of both the rolling and the static phase. For instance, in the case of the static profile $`\stackrel{~}{h}`$, Eq. (13) tells us that the bottom boundary condition (12) \[$`\stackrel{~}{h}(x=0,t)=0`$ at any time $`t`$\] brings progressively $`\stackrel{~}{h}`$ to zero everywhere in the cell, and also, that the propagation proceeds with a velocity $`v_{uh}`$.

Hence, we expect that the description of the avalanche should naturally split up in different temporal ‘stages’, according to the degree of extension of the different boundary effects, and that the cell should divide in several ‘regions’, according to whether it is under the influence of the top boundary condition or the bottom, or both, etc. This shall become clear as we will now go into the precise description of the avalanche.

## III Unfolding of the avalanche

### A Stage I: The avalanche grows to maturity

This stage starts at $`t=0`$ with the beginning of the avalanche. From the above considerations, we know that the boundary effects start to propagate with finite velocities from both ends of the cell. We can therefore define ‘propagation fronts’ for these effects: we call $`x_{dh}(t)`$ the position of the front originating in the boundary condition at the top of the cell \[the subscript $`dh`$ means that the motion of this front is downhill\], and $`x_{uh}(t)`$ the corresponding ‘uphill’ front, originating in the bottom boundary condition, and that we already defined earlier as $`x_{uh}(t)=v_{uh}t`$ \[Eq. (14)\]. Figure 4-a presents a typical picture of the situation during Stage I, where the fronts, after leaving their respective cell ends, move in opposite directions and one toward the other. As a consequence, they shall finally meet at a certain time, that we hereafter denote $`t_1`$. This time $`t_1`$ defines the end of what we call ‘Stage I’ (which is thus characterized as the time interval $`0tt_1`$), and the beginning of ‘Stage II’ (described in next section).

The relative positions of the fronts naturally define three spatial regions in the cell (Fig. 4-a). To the left of $`x_{uh}`$, the effects of the bottom boundary condition \[Eq. (12)\] are predominant. We call this region the bottom region. We remark that it is constantly extending uphill during Stage I \[following the motion of $`x_{uh}(t)`$\]. To the right of $`x_{dh}(t)`$, we define the top region, which extends downhill, and where the evolution of the avalanche is controlled by the upper end condition \[Eq. (11)\]. Finally, between those two regions remains a central region, where none of the boundary effects can yet be felt. This last region shrinks during Stage I, and ultimately disappears at time $`t=t_1`$ when the bottom and top region connect. We now describe the precise evolution of the avalanche, region by region \[we will only give the form of the rolling amount $`R(x,t)`$, since $`h(x,t)`$ is already known from Eq. (13)\].

#### 1 Top region

From the above definition, the top region corresponds to the spatial domain $`x_{dh}(t)xL`$. Within this domain, Eq. (16) reads:

$$\frac{R}{t}\mathrm{\Gamma }(R+d)\frac{R}{x}=\eta v_{uh}$$
(17)   
\[since $`x>x_{dh}(t)>x_{uh}=v_{uh}t`$\].

Solving this equation with the boundary condition $`R(x=L,t)=0`$ (see Appendix for details) gives the expression of $`R`$ valid in this region:

$$R(x,t)=d+\sqrt{d^2+2(Lx)\frac{\eta v_{uh}}{\mathrm{\Gamma }}}.$$
(18)   
We also obtain the precise position of the ‘downhill front’ $`x_{dh}(t)`$:

$$x_{dh}(t)=L+\frac{\mathrm{\Gamma }}{2\eta v_{uh}}d^2\frac{1}{2}\mathrm{\Gamma }\eta v_{uh}\left(t+\frac{d}{\eta v_{uh}}\right)^2.$$
(19)   
Thus, according to Eq. (18), $`R`$ has a stationary shape (independent of time), but on a domain that extends downhill with time. Interestingly, we note that the motion of $`x_{dh}(t)`$ is uniformly accelerated throughout the stage. This is a direct consequence of the non-linearity in Eq. (16).

#### 2 Central region

In the central region, the boundary conditions have no influence on the evolution of $`R`$ and $`\stackrel{~}{h}`$. The central region is hence spatially defined by $`x_{uh}(t)xx_{dh}(t)`$, and shrinks at both ends to disappear at the end of the stage.

The evolution equation for $`R`$ is the same as in the top region \[Eq. (17)\], but now we must impose the initial condition (9) (no boundary condition). The solution (see Appendix) writes:

$$R(x,t)=\eta v_{uh}t.$$
(20)   
In the central region the rolling phase grows linearly with time and uniformly in space, thus forming a plateau (see Fig. 4-b). This constant growth rate is a consequence of the saturated form of the BCRE equations we have used \[Eq. (3)\]. The uniformity of the solution, on the other hand, stems from the fact that since none of the boundaries is at work, and since the initial static profile was uniform, the central region behaves like an infinite medium for which translational invariance is to be obeyed.

We finally remark that solutions (18) and (20) connect continuously at $`x=x_{dh}(t)`$.

#### 3 Bottom region

The bottom region is controlled by the bottom boundary condition, and spreads over the spatial interval $`0xx_{uh}(t)`$. In this region the evolution equation for $`R`$ displays no more amplification (because $`x<x_{uh}(t)=v_{uh}t`$): $`R/t\mathrm{\Gamma }(R+d)R/x=0`$.

Since there is no constraint on $`R`$ at the bottom of the cell, the condition on $`R`$ is given by the physical assumption that it should be continuous across the border of the central and bottom regions, i.e. $`R(x=x_{uh}(t),t)=\eta v_{uh}t`$ for $`t0`$.

This leads to the following expression for $`R`$:

$`R(x,t)={\displaystyle \frac{1}{2}}\left(d+{\displaystyle \frac{v_{uh}}{\mathrm{\Gamma }}}\eta v_{uh}t\right)`$ (22)   
$`+{\displaystyle \frac{v_{uh}}{2\mathrm{\Gamma }}}\sqrt{\left({\displaystyle \frac{\mathrm{\Gamma }d}{v_{uh}}}+1\mathrm{\Gamma }\eta t\right)^2+4{\displaystyle \frac{\mathrm{\Gamma }\eta }{v_{uh}}}(x+\mathrm{\Gamma }dt)}.`$   
In this region also the height of rolling grains increases with time, due to an increasing input of material at the frontier with the central region.

#### 4 Derivation of time $`t_1`$

Stage I ends when the top and bottom regions meet, at $`t=t_1`$ defined by $`x_{uh}(t=t_1)=x_{dh}(t=t_1)`$. Using Eqs. (14) and (19), we easily obtain

$$t_1=\left(\frac{d}{\eta v_{uh}}+\frac{1}{\mathrm{\Gamma }\eta }\right)+\sqrt{\left(\frac{d}{\eta v_{uh}}+\frac{1}{\mathrm{\Gamma }\eta }\right)^2+\frac{2L}{\mathrm{\Gamma }\eta v_{uh}}}.$$
(23)   
As will be shown later, the maximum thickness of the avalanche, $`R_{max}`$, is actually reached for $`t=t_1`$ and $`x=x_{uh}(t=t_1)=x_{dh}(t=t_1)`$. We can clearly see on Fig. 4 that the $`R`$-profile at time $`t=t_1`$ displays a cusp. As the prediction of the maximum amplitude $`R_{max}`$ is an important result of our analysis, we shall devote Section IV A to it, and defer the analytical derivation of $`R_{max}`$ and its application to physical examples until there.

Fig. 4-b presents successive ‘snapshots’ of the rolling phase profile during Stage I.

### B Stage II: The static profile reaches its final state

Stage II starts at $`t=t_1`$. At time $`t_1`$, the two ‘propagation fronts’ of the boundary effects $`x_{uh}(t)`$ and $`x_{dh}(t)`$ pass each other, and then pursue their respective motions towards the opposite cell edge. Figure 5-a illustrates this situation. As in Stage I, it appears that the cell is naturally divided in three spatial regions: a top region \[defined spatially as $`x_{uh}(t)<xL`$\], under the sole influence of the upper edge of the cell; a bottom region \[$`0x<x_{dh}(t)`$\], under the influence of the bottom edge; and finally, a central region \[$`x_{dh}(t)xx_{uh}(t)`$\], where in contrast with Stage I, the effects of both boundaries now combine. As another difference with the situation described in Stage I, the top and bottom regions progressively shrink, whereas the central one grows in extension (Fig. 5-a).

Due to their motion, the fronts $`x_{dh}`$ and $`x_{uh}`$ are bound to reach, sooner or later, the bottom and top end of the cell (respectively). At time $`t_2`$ \[Eq. (15)\], the uphill front reaches the upper limit of the cell ($`x_{uh}=L`$). The static profile is then in its relaxed final state, with a uniform slope $`\theta _n`$. This is the end of Stage II, which is thus defined as the time interval $`t_1<tt_2`$. In most cases, as is discussed below, we expect the downhill front to reach the bottom edge before $`t=t_2`$.

#### 1 Top region

In this region, we have $`x>x_{uh}`$, so that the evolution equation (17) still holds, and we still have to solve with respect to the upper boundary condition of Eq. (11). Therefore, as in Stage I, $`R`$ is given by Eq. (18) \[but now, the lower limit of the domain on which this solution is valid is $`x_{uh}(t)`$, not $`x_{dh}(t)`$\]. This top region shrinks, until finally disappearing when the uphill front reaches the upper end of the cell ($`t=t_2`$).

#### 2 Central region

In this part, since $`xx_{uh}(t)`$, there is no amplification of the rolling amount, so that the right-hand side of the evolution equation of $`R`$ vanishes: $`R/t\mathrm{\Gamma }(R+d)R/x=0`$. Now, we further impose that $`R`$ shall be continuous at the border with the top region, i.e.:

$`R(x=x_{uh}(t),t)=d+\sqrt{d^2+2(Lx){\displaystyle \frac{\eta v_{uh}}{\mathrm{\Gamma }}}}.`$

Solving these two equations together leads to look for $`R`$ as one of the roots of the third-degree equation $`R^3+a_2R^2+a_1(t)R+a_0(x,t)=0`$, with the following coefficients:

$`a_2`$ $`=`$ $`\left({\displaystyle \frac{v_{uh}}{\mathrm{\Gamma }}}+3d\right)`$   
$`a_1(t)`$ $`=`$ $`2{\displaystyle \frac{v_{uh}}{\mathrm{\Gamma }}}\left(d+{\displaystyle \frac{\mathrm{\Gamma }d^2}{v_{uh}}}\eta (Lv_{uh}t)\right)`$   
$`a_0(x,t)`$ $`=`$ $`2\eta {\displaystyle \frac{v_{uh}}{\mathrm{\Gamma }}}\left({\displaystyle \frac{v_{uh}}{\mathrm{\Gamma }}}(Lx)+2d(Lv_{uh}t)\right).`$   
The solution, given by Cardano formulas \[\], writes:

$$R(x,t)=\frac{a_2}{3}+S\frac{Q}{S},$$
(24)   
with the auxiliary quantities \[\]:

$`S`$ $``$ $`\sqrt[3]{P+\sqrt{D}}`$   
$`D`$ $``$ $`Q^3+P^2`$   
$`Q`$ $``$ $`{\displaystyle \frac{1}{9}}(3a_1a_2^2)`$   
$`P`$ $``$ $`{\displaystyle \frac{9a_2a_127a_02a_2^3}{54}}.`$   
We saw in the previous section that the crest of the avalanche $`R_{max}`$ appeared at the end of Stage I. What happens to this crest during Stage II? It is easy to prove that the crest remains located on the downhill front $`x_{dh}`$. Besides, $`x_{dh}`$ now moves at constant speed (in contrast with Stage I where it accelerated): $`x_{dh}(t)=v_{uh}t_1\mathrm{\Gamma }(R_{max}+d)(tt_1)`$.

We can also prove that the height of the crest remains constant (equal to $`R_{max}`$) as it travels downhill, until it finally comes out of the cell. The exit time $`t_{exit}`$ of this crest $`R_{max}`$ is obtained by solving $`x_{dh}(t)=0`$:

$$t_{exit}=t_1+\frac{v_{uh}t_1}{\mathrm{\Gamma }(R_{max}+d)}=t_1\left(1+\frac{v_{uh}}{\mathrm{\Gamma }(R_{max}+d)}\right).$$
(25)   
#### 3 Bottom region

The bottom region is defined as the region where $`0xx_{dh}(t)`$. The evolution equation for $`R`$ is the same as in the central region, and we impose continuity at the border with the central region. We find that $`R(x,t)`$ is given by Eq. (22) as in Stage I. Physically, in this region, we simply observe the convection of what was left in the bottom region at the end of Stage I.

The bottom region disappears when the downhill front reaches the bottom end: $`x_{dh}(t)=0`$, that is by definition at time $`t=t_{exit}`$ \[Eq. (25)\]. To determine precisely the subsequent evolution of the avalanche, we must discuss whether the disappearance of the bottom region occurs before the end of Stage II or not (that is, whether $`t_{exit}t_2`$ or $`t_{exit}t_2`$). Using Eq. (25), we form the ratio:

$`{\displaystyle \frac{t_{exit}}{t_2}}={\displaystyle \frac{t_1}{t_2}}\left(1+{\displaystyle \frac{v_{uh}}{\mathrm{\Gamma }(R_{max}+d)}}\right).`$

Provided that the cell dimension $`L100d`$ and that $`v_{uh}\mathrm{\Gamma }d`$ (these requirements being usually satisfied for common experiments), we have $`\mathrm{\Gamma }R_{max}v_{uh}`$, and $`t_1/t_21`$. Hence we generally expect $`t_{exit}t_2`$, and, consequently, as claimed earlier, in most cases the bottom region disappears before the end of Stage II.

To resume, during Stage II, the central region extends both downhill and uphill with constant (though different) velocities at each end, progressively invading the whole cell. It reaches the bottom edge at time $`t_{exit}`$ (in situations where $`t_{exit}t_2`$), then the upper edge at time $`t_2`$. At the end of Stage II, the central region occupies the entire cell, and $`R(x,t=t_2)`$ is everywhere given by Eq. (24). We present successive calculated ‘snapshots’ of the rolling amount $`R`$ during this stage in Fig. 5-b.

### C Stage III: The last grains are evacuated

This stage lasts from $`t=t_2`$ until the end of the avalanche, at $`t=t_{end}`$. Both fronts have reached the edges of the cell at the end of Stage II ($`x_{dh}=0,x_{uh}=L`$), and there is only one region (see Fig. 6-a). Moreover, as $`t>t_2`$, the slope of the static part is everywhere $`\theta _n`$ and no amplification of the rolling grains can take place; the rolling phase is simply convected downwards. We now have to solve the evolution equation $`R/t\mathrm{\Gamma }(R+d)R/x=0`$, with respect to the initial condition that Stage III evolves from what has been left by Stage II \[i.e. $`R(x,t=t_2)`$ as given by Eq. (24)\].

Solving with the method of characteristics gives the following implicit solution

$$R(x,t)=R(\xi ,t_2),$$
(26)   
where

$$\xi =x\mathrm{\Gamma }\left[R(\xi ,t_2)+d\right](tt_2)$$
(27)   
(and $`0\xi L`$).

The physical interpretation of these equations is actually very simple: Eq. (26) states that the quantity of rolling species found in $`x`$ at time $`t`$ was previously located in $`\xi `$ at the beginning of Stage II $`(t=t_2)`$. Equation (27) gives a determination of this initial position $`\xi `$, by stating that from $`t_2`$ until the considered instant $`t`$, the quantity of grains moved with a constant speed $`\mathrm{\Gamma }[R(\xi ,t_2)+d]`$, dependent on the local height. In other words, during Stage III, the $`R`$-profile left by Stage II is convected downhill, but each vertical slice rolls with its own velocity, which is a function of its height. The grains that were near the top edge of the cell at the end of Stage II are convected the most slowly, since there $`R`$ was close to zero. The profile inherited from Stage II thus dilates upon rolling, under the effect of velocity inhomogeneities (Fig. 6-b) \[\].

The last grains to fall off the cell are those that leave the top end of the cell at the beginning of Stage III, at time $`t_2`$. Since at the top edge we have $`R=0`$, these grains move with a constant speed $`v_0=\mathrm{\Gamma }d`$. At time $`t`$, they are located at $`x_{last}(t)=L\mathrm{\Gamma }d(tt_2)`$, and the avalanche is extinct uphill: $`R=0`$ for $`x>x_{last}(t)`$.

Finally, the avalanche ends when the last grains reach the bottom limit of the cell ($`x_{last}=0`$), that is at time $`t_{end}=t_2+L/(\mathrm{\Gamma }d)`$.

## IV Discussion and simple checks

### A Predictions for the maximum amplitude of the avalanche

#### 1 Linear velocity profile

Up to now, we have focused on flows displaying linear velocity profiles. For such flows, as we saw in Section III A, the avalanche reaches its maximum amplitude $`R_{max}`$ at the end of Stage I, at time $`t=t_1`$ \[Eq. (23)\]. The exact analytical expression of $`R_{max}`$ is easily found by using Eq. (20) at time $`t_1`$ (that is, the value of $`R`$ given by the central region at the very moment it disappears):

$$R_{max}=d\frac{v_{uh}}{\mathrm{\Gamma }}+\sqrt{\left(d+\frac{v_{uh}}{\mathrm{\Gamma }}\right)^2+\frac{2L\eta v_{uh}}{\mathrm{\Gamma }}}.$$
(28)   
For large values of L, $`R_{max}`$ scales as:

$$R_{max}\sqrt{2\eta \frac{v_{uh}}{\mathrm{\Gamma }}L},$$
(29)   
that is, as the square root of the system size $`L`$.

Let us give a couple of numerical applications of this last expression. For the case of a standard laboratory experiment, with $`L=1`$ m, $`d=1`$ mm, $`v_{uh}/\mathrm{\Gamma }=3d`$ and $`\eta 0.1`$ rad, we find $`R_{max}=2.45`$ cm. In the case of a system at the scale of a desert dune, made of fine sand, we take $`L=10`$ m, $`d=0.2`$ mm and, with others parameters unchanged, we get $`R_{max}=3.46`$ cm. One has to notice that $`R_{max}`$ is quite small, even for large systems as a sand dune.

It is interesting to contrast this result with the work of Boutreux et al. \[\], who carried the same calculation in an open cell configuration, but with a constant downhill convection velocity $`v(R)v_0`$ \[instead of Eq. (4)\]. They found $`R_{max}\eta L`$. For the two above examples, this formula leads to maximum amplitudes of respectively $`10`$ cm and $`1`$ m. The effect of the velocity gradient is thus to considerably limit the amplitude of avalanches, especially for large systems.

#### 2 Generalization for a power-law dependency

In the beginning of this article, we quoted the work of Azanza et al. \[\] and of Pouliquen \[\] who find that the average speed of a chute flow of granular material on a rough plane is related to its thickness through a power-law relation $`v(R)\mathrm{\Gamma }R^\alpha `$ with $`\alpha `$ close to 3/2. However, as pointed out by Pouliquen \[\], the influence of the rough underlying bed plane on the rheology of chute flows is complex and not clearly understood, and might not be comparable to situations where the flow occurs on a free granular bed as has been considered in this paper. Since the question is still open, we will here present an intuitive derivation of $`R_{max}`$ valid for any power-law (undetermined exponent $`\alpha `$). To check the validity of this simple derivation, we first present it in the linear case $`\alpha =1`$, the generalization being then straightforward.

Let us consider a point initially at the top edge of the cell. At $`t=0`$, it starts being swept along by the granular flow and we assume that this point travels with the local surface velocity of the flow $`v=\mathrm{\Gamma }R`$. We are now interested in the temporal evolution of the rolling height $`R`$ at this travelling point, which shall be computed from the Lagrangian derivative $`dR/dt=R/t+vR/x`$. As long as the amplification process takes place, we have with the use of Eq. (16): $`dR/dt=\eta v_{uh}`$. This implies

$$R(t)=\eta v_{uh}t.$$
(30)   
Hence, $`R(t)`$ at the travelling point increases with time, as long as the amplification process lasts. After the amplification has stopped, $`R`$ at the travelling point keeps constant (since $`dR/dt=0`$). Thus, $`R`$ reaches its maximum value $`R_{max}`$ at the end of the amplification. Let us call $`t_{amp}`$ the duration of the amplification. We compute $`t_{amp}`$ in the following way: the distance that the travelling point goes over during the amplification is of order $`L`$, so that $`t_{amp}`$ must verify:

$`L`$ $``$ $`{\displaystyle _0^{t_{amp}}}𝑑x={\displaystyle _0^{t_{amp}}}v𝑑t`$   
$`\text{i.e.}L`$ $``$ $`{\displaystyle _0^{t_{amp}}}\mathrm{\Gamma }\eta v_{uh}t𝑑t={\displaystyle \frac{1}{2}}\mathrm{\Gamma }\eta v_{uh}t_{amp}^{}{}_{}{}^{2},`$   
\[in this calculation, we used $`v\mathrm{\Gamma }R(t)`$\]. We finally find: $`t_{amp}\sqrt{2L/(\mathrm{\Gamma }\eta v_{uh})}`$.

Inserting this last expression into Eq. (30) gives the value of $`R_{max}`$:

$$R_{max}\sqrt{2\eta \frac{v_{uh}}{\mathrm{\Gamma }}L}.$$
(31)   
This is exactly Eq. (29) found analytically, which was itself the limit of the complete expression of $`R_{max}`$ \[Eq. (28)\] for large values of L (greater than a hundred $`d`$), and $`v_{uh}`$ of order $`\mathrm{\Gamma }d`$.

The strongest assumption in the above simple derivation is that the amplification takes place over a distance $`L`$. Rigorously, this distance is $`Lx_{dh}(t=t_1)`$; but what makes our simple derivation successful is that the position of the downhill front at time $`t_1`$, $`x_{dh}(t=t_1)`$, is generally quite close to zero (for $`L`$ greater than a hundred $`d`$ and $`v_{uh}`$ of order $`\mathrm{\Gamma }d`$).

We may now generalize the above results to a power-law dependency of the velocity $`v(R)\mathrm{\Gamma }R^\alpha `$. The same derivation leads us to the result:

$$R_{max}\left((\alpha +1)\frac{\eta v_{uh}}{\mathrm{\Gamma }}L\right)^{1/(\alpha +1)}$$
(32)   
Note that $`R_{max}`$ diminishes as $`\alpha `$ increases.

In particular, for $`\alpha =3/2`$, Eq. (32) can be rewritten as $`R_{max}(5\eta Lv_{uh}/\mathrm{\hspace{0.17em}2}\mathrm{\Gamma })^{\mathrm{\hspace{0.25em}2}/5}`$.

### B Possible experimental checks

The loss of material at the bottom edge of the cell might be measured experimentally, and could be compared to the following theoretical prediction. This loss corresponds to the flow rate at the bottom of the cell $`Q(x=0,t)=_0^{R(x=0,t)}v(z)𝑑z`$, and is given by

$$Q(x=0,t)=\frac{\mathrm{\Gamma }}{2}R(x=0,t)^2+\mathrm{\Gamma }dR(x=0,t),$$
(33)   
where $`R(x=0,t)`$ is given by Eq. (22) during Stages I and II, and by Eq. (26) during Stage III. Figure 7 shows the predicted shape of $`Q(x=0,t)`$ as a function of time (solid curve). The curve displays a maximum at time $`t=t_{exit}`$, corresponding to the moment when the maximum amplitude $`R_{max}`$ rolls out of the cell. The maximum flow rate is obtained by replacing $`R(x=0,t)`$ by $`R_{max}`$ in Eq. (33).

It is of interest to compare our prediction for the loss of material with that of Boutreux et al. \[\], who assumed a constant downhill velocity $`v`$ in the rolling phase. This comparison, however, requires some caution: in our approach, the granular flow is characterized by a constant velocity gradient $`\mathrm{\Gamma }`$, whereas, in Boutreux et al., the description is based on a typical downhill convection velocity of the grains $`v`$ (see Section I C). Figure 7 compares the results of both approaches for the loss of material, assuming $`vv_{uh}3\mathrm{\Gamma }d`$ (see Ref. \[\]).

### C Concluding remarks

#### 1 Regions of small $`R`$

We notice that, during the avalanche, we had several spatial zones in the cell where R was close to 0 ($`R<\xi ^{}`$), e.g. at the upper edge of the cell, or at the end of the avalanche. Thus in these zones, the use of the saturated Eqs. (7) and (8) is not fully justified. In order to obtain a continuous description between the saturated case and the thin one, we could use the interpolated equations that have been proposed by de Gennes \[\] and studied in a model case by Boutreux and Raphaël in Ref. \[\]:

$`{\displaystyle \frac{\stackrel{~}{h}}{t}}`$ $`=`$ $`\gamma {\displaystyle \frac{R\xi ^{}}{R+\xi ^{}}}{\displaystyle \frac{\stackrel{~}{h}}{x}}`$   
$`{\displaystyle \frac{R}{t}}`$ $`=`$ $`\mathrm{\Gamma }(R+d){\displaystyle \frac{R}{x}}+\gamma {\displaystyle \frac{R\xi ^{}}{R+\xi ^{}}}{\displaystyle \frac{\stackrel{~}{h}}{x}}.`$   
The results of \[\] show however that the physical behaviour is not dramatically changed, and that the description in the zones of small R with saturated equations might be slightly wrong but qualitatively verified.

#### 2 Effects of polydispersity

It is of common knowledge that real granular materials are generally intrinsically polydisperse. This may have drastic effects on the behavior of the flow, and capturing more precisely the physics of real avalanches would certainly suppose to take polydispersity into account. However, the treatment of full polydispersity is a difficult task. Yet, the BCRE equations have been extended to the case of binary mixtures \[\], and it could be interesting to study the changes brought up in this case by a velocity gradient in the flow.

#### 3 Domain of validity of the BCRE approach

The general approach introduced by Bouchaud et al. to describe surface flows is rather phenomenological, and as pointed out by Bouchaud and Cates \[\], we still lack criteria to determine the range of physical situations to which it can be successfully applied. In a recent work, Douady et al. \[\] proposed a justification of the BCRE modelization on the basis of hydrodynamic conservation laws. According to these authors, the form of the BCRE equations should not remain invariant when different laws are chosen for the velocity profile in the flow. Douady et al. argue that only for a velocity linear in $`R`$ (or constant) shall the equations take the simple form of our equations \[Eqs. (7) and (8)\]; in other cases, they find that a supplementary term coupling $`R`$ and $`h`$ should add in Eq. (7). Certainly, more work needs to be done in this direction in order to exactly assert the domain of validity of the BCRE analysis.

###### Acknowledgements.

We thank T. Boutreux, F. Chevoir, A. Daerr, S. Douady, J. Duran and O. Pouliquen for oral and/or written exchanges.      

## 

### 1 Method for solving the evolution equation of $`R`$

Eq. (16) is a first-order partial differential equation, of the quasi-linear class, that is, linear in the first derivatives. Such equations can be solved by the well-known method of characteristics. See, for example, Ref. \[\].

More specifically, we will solve Eq. (16) along characteristic curves given in the parametric form $`\{t(s),x(s),R(s)\}`$, with $`s`$ the parameter. The functions $`t(s)`$, $`x(s)`$ and $`R(s)`$ are derived from the set of coupled ordinary differential equations:

$`{\displaystyle \frac{dt}{ds}}`$ $`=`$ $`1`$ (34)   
$`{\displaystyle \frac{dx}{ds}}`$ $`=`$ $`\mathrm{\Gamma }(R+d)`$ (35)   
$`{\displaystyle \frac{dR}{ds}}`$ $`=`$ $`\eta v_{uh}H(xv_{uh}t).`$ (36)   
By integration, one founds the equations for the characteristics with unspecified integration constants. One then imposes the boundary and/or initial conditions to identify these constants.

We here give the detailed calculations only for the first stage of the avalanche. The derivations are separated into the different spatial regions that were defined earlier, and we will show how they naturally emerge from the derivations.

### 2 Top region

In this region, Eqs. (35) become:

$`{\displaystyle \frac{dt}{ds}}`$ $`=`$ $`1`$ (37)   
$`{\displaystyle \frac{dx}{ds}}`$ $`=`$ $`\mathrm{\Gamma }(R+d)`$ (38)   
$`{\displaystyle \frac{dR}{ds}}`$ $`=`$ $`\eta v_{uh}.`$ (39)   
We also use the boundary condition $`R(x=L,t)=0`$ for $`t0`$, which we parametrize with the parameter $`\xi `$. For simplicity’s sake, on each characteristic crossing the boundary curve we arbitrarily choose the value of the parameter $`s`$ to be zero at the crossing point. This determines the integration constants to be:

$`t(s=0)`$ $`=`$ $`\xi `$ (40)   
$`x(s=0)`$ $`=`$ $`L`$ (41)   
$`R(s=0)`$ $`=`$ $`0.`$ (42)   
Note that $`\xi 0`$, since $`t0`$ (the experiment started at time $`t=0`$ on).

Solving for Eqs. (38) together with Eqs. (41) gives the equations of the characteristic curves:

$`t(s)`$ $`=`$ $`s+\xi `$ (43)   
$`x(s)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Gamma }\eta v_{uh}s^2\mathrm{\Gamma }ds+L`$ (44)   
$`R(s)`$ $`=`$ $`\eta v_{uh}s.`$ (45)   
We now want to write the solution $`R`$ explicitly in terms of $`x`$ and $`t`$, so that we have to eliminate $`\xi `$ and $`s`$. Eq. (44) can be solved to give $`s`$ as a function of $`x`$, and replacing into Eq. (45) brings the analytical solution:

$$R(x,t)=d+\sqrt{d^2+2(Lx)\frac{\eta v_{uh}}{\mathrm{\Gamma }}},$$
(46)   
which is the same as Eq. (18).

We now have to verify the condition that $`\xi 0`$. By combining Eq. (43) with (45) this condition can be rewritten as $`tR(x,t)/\eta v_{uh}`$. Replacing $`R`$ into this inequality \[by (46)\] gives us a spatial condition for solution (46) to be valid: we must have $`xx_{dh}(t)`$, where

$`x_{dh}(t)L+{\displaystyle \frac{\mathrm{\Gamma }}{2\eta v_{uh}}}d^2{\displaystyle \frac{1}{2}}\mathrm{\Gamma }\eta v_{uh}(t+{\displaystyle \frac{d}{\eta v_{uh}}})^2.`$

This is the mathematical origin of the ‘downhill front’ that we described intuitively in the main text as the limit of extension of the boundary effects originating in the upper edge of the cell.

### 3 Central region

The evolution equation for $`R`$ is the same as in the top zone, so that the differential equations giving the characteristics also are the same \[Eqs. (38)\]. But now we must impose the initial condition $`R(x,t=0)=0`$, which gives the following set of initial conditions for the characteristics: $`t(s=0)=0,x(s=0)=\xi ,R(s=0)=0`$ (and $`0\xi L`$). We obtain:

$`t(s)`$ $`=`$ $`s`$ (47)   
$`x(s)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Gamma }\eta v_{uh}s^2\mathrm{\Gamma }ds+L`$ (48)   
$`R(s)`$ $`=`$ $`\eta v_{uh}s.`$ (49)   
Combining (47) and (49) gives an explicit solution for $`R`$: $`R(x,t)=\eta v_{uh}t`$.

This solution is valid in a certain spatial domain. It is limited upwards by the top region \[i.e. $`xx_{dh}(t)`$\]. It is also limited downwards by $`x_{uh}(t)`$, because at this point the form of the evolution equation of $`R`$ changes (the amplification term vanishes), and consequently does the form of the differential equations that give the characteristics.

### 4 Bottom region

In this region, Eqs. (35) are given by: $`dt/ds=1,dx/ds=\mathrm{\Gamma }(R+d),dR/ds=0`$.

Here, the boundary condition is given by the continuity of $`R`$ at the border of the central and the bottom zones: $`R(x=x_{uh}(t),t)=\eta v_{uh}t`$ for $`t0`$. This gives the initial conditions: $`t(s=0)=\xi ,x(s=0)=v_{uh}\xi ,R(s=0)=\eta v_{uh}\xi `$. Solving and rewriting $`R`$ explicitly in terms of $`x`$ and $`t`$ leads to the solution:

$`R(x,t)={\displaystyle \frac{1}{2}}\left(d+{\displaystyle \frac{v_{uh}}{\mathrm{\Gamma }}}\eta v_{uh}t\right)`$ (51)   
$`+{\displaystyle \frac{v_{uh}}{2\mathrm{\Gamma }}}\sqrt{\left({\displaystyle \frac{\mathrm{\Gamma }d}{v_{uh}}}+1\mathrm{\Gamma }\eta t\right)^2+4{\displaystyle \frac{\mathrm{\Gamma }\eta }{v_{uh}}}(x+\mathrm{\Gamma }dt)},`$   
valid for $`0xx_{uh}(t)`$.