## Abstract

A non-local theory is proposed to model dense granular flows. The idea is to describe the rearrangements occurring when a granular material is sheared as a self-activated process. A rearrangement at one position is triggered by the stress fluctuations induced by rearrangements elsewhere in the material. Within this framework, the constitutive law, which gives the relation between the shear rate and the stress distribution, is written as an integral over the entire flow. Taking into account the finite time of local rearrangements, the model is applicable from the quasi-static regime up to the inertial regime. We have checked the prediction of the model in two different configurations, namely granular flows down inclined planes and plane shear under gravity, and we show that many of the experimental observations are predicted within the self-activated model.

## 1. Introduction

Although granular flows can be found in many industrial applications or geophysical situations, their description still represents a challenge. At first glance, the problem does not look very complex. Non-Brownian particles move under external forces, e.g. gravity, and interact through contact forces described by friction and collision rules. However, the dynamics of such an apparent simple system is extremely rich and still resists a unified description (Aranson & Tsimring 2006). In the absence of a unified framework, granular flows are commonly divided into three different regimes. At low velocity, the flow belongs to a quasi-static or plastic regime, and the response in terms of velocity profiles or concentration profiles is independent of the shear rate (Roux & Combes 2002). When the deformation rate increases, a viscous-like regime is observed and the material flows more like a liquid with a shear-rate-dependent stress (Forterre & Pouliquen 2008). In this intermediate regime, the particles experience multi-contact interactions. At very high velocity, a transition occurs towards a gaseous regime, in which the particles interact through binary collisions (Goldhirsch 2003).

From a theoretical point of view, the collisional regime is certainly the one where models are the most advanced. A kinetic theory of granular gases has been developed, which provides hydrodynamic equations suitable for describing dilute and highly agitated granular flows (Jenkins & Savage 1983; Goldhirsch 2003). The dense liquid regime, when the binary collision assumption breaks down (Lois *et al*. 2006), has attracted a lot of interest in the last 10 years, and progress has been made in our understanding of major properties (GDR Midi 2004; Forterre & Pouliquen 2008). A first important result comes from simple dimensional arguments. Unlike Brownian systems, the lack of an internal time scale for rigid grains imposes specific scalings for the constitutive laws (Da Cruz *et al*. 2005; Lois *et al*. 2006). Based on this observation, constitutive laws can be proposed under the assumption that the rheology is local, namely that the shear stress depends only on the local shear rate. The local rheology is then expressed in terms of a friction coefficient and a volume fraction depending on a single dimensionless number called the inertial number (GDR Midi 2004). The approach has been generalized and a tensorial version has been proposed that takes the form of a frictional visco-plastic rheology (Jop *et al*. 2006). It provides a relevant framework to describe dense granular flows in different configurations. Quantitative predictions have been obtained for flows down inclined planes (Da Cruz *et al*. 2005), for flows down piles (Jop *et al*. 2005), for the formation of roll waves on slopes (Forterre 2006) and for describing the collapse of granular columns (Lacaze & Kerswell 2009). However, limits exist within the visco-plastic description, which mainly concern transition to the quasi-static regime, when the flow becomes slower. First, the localization of deformation in narrow shear bands, observed, for example, in Couette flows or in silos, is not well predicted (GDR Midi 2004). Second, the exponential tail observed in velocity profiles for flow on a pile or in rotating drums is not predicted. Last, the stopping and starting properties of avalanches are not correctly described. An example is the minimum angle necessary for a flow to occur on a rough inclined plane. Experimentally, this angle is found to depend on the thickness of the granular layer, whereas the local rheology predicts a unique critical angle (Pouliquen 1999; Borzsonyi *et al*. 2008).

The visco-plastic constitutive law is thus not sufficient to capture the flow properties when the material approaches the quasi-static regime. The major problem comes from the assumption of a local rheology. There are several evidences in the literature that non-local effects are present. A first sign is the existence of large spatial correlations in dense granular flows. When approaching the quasi-static regime, correlations in the force network (Radjai & Roux 2003; Lois *et al*. 2006) and in the velocity fluctuations (Bonamy *et al*. 2002; Pouliquen 2004; Baran *et al*. 2006; Staron 2008) are observed, which become larger when the flow is slower. Another indirect but spectacular evidence of non-local effects is given by the experiment carried out by M. van Hecke (2008, personal communication). A ball resting on a sand pile sinks as soon as a shear band is created in the media, although it is created far from the position of the ball. This observation shows that the presence of a sheared region somewhere in the granular medium modified the rheological properties of the sample, even far from it. Several theoretical studies exist, which try to capture non-local effects and to describe quasi-static granular flows. Taking into account the rotation of the grains, using Cosserat-like plasticity models, is one way of introducing a non-locality (Mohan *et al*. 2002). Another approach has been developed by Kamrin & Bazant (2007), which consists of introducing fluctuations in plasticity models. Jenkins (2006) has proposed to modify kinetic theory by introducing a correlation length to model the occurrence of clusters. Other approaches consist of explicitly writing non-local equations, by assuming the existence of arches in the material, which transmit stresses in a non-local way (Mills *et al*. 1999).

The study presented in this paper is an attempt to propose an alternative non-local description of dense granular flows that is able to describe both the quasi-static and liquid regime. The idea is based on the picture of a self-activated process (Pouliquen *et al*. 2001). Shear motion somewhere in the material induces stress fluctuations, which can trigger shear motion somewhere else. We have used this idea in Pouliquen *et al*. (2001) to describe quasi-static granular flows. In this study, we improve the description by incorporating the fact that local rearrangements occur in a finite time. We then propose a non-local constitutive equation that is written as an integral equation relating the shear stress and the shear rate. This rheological description is then applied to two configurations, namely the flow down an inclined plane and the plane shear configuration under gravity. We show that important experimental observations are predicted within this approach.

The paper is organized as follows. In §2, we recall the basis of local visco-plastic rheology and discuss some of the limits in more detail. The self-activated model is presented in §3, and its predictions for flows down an inclined plane and plane shear under gravity are discussed. A discussion and conclusion are given in §4.

## 2. Local visco-plastic rheology

### (a) Plane shear

We first consider a granular material made of beads of density *ρ*_{p}, of diameter *d*, confined under a normal stress *P* and sheared at a constant shear rate (figure 1*a*). The beads are assumed to be rigid. Under such conditions, it has been shown that, for a large enough sample, the friction coefficient *μ*=*τ*/*P* and the volume fraction of the sample *Φ* depend on a single dimensionless number *I*, called the inertial number,
2.1The shapes of the functions *μ*(*I*) and *Φ*(*I*) are given in figure 1*a*,*b*. The friction coefficient *μ*(*I*) is an increasing function of the inertial number starting at a value *μ*_{s} at *I*=0 and saturating at a value *μ*_{2} at high *I*. The volume fraction decreases linearly with *I*. The inertial number can be interpreted as the ratio between a microscopic inertial time scale controlled by the pressure , which gives the typical duration of a rearrangement, and a macroscopic time scale . The inertial number is then the relevant parameter controlling the transition between the different flow regimes. Low values of *I* correspond to a quasi-static regime, namely the macroscopic deformation occurs on a time scale longer than the typical time of local rearrangements, whereas high values of *I* give rise to the transition to the kinetic regime.

If we further assume that the rheology is local, namely that the shear stress depends on the local shear rate, relations (2.1), obtained for the plane shear configuration, can be used as constitutive equations to describe other configurations. This approach gives rather good predictions for the velocity profile of flows down inclined planes and for flows on a pile (Jop *et al*. 2005). A tensorial generalization of equation (2.1) also gives quantitative predictions for three-dimensional flows, e.g. for flows between rough walls (Jop *et al*. 2006), for describing the roll waves instability observed on inclined planes (Forterre 2006) and for describing the collapse of a granular column (Lacaze & Kerswell 2009).

However, limits exist, which we discuss in detail in the next section for the case of flows down inclined planes and for the case of the formation of shear bands.

### (b) Flows down inclined planes

We consider a granular layer of thickness *h* flowing down a rough plane inclined at an angle *θ* from the horizontal (figure 2*a*). In a steady uniform regime, the equilibrium imposes the following stress distribution across the layer for the normal stress perpendicular to the plane *P* and the shear stress *τ*,
2.2The predictions of the local visco-plastic approach for this configuration are the following. First, from equation (2.2), the ratio between shear stress and normal stress is constant at . From the friction law in equation (2.1), we conclude that no flow is possible for inclinations such that . Second, for inclinations larger than *μ*_{s}, the fact that *τ*/*P* is constant implies that the inertial number *I* is constant across the layer. This implies that the volume fraction is also constant across the layer. Third, with *I* being constant, the shear rate has to vary as the square root of the pressure, by definition of the inertial number. The velocity profile is then given by a so-called Bagnold profile,
2.3The three main predictions of the local rheology are then the existence of a critical angle, a constant volume fraction across the layer and a Bagnold profile.

When the predictions of the model are compared with the experiments or results from molecular dynamics simulation, a good agreement is found for flows far from the flow threshold. Bagnold profiles are observed, and the volume fraction is constant (Silbert *et al*. 2001). However, when approaching slow flows, problems arise. First the threshold itself is not simply given by a critical angle, as predicted by the theory. The critical angle is a function of the thickness of the layer. In other terms, for a given inclination, there exists a critical thickness *h*_{stop}, below which no flow is possible (Pouliquen 1999). This is the first sign of a non-locality, as it means that the whole layer feels the presence of a rough rigid boundary. A second discrepancy is that the velocity profile close to the flow threshold no longer follows a Bagnold profile, but changes its shape, eventually becoming concave (Silbert *et al*. 2003). This is the second sign of a non-locality, as, for those profiles, the shear rate at the free surface does not vanish, whereas the shear stress does vanish (Rajchenbach 2003). Once again, this is incompatible with the simple local relation between the shear stress and shear rate postulated in Jop *et al*. (2006).

Therefore, although it gives correct predictions far from the flow threshold, the local rheology seems to be inappropriate to describe slow flow down inclined planes at low inclinations or for thin layers.

### (c) Shear bands in inhomogeneous configurations

Different flow configurations have been studied that lead to the formation of shear bands. We briefly discuss two very similar configurations. The first one is cylindrical Couette geometry, and the second one is plane shear under gravity (figure 2*b*).

A Couette cell consists of two concentric cylinders. The inner cylinder rotates and imposes a shear in the material. Under the assumption that no normal stress difference exists and neglecting the centrifugal force, the stress distribution in this cylindrical geometry is given by a constant normal stress *P* and a shear stress *τ*, which decreases with the distance to the centre of the inner cylinder,
2.4and
2.5where *R*_{i} is the radius of the inner cylinder, *r* is the radial position and *τ*_{0} is the shear stress applied on the inner cylinder.

The second geometry is a plane shear as in figure 1, but in the presence of gravity. The stress distribution is then given by
2.6with *ρ*=*ρ*_{p}*Φ*.

In both cases, the ratio *τ*/*P* between the shear stress and the normal stress decreases when going away from the moving wall. From the stress distribution, we can analyse the predictions of the local visco-plastic approach. First, in both cases, no flow is possible if the applied shear stress is not high enough. The material is sheared only if *τ*_{0}/*P*_{0}>*μ*_{s}. When this condition is fulfilled, the local approach predicts the development of a shear band, which spreads between the moving wall and a critical distance *z*_{c} (respectively, *r*_{c} for the Couette geometry) given by the condition *τ*(*z*_{c})/*P*(*z*_{c})=*μ*_{s}. Below *z*_{c} (respectively, above *r*_{c}), the material is static. The important point is that, within the local rheology framework, the width of the predicted shear band vanishes when approaching the flow threshold. This means that in the limit of quasi-static flows, the local rheology predicts an infinitely thin shear band, which reduces to a slip condition at the boundary with the moving wall. This is in contradiction with experimental observations (GDR MiDi 2004 and references therein). In both geometries, a shear band of finite thickness is observed, even in the quasi-static limit. The width of the shear band is independent of the shear velocity for slow flows, equal to 5–10 particle diameters. The size of the shear bands increases when increasing the velocity of the moving boundary when entering the inertial regime.

In conclusion, the local rheology does not capture the existence of shear bands in the quasi-static limit. In the next section, we present a simple model that tries to capture both the inertial and the quasi-static regime.

## 3. Towards a non-local model

In this section, we present a non-local model that is able to capture both the quasi-static and the liquid regime of granular flows. The theory is an extension of the approach introduced in Pouliquen *et al*. (2001), which is based on the idea of a self-activated process. We show how to develop a self-activated model in which the finite time of the local rearrangements is taken into account, which gives the possibility of describing the inertial regime. The constitutive law obtained in this approach is applied to the two geometries discussed previously, namely the flow down inclined planes and the plane shear under gravity.

### (a) A self-activated process

The non-local model is based on the idea of a process activated by fluctuations. Thermally activated processes (Eyring 1936) have been successfully applied to explain many physical phenomena, including plasticity (Feda 1989), friction (Heslot *et al*. 1994) or the rheology of complex fluids (Fielding *et al*. 2000). The basic idea is that thermal fluctuations help microscopic elements of the material to yield and overcome the energy barrier in which they are trapped. This picture is attractive for granular materials as, in a dense regime, grains are trapped in cages formed by their neighbours. However, for granular material, thermal fluctuations are negligible. Several studies (Pouliquen & Gutfraind 1996; Debregeas *et al*. 2001; Behringer *et al*. 2008) have proposed that stress fluctuations could play the role of temperature.

The idea of the model is presented in figure 3. We consider a unidirectional, stationary shear flow of particles with a mean velocity *u*(*z*) along *x*. We call *P*(*z*) and *τ*(*z*), respectively, the pressure and the shear stress at level *z* and the shear rate. The flow takes place between *z*=*a* and *b*. The particle diameter is *d* and the particle density *ρ*_{p}. The main idea is the following: a shear motion occurring at level *z*′ between two layers induces shear stress fluctuations *δσ*_{z′→z} at level *z* , which can help the particle in *z* to jump to the next hole. The motion in *z* is then activated by the fluctuations created by *z*′. This picture is non-local in the sense that one has to take into account the fluctuations coming from all levels *z*′ to know the motion in *z*. From this idea, one can write the shear rate in *z* as the sum of jump probabilities induced by all levels *z*′,
3.1where *f*_{z′→z} is the frequency of fluctuations sent by *z*′ on *z*, and and are the probabilities for *z* to jump to the right or to the left because of the perturbation from *z*′. The frequency and probabilities have to be expressed in terms of the stresses and the shear rate.

The frequency *f*_{z′→z} at which fluctuations are sent by *z*′ is simply given by the shear rate ,
3.2

To find the expression for the probability of jump , we follow the argument given in Pouliquen *et al*. (2001), based on stress fluctuations. However, we introduce a new ingredient to take into account the finite time it takes for a particle to jump. The argument is the following.

When particle *z* receives the stress fluctuation sent by particle *z*′, two conditions are necessary for a jump to be triggered. First, particle *z* has to be at rest, which means it has to have finished its previous jump, otherwise the received fluctuation is assumed to have no effect. Second, the stress fluctuation has to be strong enough for the local yield condition to be reached. We assume that these two conditions are independent. As a result, the jump probability is written as the product of two probabilities and , corresponding to the two conditions: .

The probability that particle *z* is at rest and ready to receive a fluctuation is simply related to the ratio between the time spent during the jumps and the total time. As discussed in §2*a*, the typical time of rearrangement is given by . This implies that has to be a function of the inertial number . The probability has to be equal to one at low *I* when the jumps are very rapid compared to the deformation time, and tends to zero when *I* tends to infinity, i.e. when the particle spends most of the time experiencing rearrangements. We then choose the following expression for the probability :
3.3where *I*_{0} is a constant.

The probability concerning the level of the stress fluctuation is chosen by introducing a local yield criterion given by a friction law with a coefficient of friction *μ*_{2}. The probability (respectively, ) that the yield condition is reached is then equal to the probability that the shear stress fluctuation *δσ*_{z′→z} is higher than *μ*_{2}*P*(*z*)−*τ*(*z*) to induce a jump to the right (respectively, *μ*_{2}*P*(*z*)+*τ*(*z*) to induce a jump to the left).

Under these assumptions, equation (3.1) can be written as 3.4

Note that the above description is not valid when |*τ*|/*P*>*μ*_{2}. In this case, the stress barrier vanishes and the description no longer holds.

The last step consists of writing the distribution of stress fluctuations in order to be able to describe the probability of jumps. To do so, we first assume that the stress fluctuation *σ*_{z′→z} is a decreasing function of the distance (*z*−*z*′). We choose the following function:
3.5where *δσ*_{0}(*z*′) is the maximum value of the fluctuations and *β* is a dimensionless constant of order unity, which measures the spatial extent of the stress fluctuation. The choice of the decreasing function is arbitrary, but we have checked that the results do not qualitatively depend on the exact shape of the function.

The amplitude *δσ*_{0}(*z*′) is a random variable with a mean, which is of the order of the local pressure *P*(*z*′). Several studies of the distribution of forces in granular media suggest that the distribution is exponential over a large range of forces (Mueth *et al*. 1998; Radjai *et al*. 1999; Snoeijer *et al*. 2004; Majmudar & Behringer 2005). We then assume that *δσ*_{0}(*z*′) follows an exponential distribution, i.e. the probability *pd*(*δσ*_{0}) to have the value *δσ*_{0} within a range *d*(*δσ*_{0}) is equal to
3.6

We can then easily show that the probability of jumps, equation (3.1), can also be expressed in terms of an exponential 3.7

It is interesting to note that the above expression is similar to an activated process, where the role of the energy barrier is played by the distance to the local yield criterion, *μ*_{2}*P*(*z*)−*τ*(*z*), and the role of the temperature is played by the mean stress fluctuation. The possibility of expressing the jump probability using a Boltzman-like exponential term comes from the exponential distribution of the stress fluctuations.

Finally, transforming the discrete sum in equation (3.4) by an integral, we can write the final expression for the shear rate as a function of the shear stress 3.8

This integral equation can be seen as a non-local rheological law relating the stresses *P*(*z*) and *τ*(*z*) and the shear rate . Knowing the stress distribution in different configurations, we are able to predict the velocity profile. Three parameters have been introduced in the model: the microscopic friction coefficient *μ*_{2}; the characteristic inertial number *I*_{0}; and the spatial extension of the stress fluctuation *β*. In the following, we analyse the prediction of this simple self-activated model for plane shear, for flows down inclined planes and for plane shear under gravity.

### (b) Plane shear

We consider the simple rheological test of an infinite medium without gravity, sheared at a constant shear rate and confined under a pressure *P*. From equation (3.8) with the limits of the integral going to infinity, one obtains the following relation between *τ*, *P* and :
3.9

The macroscopic coefficient of friction *τ*/*P* is a function of the inertial number *I* and varies between two values when increasing the shear rate, as plotted in figure 4. The minimum value *μ*_{s}, which gives the flow threshold, depends on the microscopic friction coefficient *μ*_{2} and on the coefficient *β*. The upper limit, which is reached for very rapid flows and high values of *I*, is simply given by the microscopic friction coefficient *μ*_{2}. It is interesting to note that, within this model, the macroscopic coefficient of friction *μ*_{s} is lower than the local friction coefficient and results from a non-local self-activated process. The shape of the friction law predicted by the model is thus very similar to the one observed in the experiments and simulations. From a quantitative comparison with data from three-dimensional flows down inclined planes, we can estimate the three parameters we have introduced in the model. The best fit gives (figure 4) *μ*_{2}=0.62, *I*_{0}=0.26 and *β*=3.3. For these parameters, one obtains *μ*_{s}=0.34.

Once the model has been calibrated on the simple plane shear geometry, we can now study the prediction for other flow configurations.

### (c) Flows down inclined planes

We now address the problem of the flow down a rough inclined plane. The two control parameters are, in this case, the inclination *θ* of the plane and the thickness of the flow *h* (inset of figure 5). In this configuration, the material is bounded by a free surface *z*=*h* and a rough bottom at *z*=0, which in the model is simply the first row of particles. We have seen that the stress distribution is given by and . In order to solve equation (3.8), we discretize the integral and use an iterative method to find the solution. As a guess value, we use the solution predicted by the local rheology equations (2.1). This procedure can be applied for different thicknesses *h* and different inclinations *θ*. The solution gives , and the velocity profile *u*(*z*) is obtained by integrating with the boundary condition *u*=0 at *z*=0.

The first result of the model is that a solution exists only in a limited region of the parameter space (*h*, *θ*), as shown in figure 5. For a given inclination *θ*, there exists a critical thickness *h*_{stop}, below which no flow is possible. In order to compare with the experiments (figure 5*a*), where the thickness *h*_{stop} is defined as the thickness without counting the glued layer, we have plotted in figure 5*b* the critical thickness predicted by the model minus one layer. We observed that the non-local model correctly reproduces the observed critical thickness. The agreement is not perfect, but the prediction gives the right order of magnitude. A thin layer, has more difficulty to flow than a thick layer. In the model, this effect simply comes from the truncature of the integral in equation (3.8). In other words, in a thin layer, the amount of fluctuations received by the grains is less, which implies that the system has to be closer to the local yield criterion, i.e. at a higher angle, in order to flow.

The second result concerns the flow regime. In the region where solutions exist, we have computed the depth-averaged velocity . In figure 6*a*,*c*, the mean velocity is plotted as a function of the thickness for different inclinations for both the experiments of Pouliquen (1999) and the model. The agreement is quantitatively correct. Experimentally, it has been observed that plotting the Froude number as a function of the thickness *h* rescaled by *h*_{stop} allows one to collapse all the data, as shown in figure 6*b*. This flow rule has been observed in other numerical or experimental studies (Forterre & Pouliquen 2003; Silbert *et al*. 2003; Borzsonyi *et al.* 2008). In order to test the model, we have plotted the same scaling using the *h*_{stop} previously determined. Although the scaling has a tendency to gather the data, the collapse is not perfect and is less striking in the model than in the experiments. We will discuss this point in the last section of the paper.

Finally, in figure 7, we compare the shape of the velocity profile predicted by the model with results from three-dimensional molecular dynamics simulations (Silbert *et al*. 2003). The different profiles are obtained for a fixed angle (24°) and for different thicknesses *h*, varying from 38*d* down to 4*d*. One clearly observes that the profile follows the Bagnold profile for thick layers, but exhibits a different shape when decreasing the thickness. For the thinnest flow, a concave profile is observed. The non-local theory qualitatively reproduces this tendency. The profile changes shape when decreasing the thickness (figure 7*b*). However, for the thinnest flow, the predicted profile is less concave than the one observed in the simulation.

From this study on the flow down inclined planes, we can conclude that many of the features observed experimentally are well captured by the model. However, some details are still not properly described.

### (d) Plane shear under gravity

In this section, we analyse the prediction of the model in the case of plane shear in the presence of gravity. A pressure *P*_{0} is applied on a rough top plate in contact with a semi-infinite granular material. A shear is then imposed by either a velocity *u*_{s} of the plate or a shear stress *τ*_{0}. Under gravity, the stress distribution within the medium is given by equations (2.6). In this configuration, the flow is controlled by two parameters: the level of confining pressure compared to the gravity *P*_{0}/*ρgd* and the ratio between shear stress to normal stress *τ*_{0}/*P*_{0}. Knowing the stress distribution, we can predict the velocity profile from equation (3.8). The results are summarized in figure 8.

The first result is that a minimum shear stress is necessary to create a flow. Interestingly, the theory predicts that the minimum friction coefficient *τ*_{0}/*P*_{0} for the top plate to move depends on the confining pressure *P*_{0}/*ρgd*. At low confining pressures compared with gravity, it is more difficult to shear the material. We are not aware of any experiments in which the confining pressure has been significantly varied to observe such an effect. When applying a shear stress higher than the threshold value, the model predicts that the top plate starts to move at a velocity *u*_{s} and a shear band develops in the material, as shown by the velocity profile plotted in figure 8*a*. The shear band is localized at the top plate, the velocity being maximum at the top and decreasing rapidly to zero. To characterize how the velocity profile varies when varying the parameters, we study how the plate velocity *u*_{s} and the typical thickness *λ* of the shear band vary when varying the applied shear stress. The relation between *u*_{s} and the friction coefficient *τ*_{0}/*P*_{0} is plotted in figure 8*c* for different levels of confining pressure. It is observed that the velocity increases from zero to a finite value when *τ*_{0}/*P*_{0} increases from its threshold value to the maximum value *μ*_{2}. The existence of a critical velocity when *τ*_{0}/*P*_{0} approaches *μ*_{2} indicates that no solution exists in the model with a plate moving faster than . The maximum speed corresponds to the development of a velocity profile with an infinite shear rate at the top. This prediction is not linked to the non-local rheology and is also obtained within the local visco-plastic rheology, as discussed by Jop *et al*. (2005) for flows on a pile. Jop *et al*. (2005) have observed that going faster than this critical velocity seems to correspond to transition to a dilute kinetic regime.

However, the non-locality introduced in the model plays a crucial role in the shape of the velocity profile and in the variation of the width of the shear band. First, the velocity profile exhibits a creeping tail. Unlike the prediction of the local rheology, the velocity is never strictly equal to zero, as shown by the inset of figure 8*a* presenting the velocity profile on a semi-log scale. An exponential tail is predicted, as observed in many experiments. Second, the non-locality plays a role in the width *λ* of the shear band. We define *λ* as the position where the velocity is equal to 10 per cent of the plate velocity *u*_{s}. The variation of *λ* with the plate velocity is plotted for different confining pressures in figure 8*d*. The model predicts that the shear band has a finite thickness in the quasi-static regime when *u*_{s} tends to zero, and increases when increasing the plate velocity. The non-local rheology then predicts a transition between a quasi-static regime where the shape of the velocity profile is independent of the velocity, to an inertial regime where the shear rate dependence induces the thickening of the shear band. We have not been able to compare these predictions with experimental data, as, to our knowledge, all experiments in this configuration take place in the limit of very slow flows (Losert *et al*. 2000; Siavoshi *et al*. 2006; Divoux & Géminard 2008).

## 4. Discussion and conclusion

In this paper, we have presented a simple approach to describe the flow of granular media as a self-activated process. Based on the picture that a shear motion somewhere induces stress fluctuations, which can trigger a rearrangement somewhere else, we have been able to write a non-local constitutive equation. In this formulation, the shear rate is written as an integral over space of the shear rate times a probability, which depends on the stress distribution. Our purpose was to provide a framework capable of unifying the quasi-static flows and the inertial flows and to go beyond the simple visco-plastic rheology, which has been successfully applied in the inertial regime.

This goal has been partially achieved, and many experimental observations in different flow configurations are predicted by the model. First, the quasi-static regime and the shear bands observed in plane shear under gravity, in silos or in the Couette cell are well predicted by the theory. Second, the model correctly describes the shear-rate dependence of the friction coefficient for plane shear in the inertial regime. The friction coefficient is found to be an increasing function of the inertial number and has the same shape as the one proposed by the visco-plastic local approach. Third, the transition between the quasi-static and the inertial regime can also be described. We have tested this intermediate flow regime for flows down inclined planes. In this configuration, the self-activated model significantly improves the predictions of the local rheology. The flow threshold is correctly reproduced, and the existence of a minimum thickness below which no flow is possible is captured by the model. The existence of velocity profiles that differ from the Bagnold profile when flow occurs close to the flow threshold is also predicted by the model. Finally, for confined flows, such as plane shear under gravity or Couette flow, the model correctly predicts that the shear band becomes thicker and thicker when increasing the velocity.

However, discrepancies exist, revealing weaknesses of the approach. A first problem concerns the result for flows down inclined planes. The model does not reproduce the collapse of the flow rule when it is expressed in terms of the Froude number as a function of the thickness over the stopping height. Whether this collapse observed experimentally is a coincidence or a signature of a more profound physical meaning is still a matter of debate. In our model, the quality of the collapse strongly depends on the parameters. A better collapse than in figure 6*d* can be obtained, but for another set of parameters, which then does not correctly fit the velocities. This shows that, in the self-activated model, the collapse with the stopping height (when it exists) is more a coincidence than an intrinsic mathematical property.

A perhaps more serious limit of the present approach concerns the hysteresis observed in granular flows between the initiation and the stopping properties of the flow. It is well known that, in order to flow, an initially static layer has to be inclined at an inclination angle higher than the angle at which it stops. In our self-activated approach, a hysteresis exists. The zero shear-rate state is a solution of equation (3.8) as long as . This means that the model predicts that a static granular layer will start to flow only if the inclination reaches the critical inclination . However, experimentally, it is observed that the starting angle is not constant, but also depends on the thickness of the layer and follows a similar dependence to the stopping angle (Pouliquen & Renaut 1996; Daerr & Douady 1999). This indicates that the same non-local effect influences the static, as well as the dynamic, critical angle.

This discrepancy suggests that the non-locality introduced in the self-activated model is perhaps not sufficient. In our model, non-local effects take place through the integral of the stress fluctuations and are observed when a truncature of the integral takes place. However, the existence of finite-size effects, starting from a static pile, strongly suggests that the energy barrier itself and the local yield criterion are perhaps controlled by non-local effects. New ideas based on recent concepts of soft modes could be relevant to better understand this issue (Wyart 2009).

A last limit of the approach presented in this study concerns the volume fraction. The model is written without taking into account variations of the volume fraction. The volume fraction is considered as a slave variable, which can be computed afterwards, once the velocity profile is known, using the local law *Φ*(*I*). However, dilatancy effects could be important and could perhaps participate in non-local effects.

In conclusion, the self-activated approach developed in this paper is a first step towards a non-local rheology for dense granular flows. The derivation of the model remains phenomenological and contains many assumptions. However, the success obtained in describing many experimental observations represents an encouragement to pursue this direction.

## Acknowledgements

This work has been supported by the Indo-French Centre for the Promotion of Advanced Research (IFCPAR) and by the Agence National de la Recherche (ANR).

## Footnotes

One contribution of 12 to a Discussion Meeting Issue ‘Colloids, grains and dense suspensions: under flow and under arrest’.

- © 2009 The Royal Society