# Exact Solutions of a Model for Granular Avalanches

## 1 Introduction and Model

The study of avalanches and surface flows in granular materials has attracted much attention recently, both from a theoretical  and an experimental point of view . A simple model, thought to capture some of the essential phenomena, has been proposed in . It is based on the assumption that a strict separation between rolling grains and static grains can be made. Coupled dynamical equations for these two species, based on phenomenological arguments, can then be written. Calling $`R`$ the local density of rolling grains and $`h`$ the height of static grains, the simplest form of the bcre equations read:

$`H_t`$ $`=`$ $`\gamma RH_x,`$ (1)   
$`R_t`$ $`=`$ $`R_x+\gamma RH_x,`$ (2)   
where $`H`$ is the height of static grains, counted from the repose slope of angle $`\theta _r`$: $`h(t,x)=H(t,x)+x\mathrm{tan}(\theta _r)`$ (the heap is sloping upwards from left to right). In the above equations, the units of lengths and time are chosen such that the (downhill) velocity of grains is $`v=1`$, while $`H`$ and $`R`$ are counted in units of the grains diameter. The term $`\gamma RH_x`$ describes the conversion of static grains into rolling grains if $`H_x>0`$, or vice versa if $`H_x<0`$. $`\gamma `$ is a grain collision frequency, typically of the order of $`100`$ Hz.

Many important phenomenon are left out from the above description, and can be included by adding more terms. For example, diffusion terms (such as $`D_1R_{xx}`$ or $`D_2H_{xx}`$, describing, e.g. non local dislodgement effects) will generically be present, and qualitatively change the structure of the solutions . Another aspect not described by the linear form of the conversion term above is the expected saturation of rolling grains with time, rather than the exponential growth predicted by Eq. (2) for a constant positive slope $`H_x`$. Non linear saturation terms, as well as a dependence of the velocity of the rolling grains on $`R`$, are thus expected in general, and can lead to important differences with the above equations .

Recently, these equations has been studied by Mahadevan and Pomeau (mp) . They found a conservation law, which relates the solutions $`R(t,x)`$ and $`H(t,x)`$ in a frame moving with the velocity of the grains. From this law, they concluded that the bcre equations have characteristics that are straight lines, along which both $`R(t,x)`$ and $`H(t,x)`$ are constant. Independently of the initial profile $`H_0(x)`$, they found that a shock forms at time $`t_s=1/(\gamma R_{0,\mathrm{max}}^{})`$ with $`R_{0,\mathrm{max}}^{}`$ is the maximum (in absolute value) of the initial gradient of rolling grains. Whereas our exact solution fulfills the same conservation law, our results for the characteristics and the shock time disagree with the results of mp. As we will discuss below, the reason for this disagreement is their implicit assumption of a very restrictive relation between the initial profiles $`R_0(x)`$ and $`H_0(x)`$.

## 2 Characteristic coordinates

The general basis of the method  we used to solve Eqs. (1,2) consists in a replacement of the original equations by an equivalent system of four partial differential equations for the functions $`t`$, $`x`$, $`R`$ and $`H`$, but now considered as functions of new coordinates $`\mu `$ and $`\nu `$, which will be defined below.<sup>1</sup><sup>1</sup>1The theory used here is actually more general and can be used in the presence of non-linear saturation terms or for ripple models . These new equations will be particularly simple inasmuch as each equation has derivatives with respect to either $`\mu `$ or $`\nu `$, though the mapping between the coordinates $`(t,x)`$ and $`(\mu ,\nu )`$ will be in general complicated. To define the characteristic coordinates $`(\mu ,\nu )`$, we have to specify first the characteristic curves of the system (1,2). For practical reasons, we introduce new functions $`u(t,x)=1R(t,x)/\alpha `$ and $`v(t,x)=(\alpha +xR(t,x)H(t,x))/\alpha `$ instead of $`H(t,x)`$ and $`R(t,x)`$. For this new functions the differential expressions become

$`L_1[u,v]`$ $`=`$ $`u_t\gamma \alpha (1u)u_x+v_t+\gamma \alpha (1u)v_x\gamma (1u)=0,`$ (3)   
$`L_2[u,v]`$ $`=`$ $`u_t+[1+\gamma \alpha (1u)]u_x\gamma \alpha (1u)v_x+\gamma (1u)=0.`$ (4)   
Both operators $`L_1`$ and $`L_2`$ contain linear combinations of the type $`au_t+bu_x`$ of the derivatives of $`u`$ (and the same holds for $`v`$). This combination means that $`u`$ is differentiated in the direction given by the ratio $`t/x=a/b`$. Since the coefficients $`a`$ and $`b`$ differ for $`u`$ and $`v`$ and also for $`L_1`$ and $`L_2`$, the functions $`u`$ and $`v`$ are differentiated in each of the operators in different directions in the $`(t,x)`$ plane. Notice that the directions depend also on $`u`$ itself, and therefore on the solution under consideration, which is a typical feature of non-linear systems. As noted above, our goal is to find equivalent differential equations of which each contains derivatives in only one (local) direction corresponding to one of the new coordinates $`\mu `$ and $`\nu `$. Therefore we take a linear combination $`L=\lambda _1L_1+\lambda _2L_2`$ of the operators in Eqs. (3,4) such that the derivatives of $`u`$ and $`v`$ in $`L`$ combine to derivatives in the same direction, which is called a characteristic direction. Moreover we assume that these local directions change smoothly as functions of $`t`$ and $`x`$, and are given by the tangential vectors $`(t_\sigma (\sigma ),x_\sigma (\sigma ))`$ of a smooth path $`(t(\sigma ),x(\sigma ))`$ with $`\sigma `$ as parameter. Considering the functions $`u`$ and $`v`$ along this path, they depend only on $`\sigma `$ and we have, e.g., $`u_\sigma =u_tt_\sigma +u_xx_\sigma `$. Using these conditions, we obtain four homogeneous linear equations for the coefficients $`\lambda _1`$ and $`\lambda _2`$ with coefficients depending on $`t`$, $`x`$, $`u`$, $`v`$ and their derivatives with respect to $`\sigma `$. For non-trivial solutions all possible determinants of the matrix of these coefficients have to vanish, leading to three independent equations or characteristic relations (cr). The first one can be written as a quadratic equation for the local direction $`\zeta =x_\sigma /t_\sigma `$ of differentiation, the solution of which are: $`\zeta _+=1`$ and $`\zeta _{}=\gamma \alpha (1u)`$. Now, for a fixed solution $`u`$, the equations $`dx/dt=\zeta _+`$ and $`dx/dt=\zeta _{}`$ are ordinary differential equations, which define two families of paths with the starting position $`x_0`$ at $`t=0`$ as parameter. These families of paths are the characteristics $`C_+`$ and $`C_{}`$ of the system (1,2). From a physical point of view, they are simply the paths along which $`R(x,t)`$ ($`\zeta _+`$) and $`H(t,x)`$ ($`\zeta _{}`$) evolves with time.

The new curved coordinate frame $`(\mu ,\nu )`$ is now defined such that the two one-parametric families of characteristics are mapped by the coordinate transformation on an usual Cartesian coordinate frame in the $`(\mu ,\nu )`$-plane, i.e., along the characteristics the coordinate functions $`\mu (t,x)`$ and $`\nu (t,x)`$, respectively, are constant. Here we have chosen to map the line $`t=0`$ on the line given by $`\mu =\nu `$. In terms of the new coordinates we find

$$x_\nu +t_\nu =0,x_\mu \gamma \alpha (1u)t_\mu =0.$$
(5)   
Now we make use of another cr, which evaluated along $`C_+`$ and $`C_{}`$ by identifying $`\sigma `$ with $`\nu `$ and $`\mu `$, respectively, yields the conditions

$$u_\nu +\gamma \alpha (1u)v_\nu +\gamma (1u)t_\nu =0,u_\mu v_\mu +\gamma (1u)t_\mu =0.$$
(6)   
These equations together with the Eqs. (5) form the desired set of four equations mentioned before. Every solution of this new system satisfies the original Eqs. (1,2), since the Jacobian $`t_\nu x_\mu t_\mu x_\nu 1+\gamma R(\mu ,\nu )`$ of the coordinate map does not vanish due to $`\gamma R(\mu ,\nu )>0`$.

## 3 General solution

Before we can construct a solution to the equivalent system (5,6), we have to specify initial data along the line $`\mu =\nu `$ corresponding to $`t=0`$. We choose an general profile $`H_0(x)`$, perturbed at $`t=0`$ by a uniform ‘rain’ of rolling grains: $`R_0(x)=\alpha `$. In terms of the new coordinates, the initial conditions become $`t_0(\mu )=0`$, $`x_0(\mu )=\mu `$, $`u_0(\mu )=0`$, $`v_0(\mu )=(\mu +H_0(\mu ))/\alpha `$. By introducing the function $`\mathrm{\Delta }(\mu ,\nu )=1\gamma \alpha (1u(\mu ,\nu ))`$, one can show that the problem of solving the system given by Eqs. (5,6) can be reduced to the task of finding a solution to the equation $`\mathrm{\Delta }_\nu =\gamma H_0^{}(\nu )(1+1/\mathrm{\Delta })`$, with initial condition $`\mathrm{\Delta }(\mu ,\mu )=1\gamma \alpha `$. The solution of this equation can be simply expressed in terms of the so-called Lambert function $`W`$ :

$$\mathrm{\Delta }(\mu ,\nu )=1W\left\{\alpha \gamma \mathrm{exp}[\alpha \gamma +\gamma (H_0(\mu )H_0(\nu ))]\right\}.$$
(7)   
With this solution at hand, the solution to the system (5,6) is determined by

$`t(\mu ,\nu )={\displaystyle _\nu ^\mu }{\displaystyle \frac{ds}{\mathrm{\Delta }(s,\nu )}}=\mu \nu +{\displaystyle _\nu ^\mu }{\displaystyle \frac{\mathrm{\Delta }_\mu (s,\nu )ds}{\gamma H_0^{}(s)}}`$ , $`x(\mu ,\nu )=\mu t(\mu ,\nu )`$   
$`R(\mu ,\nu )={\displaystyle \frac{1+\mathrm{\Delta }(\mu ,\nu )}{\gamma }}`$ , $`H(\mu ,\nu )=H_0(\nu ),`$ (8)   
where we have expressed already the original fields $`R(\mu ,\nu )`$ and $`H(\mu ,\nu )`$ in terms of the functions $`u`$ and $`v`$. To get the fields as functions of $`t`$ and $`x`$, one has to invert the coordinate map. This can be done by using $`\mu (t,x)=tx`$ and integrating the equation for $`t(\mu ,\nu )`$ to obtain also $`\nu (t,x)`$. As announced before, the height profile $`H(t,x)=H_0(\nu (t,x))`$ turns out to be constant along the characteristics $`C_{}`$.

## 4 Generic shape for $`H(t,x)`$

In the following we will consider a situation which is generic for sandpile surfaces. Suppose that one starts with a sandpile profile, which consists of two regions with constant but different slopes matching with a kink at $`x=0`$, and again with a constant amount of rolling grains. The slopes may be either larger or smaller than the angle of repose $`\theta _r`$. If we denote the slope to the right (left) by $`\theta _r+\theta _>`$ ($`\theta _r+\theta _<`$), we have $`H_0(x)=\theta _>x`$ for $`x>0`$ and $`H_0(x)=\theta _<x`$ for $`x<0`$. In the case of a piecewise constant $`H_0^{}(x)`$ one can integrate the equation for $`t(\mu ,\nu )`$ easily as can be seen from Eq. (8). The structure of Eq. (7) suggests to distinguish between three regions given by $`\mu >0`$,$`\nu <0`$ (I), $`\mu ,\nu <0`$ (II) and $`\mu <0,\nu >0`$ (III). <sup>2</sup><sup>2</sup>2The region where $`\mu ,\nu >0`$ turns out to be mapped on the half space with $`t<0`$ and is therefore not of physical interest. In regions I and III one can find the explicit expression $`\nu (t,x)=x+\frac{\alpha }{\theta }(1e^{\gamma \theta t})`$ with $`\theta =\theta _<`$ (I) or $`\theta =\theta _>`$ (III), i.e., the characteristics $`C_{}`$ are in these regions simple exponential curves. As a consequence, no shocks can appear in these two regions and the corresponding solutions are particularly simple:

$$R_{I(III)}(t,x)=\alpha e^{\gamma \theta _{<(>)}t},H_{I(III)}(t,x)=H_0(x)+\alpha R_{I(III)}(t,x).$$
(9)   
The boundaries of the regions I and III in real space $`(t,x)`$ are given by the conditions $`x<t`$ and $`x>x_1(t)=\frac{\alpha }{\theta _>}(e^{\gamma \theta _>t}1)`$ corresponding to the $`\mu =0`$ and $`\nu =0`$ characteristics, respectively, see Fig. 1. The boundary for region I has an obvious physical meaning: The information that there is a kink at $`x=0`$ can only propagate to the left with the velocity of the moving grains, which is $`1`$ in our rescaled units. Moreover, it is important to note that the ‘uphill’ velocity with which the kink moves is only equal to $`\gamma \alpha `$ at small times, before growing exponentially. As discussed in the introduction, this growth eventually saturates, as does the value $`R`$, or else the characteristic $`C_{}`$ quickly reaches the edge of the pile.

The range of $`x`$ in between the above two regions corresponds to the intermediate region II. Within this range one can obtain only an implicit solution for the coordinate map $`\nu (t,x)`$. It reads

$$\nu =x+\frac{1}{\gamma }\left[\frac{\mathrm{\Delta }(xt,\nu )\mathrm{\Delta }(0,\nu )}{\theta _>}+\frac{\mathrm{\Delta }(0,\nu )+1+\alpha \gamma }{\theta _<}\right],$$
(10)   
where $`\mathrm{\Delta }(\mu ,\nu )=1W\left\{\alpha \gamma \mathrm{exp}[\alpha \gamma \gamma (\theta _>\mu +\theta _<\nu )]\right\}`$ as follows from Eq. (7). The shape of $`R(t,x)`$ and $`H(t,x)`$ can be obtained directly from the last two equations of (8). In general, Eq. (10) has to be solved numerically although several results can be obtained in an analytic way. It turns out that the solutions of Eq. (10) fall into two qualitatively different classes, according to the values of $`\beta =\theta _>/\theta _<`$ and $`\theta _<`$: for $`\beta >1\alpha \gamma `$ or $`\theta _<<0`$, both $`R(t,x)`$ and $`H(t,x)`$ remain continuous for all times, while for $`\beta <1\alpha \gamma `$ and $`\theta _<>0`$, the solutions develop a discontinuity in $`R(t,x)`$ and $`H(t,x)`$ beyond a finite shock time $`t_s`$. This must be contrasted with mp, since in the present case $`R_0(x)=\alpha `$, they predict that shocks are absent for all times.

## 5 Examples

The characteristics resulting from numerical solutions of Eq. (10) have been plotted in Fig. 1. The left part of this figure has been obtained for $`\theta _>>0`$ and $`\theta _<<0`$, corresponding to $`\beta <0`$. In this case, the characteristics are more and more ‘diluted’ as time increases, and therefore never cross – no shock. In the limit of large times, the argument of the Lambert $`W`$ function becomes very large. Using the first two terms of the asymptotic expansion of $`W`$  we get $`\nu (t,x)=[\beta t+\mathrm{ln}(x+\frac{\beta t}{\beta 1})/(\gamma \theta _<)]/(\beta 1)`$. The corresponding expression for $`R(t,x)`$ and $`H(t,x)`$ can be obtained from Eq. (8). A particularly interesting quantity to look at is the local slope at, say, $`x=0`$. In this limit the slope is negative and decays with time as $`H_x(t,x=0)=1/(\gamma \beta t)`$.<sup>3</sup><sup>3</sup>3Note that this $`t^1`$ relaxation of the slope has also been obtained in  within a very different model. It means that the ‘true’ slope $`h_x`$ actually relaxes to the angle of repose $`\theta _r`$ for very large time. If $`L`$ is the size of experimental system, then $`C_{}`$ reaches the boundary of the system at a time $`t^{}`$ such that $`L\frac{\alpha }{\theta _>}e^{\gamma \theta _>t^{}}`$. One should therefore measure a final slope $`h_x\theta _r+\theta _</\mathrm{ln}(\theta _>L/\alpha )`$ smaller than the repose angle. This result is consistent with the qualitative discussion of Boutreux and de Gennes for a similar situation .

Another experimentally important quantity is the velocity $`v_R`$ of the “active” region. Following , this region can be defined by the condition $`R(t,x)>R_{\mathrm{min}}`$, where $`R_{\mathrm{min}}`$ is a small threshold. $`v_R`$ is then given by the slope of the curves of constant $`R(t,x)`$, which tends to a constant in the large $`t`$ limit as can be seen in Fig. 1(a). The asymptotic analysis yields $`v_R=\beta /(1\beta )`$. Since $`\beta <0`$, $`1<v_R<0`$, and the avalanche proceeds downhill, but slower than the grains themselves. This is an effect of the non-linear term in the bcre equations since the linearized theory yields $`v_R=1`$ .

The situation where $`\theta _><0`$ and $`\theta _<>0`$ is qualitatively different. In this case, the characteristics cross at some finite time: a shock occurs – see Fig. 1(b). A crossing point of two characteristics means indeed that at this point, two different values of $`R`$ (or $`H`$) are possible and these functions then become discontinuous. Strictly speaking, the Eqs. (1,2) are no longer valid, and the diffusion terms left of from the analysis become important to smooth out this discontinuity. In Fig. 2, we plotted snapshots of the $`h`$ and $`R`$ profiles at different times, for both situations (with and without the occurrence of a shock).

One can calculate the time $`t_s`$ and location $`x_s`$ at which the shock occurs. For that purpose, let us introduce the envelop of the characteristic curves $`x(t,\nu )`$, where $`\nu `$ is a label. The envelop can be represented in a parametric way as $`(t_e(\nu ),x_e(\nu ))`$. It has the property that for each of its points exists a characteristic, which touches it tangentially. It has then to fulfill the conditions $`x(t_e(\nu ),\nu )=x_e(\nu )`$, $`x_\nu (t_e(\nu ),\nu )=0`$. After some calculations, one can find the explicit expression for the envelop,

$`x_e(\nu )`$ $`=`$ $`\nu {\displaystyle \frac{1}{\gamma \theta _<}}\left[1+\alpha \gamma +\mathrm{\Delta }(0,\nu )\left(1+{\displaystyle \frac{1}{1\beta }}\right)\right]`$ (11)   
$`t_e(\nu )`$ $`=`$ $`x_e(\nu )+{\displaystyle \frac{\nu }{\beta }}{\displaystyle \frac{1}{\gamma \theta _>}}\left[\alpha \gamma +\mathrm{ln}(\alpha \gamma )+1+{\displaystyle \frac{\mathrm{\Delta }(0,\nu )}{1\beta }}\mathrm{ln}\left(1{\displaystyle \frac{\mathrm{\Delta }(0,\nu )}{1\beta }}\right)\right].`$ (12)   
This envelop has two branches, separated by a kink, see Fig. 1(b), given by $`\nu =\nu _s=[\alpha \gamma +\mathrm{ln}(\alpha \gamma )1+\beta \mathrm{ln}\left(1\beta \right)]/(\gamma \theta _<)`$. Whereas the upper branch is parameterized by $`\mathrm{}<\nu <\nu _s`$, the lower one corresponds to $`\nu _s<\nu <\nu _c=[\beta +\alpha \gamma +\mathrm{ln}(\alpha \gamma /\beta )]/(\gamma \theta _<)`$. The resulting shock coordinates are

$$x_s=\frac{1}{\gamma \theta _<}\left[\mathrm{ln}(\alpha \gamma )\mathrm{ln}(1\beta )+1+\frac{1}{1\beta }\right],t_s=\frac{1}{\gamma \theta _<}\left[\left(1\frac{2}{\beta }\right)\mathrm{ln}\left(1\beta \right)\mathrm{ln}(\alpha \gamma )\right].$$
(13)   
The condition that ($`t_s,x_s`$) has to be located inside region II leads to the boundary between the classes with and without shock as mentioned above. At the shock position, the amount of moving grains is universal (independent of the initial value $`\alpha `$), and given by $`R_s=1/(\gamma (1\beta ))`$, while $`H_s=\theta _<\nu _s`$. Since typically $`v\gamma d`$ with $`d`$ the grain diameter, we have in our rescaled units $`\gamma 1`$ showing that due to $`R_s\stackrel{<}{}1`$ non linear saturation terms can be neglected at the shock if $`\beta \stackrel{<}{}1`$. The lower branch of the envelop saturates for large $`t`$ exponentially fast with a characteristic time $`1/(\gamma \theta _>)`$ at $`x_{\mathrm{}}=[1+\mathrm{ln}(\alpha \gamma /\beta )]/(\gamma \theta _<)`$, which is always larger than $`x_s`$. This means that the shock stops propagating upwards. A large time expansion in the shock free range $`t<x<x_{\mathrm{}}`$ gives, taking the two leading terms of $`W`$, $`\nu (t,x)=(\alpha /\theta _<)\mathrm{exp}[\gamma \theta _<t(\theta _</\alpha )xe^{\gamma \theta _<t}]`$. Thus the slope is non monotonous within this range: after increasing for small times it relaxes again to the initial value $`\theta _<`$ as $`H_x(t,x)=\theta _<\mathrm{exp}[(\theta _</\alpha )xe^{\gamma \theta _<t}]`$.

## 6 Discussion

Let us summarize the major results of this paper, which could be explored experimentally. Starting from an initial profile made up of two different slopes, we find that shocks can occur after a finite time, depending on the value of the two slopes and the initial density of rolling grains. When shocks are absent, we find that the evolution surface profile is characterized by different velocities: the kink moves upwards with a velocity of the order of $`\alpha \gamma `$ for early times, while the edge of the “active” region moves downwards at a velocity which only depends on the initial slopes, and is smaller than the velocity of the grains. The final slope is shown to be the angle of repose; however, for finite size systems, one expects the final slope to be smaller by an amount which varies as $`1/\mathrm{ln}L`$. When a shock appears, we predict the time and position of this shock, as well as the density of rolling grains there, which takes a universal value. The shock is found to stop progressing upwards.

Our results are in disagreement with those of mp. For the situation considered here, they predict that the initial profile is rigidly shifted along straight characteristics. Therefore, for example, the final slope would be given by $`H_x(t,x=0)=\theta _<`$, which is completely different from our prediction of a decaying slope. The reason for this discrepancy comes from their implicit assumption that $`R_0(x)+H_0(x)+\mathrm{ln}(R_0(x))/\gamma =\mathrm{const}.`$, which does not hold in the cases considered here.

The method presented here can be extended to more general situations. For example, each profile $`H_0(x)`$ can be approximated by a piecewise linear function. Therefore, our analysis can be used to obtain analytical results for more complicated situations as, e.g., bumps or sinusoidal shapes. Another interesting situation is the case where $`R_0(x)`$ is localized in space. Applications of this method to the problem of ripple formation are under way.

Two important physical phenomena have been neglected: diffusion terms, which are expected to be important in the presence of shocks or in the case of a localized initial $`R_0(x)`$ (see ), and non-linear effects, which lead to a saturation of the static/rolling grains conversion term. A simple way to account for the latter effect is to replace the characteristics by straight lines of velocity $`\gamma R_{\mathrm{}}`$ as soon as $`R=R_{\mathrm{}}`$. The influence of a dependence of the velocity of grains on their density would also be worth investigating .

\***

This research was partly supported by the Deutsche Forschungsgemeinschaft (DFG) under grant EM70/1-1.