Nonlocal modeling of granular flows down inclines

Ken Kamrin *a and David L. Henann b
aDepartment of Mechanical Engineering, MIT, Cambridge, MA, USA. E-mail:
bSchool of Engineering, Brown University, Providence, RI, USA. E-mail:

Received 19th August 2014 , Accepted 29th September 2014

First published on 6th October 2014

Flows of granular media down a rough inclined plane demonstrate a number of nonlocal phenomena. We apply the recently proposed nonlocal granular fluidity model to this geometry and find that the model captures many of these effects. Utilizing the model's dynamical form, we obtain a formula for the critical stopping height of a layer of grains on an inclined surface. Using an existing parameter calibration for glass beads, the theoretical result compares quantitatively to existing experimental data for glass beads. This provides a stringent test of the model, whose previous validations focused on driven steady-flow problems. For layers thicker than the stopping height, the theoretical flow profiles display a thickness-dependent shape whose features are in agreement with previous discrete particle simulations. We also address the issue of the Froude number of the flows, which has been shown experimentally to collapse as a function of the ratio of layer thickness to stopping height. While the collapse is not obvious, two explanations emerge leading to a revisiting of the history of inertial rheology, which the nonlocal model references for its homogeneous flow response.

1 Introduction

Nonlocal effects in granular media manifest in myriad different ways. At the origin of the nonlocality is the finite size of the grains themselves, inducing cooperative behaviors that defy local rheological description. Examples include grain-size-dependent shear features in the steady flow profiles of granular media. A local rheology can be extracted from uniform simple shearing data of a granular media;1 however, nonuniform steady flows of the same material can be seen to violate such a relation,2–4 as the grain-size sets up an internal length-scale that effectively penalizes variations in flow-rate over space. A more recently observed nonlocal manifestation is the mechanically-induced creep effect, also known as “secondary rheology”.5–7 A highlight of this phenomenon is that flow anywhere in a granular media removes the yield condition everywhere, and the rheology of probes in the body varies with distance from the source of motion. Some of the most compelling demonstrations of nonlocality in granular media can be observed in the behavior of grains on an inclined surface.8–11 Contrary to local rheological models, which predict a thickness-independent repose angle, experiments and discrete simulations verify that granular layers have a critical “stopping height” proportional to the grain size — layers thinner than this value come to a stop, whereas thicker layers admit steady flow down the incline. Examples of past continuum approaches aimed at modeling this effect include the theory of partial-fluidization,12 and extensions of the kinetic theory.13,14

Recently, a nonlocal rheological model based on the concept of “granular fluidity” has shown itself able to reconcile the issues of grain-size dependent shear features as well as the mechanically-induced creep effect in granular media.15–19 While these demonstrations pertain primarily to the way in which grain size influences spatial fields in materials that are driven to flow, nonlocality as studied in inclined plane flows is of a relatively different nature, particularly the notion of a stopping height, i.e., whether a material flows at all. Despite this distinction, in this paper we shall show that the nonlocal fluidity model also captures the phenomenology observed in the inclined plane geometry.

Herein, we compare theoretical results to known results (i.e., experimental or DEM) on the issues of (i) the presence of a stopping height and how it varies with inclination, (ii) the variation in shape of the flow profile as inclination increases above the stopping height, and (iii) the dependence of mean flow speed on the ratio of the thickness to the stopping height. These three phenomena have been well-studied (see aggregated data in MiDi11) and constitute the major manifestations of nonlocality in inclined plane flows.

2 From local to nonlocal granular rheology

We begin by describing the inertial granular rheology and the nonlocal granular fluidity model (per the review in Kamrin and Koval17). The inertial rheology for steady flow of dense granular media relates the local stress state and the strain-rate1,20,21 and is expressed via the dimensionless relationship:
image file: c4sm01838a-t2.tif(1)
In the above, μ is the ratio of shear stress τ and normal pressure P, and I is the inertial number, where [small gamma, Greek, dot above] is the shear rate, d is the mean grain diameter, and ρs is the grain density. The function for μloc(I) is empirically fit and is typically characterized by a yield coefficient μs such that μloc(I → 0) = μs. Under quasi-static circumstances, for which μ does not increase substantially above μs, a linear fit is often used,
μloc (I) = μs + bI,(2)
where b is a dimensionless material parameter. As μ increases further above μs, curvature of the function is observed with an asymptote at an upper-limiting value of μ, denoted μ2. A common fit is20
μloc (I) = μs + Δμ/(I0/I + 1)(3)
where Δμμ2μs and I0 ≡ Δμ/b.

The inertial rheology works well in describing uniform flows (e.g., planar shearing) over a wide range of flow rates.1 However, in non-uniform flow geometries in zones of slowly-flowing, quasi-static material the one-to-one inertial relation between μ and I is violated.3 The behavior displayed in these zones is definitively nonlocal; the bulk stress/strain-rate behavior at steady-state changes as the macroscopic geometry varies, even when the local kinematics are identical.

The “nonlocal granular fluidity” (NGF) model may provide the solution to the above issues. It has demonstrated the ability to quantitatively predict granular flows in many disparate geometries, with predictions verified in 2D and 3D as compared to both discrete-particle simulations15 and experiments.16 It is the first continuum model to quantitatively describe all experimentally-obtained flow data in the complex split-bottom family of geometries.16 It has also been shown to correctly capture other nonlocal phenomena such as nonlocally induced material weakening, whereby the motion of a boundary removes the yield strength of the material everywhere, permitting far-away loaded objects to creep through the grains when otherwise they would remain static.19

Our existing work has implemented the NGF equations in a reduced, steady-state-only form

image file: c4sm01838a-t3.tif(4)
where the linearized version of μloc (I), eqn (2), was used since flows of interest were all close to quasi-static. The field g is a state parameter called the granular fluidity, and ξ is the plastic cooperativity length, which is proportional to d. Note that in planar shear, flow gradients vanish and the above reduces appropriately to the local rheology, but in the presence of gradients, the Laplacian term “spreads” fluidity based on ξ. In our previous work,15,16 we have verified that ξ is, in fact, not a constant length but satisfies
image file: c4sm01838a-t4.tif(5)
roughly in-line with prediction of the kinetic-elasto-plastic (KEP) theory on which other fluidity models are based.22 The parameter A, the dimensionless nonlocal amplitude, is the only new parameter in the model, which quantifies the spatial extent of cooperativity in the flow. We note that other models based on the KEP mechanism are shown to yield size-dependent jamming in other geometries,23 which is encouraging for our current study.

2.1 Dynamical system for fluidity

In Henann and Kamrin,18 we describe how the steady-state-only NGF model – i.e., eqn (1), (2), (4) and (5) – is obtained in its entirety from the steady flow limit of a thermomechanically derivable dynamical system for g. One treats the fluidity and its gradient as kinematic state variables with separate free-energy contributions. By selecting the corresponding free-energies in a fashion that preserves the linear inertial rheology in uniform flow conditions, eqn (2), we obtain the following system:
image file: c4sm01838a-t5.tif(6)
where t0 > 0 is a constant time-scale. The steady-state-only model arises as the approximate solution of the stable, steady behavior of g in the above dynamical partial differential equation (PDE). The approximation is valid as long as g = gloc + δ for some small function δ, where gloc emerges as the stable solution when the Laplacian term above vanishes. To switch the local response to the more robust form, eqn (3), the same analysis can be reapplied giving
image file: c4sm01838a-t6.tif(7)
When the above is reduced to a steady-state only model, the nonlocal system obtained matches the aforementioned one in the appropriate limit of μ near μs, as it should, with ξ maintaining the same inverse square-root divergence behavior. The precise form of ξ in the steady-state system corresponding to eqn (7) is
image file: c4sm01838a-t1.tif(8)
Note that the limit of eqn (8) as μ2 goes to infinity corresponds to eqn (5).

The dynamics of g presented above have yet to account for bistability of granular flows nor do they account for the totality of transient effects that occur in a sheared granular media. To the former point, recent experimental evidence suggests non-monotonicity of the local response (contrary to eqn (2) and (3)) may exist giving rise to bistable flow hysteresis in certain circumstances.24 To the latter point, Reynolds dilation during shear initiation induces transient variations in material strength, which are commonly described using critical-state models25 but which are not yet included in our fluidity dynamics. Instead, our model above intends to provide is an accurate long-term dynamical behavior of the material when passing through a sequence of developed flowing states. We note that the process of obtaining a steady-state-only relation from a fully dynamical form has a history within fluidity-based modeling for other amorphous media.22 The dynamical form of the nonlocal fluidity model shows mathematical similarities to order-parameter-based rheological approaches, which also feature a diffusing state variable.12

3 Strengthening due to thinness

In the inclined plane geometry, a layer of thickness H is inclined to an angle θ. Assuming a packing fraction of ϕ and acceleration of gravity G, the stress distribution in the layer obeys
μ = tan[thin space (1/6-em)]θ, P = ϕρsGz[thin space (1/6-em)]cos[thin space (1/6-em)]θ.(9)
See Fig. 1(a) for reference. Accordingly, if a local relation were applied, either eqn (2) or (3), one would predict a universal angle of repose θr = arctan[thin space (1/6-em)]μs, i.e., any layer tilted above θr would be predicted to flow.

image file: c4sm01838a-f1.tif
Fig. 1 (a) Setup of the inclined plane flow geometry. (b) Theoretically predicted Hstop locus (−), eqn (11), as calibrated for glass beads on a fully rough base, compared to experimentally determined values (○) from Pouliquen.8

A signature of the cooperativity of granular motion is the fact that this universal repose angle is contradicted in experiments. As shown initially by Pouliquen8 and verified by others,9–11 the angle at which an initially flowing layer of grains comes to a stop, θstop, depends sensitively on the thickness of the layer when H is small. Inverting, one can extract a function Hstop (θ) for every granular media and substrate, which represents the critical thickness at which a flowing layer at a certain angle would arrest.

Unlike our past studies, which have focused on nonlocal effects in media that is necessarily flowing, the problem of size-dependent strengthening of thin layers often presents conditions that do not satisfy the g = gloc + δ approximation necessary to reduce the mathematics to steady-state-only form; i.e., material stops (g = 0) although the inclination can be steep enough for gloc to be notably greater than zero. To study inclined flows, we must revert to the primitive, dynamical form of the nonlocal fluidity model, eqn (7), to ensure accuracy of solutions, which should be valid up to and including steep inclination angles.

We compare the results of our model to the experiments of Pouliquen,8 where glass beads were used as the media. The local continuum parameters for glass beads have been calibrated in Jop et al.20 and are μs = 0.382 = tan(20.9°), μ2 = 0.643 = tan(32.76°), b = Δμ/I0 = 0.938, and ρs = 2450 kg m−3. The nonlocal amplitude for glass beads was calibrated in Henann and Kamrin16 to be A = 0.48. We can therefore apply our model in this geometry without additional parameter fitting, as long as we can identify accurate boundary conditions for the fluidity and basal slip velocity. We have tried two different fluidity boundary conditions in our past studies, chosen primarily based on simplicity, with the understanding that the choice is less important far from system boundaries due to the source-diffusion nature of the fluidity system. However, due to the thinness of the layers we wish to consider here, the accuracy of the boundary condition has a much larger influence on the flow behavior. To remove this issue, we opt instead to extract the fluidity boundary conditions directly from existing flow data. In discrete particle simulations of Silbert et al.,9 who also used spherical particles, it is apparent that adjacent to a fully rough boundary, the shear-rate (and velocity fluctuations) approach an approximately vanishing state, and the mean velocity vanishes, indicating no slip. At or near the free surface, there is usually (though not always) a zone where the shear-rate appears to level off and a vanishing shear-rate gradient is observed. Likewise, we will assume g = 0 and ν = 0 along the rigid surface at z = H and take ∂g/∂z = 0 at the free surface, z = 0. Since these boundary conditions are homogeneous, they always permit both flowing and static solutions.

The phase diagram of flowing and non-flowing states according to the nonlocal model can be obtained by determining the conditions necessary for the global g = 0 solution to be stable, as this decides if a system perturbed to flow returns to a static state. In eqn (7), a perturbation about g = 0 renders the g2 term negligible. The remaining PDE is linear and, substituting μ = tan[thin space (1/6-em)]θ, it is solved by

image file: c4sm01838a-t7.tif(10)
where C is a constant prefactor. The solution decays to g = 0 as long as the term in the exponent is negative. Hence, a perturbed layer solidifies if
image file: c4sm01838a-t8.tif(11)
A similar mathematical technique was utilized in Aranson and Tsimring.12 The above shows that for tall layers, HO(d), the material is predicted to have a common repose angle of θ1 = arctan(μs). However, for thinner layers, an inclination higher than this value is needed to stop flow. This applies up to an upper limit of θ2 = arctan(μ2) at which all layers, independent of height, are predicted to flow. The existence of two critical angles having these properties has been verified as a common trait of Hstop data in inclined plane flows in multiple experiments and particle simulations involving different materials and substrates.11 When it comes to specific functional forms for Hstop data, several fit curves have been proposed, a topic we shall return to later, but all invoke critical θ1 and θ2 values as described here. It bears noting that the derived Hstop relation, eqn (11), is linearly proportional to the cooperativity length, eqn (8), i.e., Hstop(arctan[thin space (1/6-em)]μ) = (π/2)ξ(μ). This is in line with experimental observations of Pouliquen,26 who reported that velocity correlation lengths in inclined plane flows vary with θ (and hence μ) in a similar manner as Hstop. We emphasize that this is a derivable consequence – not an underlying assumption – of the dynamical system, eqn (7). Further, one will obtain the same result when working with the simplified dynamical system, eqn (6), albeit with slightly modified functional forms.

For a direct comparison, in Fig. 1(b), the above predicted form for Hstop(θ), using the parameters for glass beads, is compared against Pouliquen's experimentally determined values using glass beads with a fully rough base. It is worth pointing out that the model's result was obtained using the same continuum parameters that were used to successfully predict steady flow fields of glass beads in split-bottom cells and other geometries.16 Yet here, the question is of a different nature, one of predicting input conditions for flow stoppage rather than velocity profiles in a flowing body.

We note that a size-dependent starting height Hstart(θ) – the height at which a static layer initiates flow – is also observed experimentally. Its curve depends on the preparation of the static layer and differs somewhat from Hstop, though it shares the same qualitative features (e.g., θ1 and θ2 values where the curve asymptotes and vanishes respectively). Our analysis above is tailored to Hstop since it reflects the limiting behavior of a flowing solutions as H (or θ) is decreased, however we note that without a bistable term in our dynamics, a distinct Hstart curve does not arise theoretically.

4 Velocity profiles

For angles between θ1 and θ2, and H exceeding Hstop, the material has a steady flow state down the incline. As reported collectively in GDR MiDi11 and explored directly in Silbert et al.,9 the shape of the flow profile varies from concave up, to roughly linear, to concave down as H is increased. The concave down shape observed for large H is well-fit by the so-called “Bagnold profile” obtainable by integrating the strain-rate predicted by the inertial rheology through the thickness of the layer. The result is of the form
νBag(z) ∝ H3/2z3/2.(12)
Fig. 2(a) shows the discrete particle simulation data of Silbert9 for inclination angle θ = 24° and various values of H > Hstop(θ = 24°). We have computed numerical solutions to our model at the same inclination angle and several layer thicknesses. We note that the stopping curves of the simulated material and glass beads are different; at θ = 24°, the simulated particles have an Hstop/d value more than twice that of the experiment. To make an appropriate albeit qualitative comparison, we select our H values so as to (roughly) match the relative heights H/Hstop used in the discrete simulations. The results are in Fig. 2(b). The model predictions are calculated by evolving eqn (7) using finite-differences in MATLAB, using Δxd and the Bagnold profile as the initial guess.

image file: c4sm01838a-f2.tif
Fig. 2 (a) Discrete particle simulation data of Silbert et al.9 for θ = 24° and many layer thicknesses (increasing from top to bottom). Dashed line is the Bagnold profile. (b) Theoretical profiles for the same relative heights, H/Hstop ≅ 1.1, 1.9, 2.7, 4.1, 6.4, as well as H/Hstop = 20.

The model shows that for H near Hstop the profiles display a concave-up feature in agreement with the particle simulations. It could be said that the concave up appearance is due to the comparative unimportance of the g2 term in eqn (7) when flows are slow, which, if neglected, gives a cosinusoidal solution for g and ν ∼ 1 − sin(zπ/2H). As H increases to be much greater than Hstop, the profiles approach the Bagnold profile, which agrees with the discrete simulations and is sensible since the relative importance of the particle-size-dependent term in eqn (7) diminishes in a tall flow; if neglected, the remaining terms give exactly the inertial rheology at steady flow. The transition between these extremes is marked by velocity profiles displaying a linear character. It can also be observed that, in agreement with the data and in contrast to the Bagnold profile, the model does not always require a zero strain-rate at the free surface. Because the value of A determines the degree of nonlocality, and hence the degree of deviation from the Bagnold profile, it is clear that as A increases, the range of heights over which the profile remains concave up will increase. However, since increasing A also increases Hstop proportionally, it is the value of H/Hstop that has the most direct influence over the shape of the flow profile.

5 Discussion of the Froude number

We now discuss the mean flow speed of granular media down an incline. Pouliquen,8 and others thereafter,9–11 have observed that the Froude number of the flow, image file: c4sm01838a-t9.tif, appears to be a relatively well-defined, one-to-one function of the relative height, H/Hstop, for all angles and heights as long as H is not close to Hstop. In this range, they find
image file: c4sm01838a-t10.tif(13)
where, for glass beads, β = 0.136. When results of the nonlocal theory are plotted in this fashion, Fig. 3(a), we do not observe a one-to-one correspondence between the Froude number and the relative height. We have two comments on this point, that may explain the discrepancy.

image file: c4sm01838a-f3.tif
Fig. 3 Predicted Froude number, image file: c4sm01838a-t19.tif, at inclinations θ = 22.5(+), 24(□), 25.5(○), 27(Δ), 28.5(∇), and 30(*) degrees while varying H from 12–80d, plotted against the relative height, utilizing Hstop(θ) as shown in (a) eqn (11) and (b) eqn (14). (c) Theoretical result when the local rheology is chosen consistent with the theoretical Hstop, per eqn (15). In (a)–(c), dashed line is the experimental collapse line, eqn (13), with β = 0.136.

5.1 A more precise Hstop curve

Collapse of the Froude number is sensitively dependent on the fit one uses for Hstop(θ). In the form obtained from our theory, we notice some deviations from the experiment which magnify, unsurprisingly, as one approaches the asymptote at θ1 (see Fig. 1). A different fit function,10 empirically matched to the data in Pouliquen,8 takes the form
image file: c4sm01838a-t11.tif(14)
where L0/d ≈ 1.65 (see the caption to Fig. 8 of Forterre and Pouliquen10). If we use the above relation for Hstop, a rather strong collapse of the Froude number versus relative height emerges (see Fig. 3(b)). Thus, it may be that the lack of a collapse of the Froude number is attributable to deviation of the predicted Hstop curve.

5.2 Consistent local rheology

An alternative explanation could be found if we revisit the history of the local inertial rheology. The rheology in eqn (3) was arrived at, coincidentally, by fitting data for glass beads flowing down inclines (see Appendix A of Jop et al.20). First, a form for Hstop was fitted empirically from experiments (i.e., eqn (14)). Then 〈ν〉 was calculated in terms of an as-yet-unknown function μloc(I) (assuming fully local rheology and, consequently, the Bagnold velocity profile, eqn (12)). This result was subsequently substituted, along with the empirical Hstop function, into eqn (13) to solve for μloc(I). The result is a local rheology having the form
image file: c4sm01838a-t12.tif(15)
which generates eqn (3), upon substituting eqn (14) for Hstop. As was done in Jop et al.,20 the term containing the square-root is replaced with a constant B (≈2.17 for glass beads) because the cos function does not vary much in the range of use. We see that, historically, the fit function chosen for Hstop gives rise to the functional form one obtains for the inertial rheology.

Let us repeat this logic in a thought experiment. Suppose, through eqn (15), we replace the local rheology in our theory with the one corresponding to our fit for the Hstop function, eqn (11). The result is

image file: c4sm01838a-t13.tif(16)
With this, the nonlocal fluidity dynamics are expressible conveniently as
image file: c4sm01838a-t14.tif(17)
Only the g2 term has changed so the previous stability argument identically produces our Hstop function, and eqn (16) is now obtained in the homogeneous flow limit as the new local response.

We apply the above system to inclined plane flows. Letting [z with combining tilde] = z/H and image file: c4sm01838a-t15.tif, writing P = ϕρsGz[thin space (1/6-em)]cos[thin space (1/6-em)]θ, recalling the definition of B, and allowing θ1θθ2, eqn (17) produces steady solutions given by

image file: c4sm01838a-t16.tif(18)
The only varying parameter in the above is the ratio H/Hstop. Hence, all solutions have the form, [g with combining tilde] = [g with combining tilde]([z with combining tilde]; H/Hstop). Likewise, the velocity field obeys
image file: c4sm01838a-t17.tif(19)
Finally, we arrive at a relation for the Froude number:
image file: c4sm01838a-t18.tif(20)
We conclude that when the same Hstop curve we derive theoretically is used to fit the local inertial relation, the emergent nonlocal theory requires a collapse between the Froude number and H/Hstop for all H > Hstop. This is an interesting result because while the local relation on its own would require such a collapse by definition, this is the full nonlocal solution. When H/Hstop is large enough, the second-derivative term in eqn (18) is negligible (excepting a rapid variation near the lower boundary), and the solution for [g with combining tilde] approaches that of a purely local material behavior. Likewise the function Fr(H/Hstop) obeys Fr(H/Hstop) ∼ βH/Hstop for HHstop, matching eqn (13). Numerical integration of eqn (17) verifies these analytical properties (see Fig. 3(c)). We observe that Fr(1) = 0, which is sensible as this corresponds to the stopping height, and our model does not include bistable terms. Although it is unclear in discrete simulations and experiments what the precise behavior of the Froude number is near H/Hstop = 1, for round grains on a fully rough surface many data trends show the Froude number exceeds zero at this height8 and other data suggests it equals zero,27 as our model purports.

A short commentary is in order. The above result comes at the cost of adjusting the μloc relation from the commonly used form of eqn (7), to a form compatible with the theory's own predicted Hstop function. The new relation, eqn (16), is qualitatively similar to the previous, with μ increasing monotonically from μs to μ2 as I increases. One difference is that the new relation has loc/dI = 0 at I = 0 whereas the former has a positive slope. There has been a debate recently over the behavior of μloc near I = 0. Some experimental data suggests, in fact, a negative initial slope for μloc,24,28 while other fits suggest a steeper-than-linear power-law behavior near I = 0.29 It bears noting that our previous work on nonlocal fluidity has focused solely on quasi-static flows, in which μ < μs almost everywhere and the only important aspect of the local rheology is the value of μs.

6 Conclusion

We have applied the theory of nonlocal granular fluidity to the canonical problem of size-dependence in granular inclined plane flows. The theory, as calibrated to glass beads based on prior data, predicts an Hstop curve that matches experimental data rather well. Further, the theory predicts flow profiles that vary in shape in a fashion consistent with existing discrete simulation data, marked by an upward curvature for layer thicknesses H close to Hstop and transitioning to the Bagnold profile for large H. While the Froude number does not exhibit a direct collapse against H/Hstop we have identified two possible explanations for this and corresponding ways in which the collapse can indeed be obtained.

Perhaps the most compelling result herein is the prediction of the Hstop function, which suggests the nonlocal fluidity concept could be used to model other size-sensitive flow stoppage phenomena, such as silo jamming. The famous Beverloo correlation,30 an empirical functional form that gives silo flow rate in terms of aperture opening size, indicates a critical opening size at which flow stops. This effect may be roughly analogous to the stopping height observed for inclined flows, and it would be useful to apply the nonlocal fluidity model in this geometry to determine if the Beverloo correlation can be obtained from the nonlocal model.


KK acknowledges funds from NSF-CBET-1253228 and the MIT Department of Mechanical Engineering, and DLH acknowledges funds from the Brown University School of Engineering.


  1. F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux and F. Chevoir, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 021309 CrossRef .
  2. P. Jop, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 032301 CrossRef .
  3. G. Koval, J.-N. Roux, A. Corfdir and F. Chevoir, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 021306 CrossRef .
  4. K. Kamrin, Int. J. Plast., 2010, 26, 167–188 CrossRef CAS PubMed .
  5. K. Nichol, A. Zanin, R. Bastien, E. Wandersman and M. van Hecke, Phys. Rev. Lett., 2010, 104, 078302 CrossRef .
  6. K. Reddy, Y. Forterre and O. Pouliquen, Phys. Rev. Lett., 2011, 106, 108301 CrossRef CAS .
  7. E. Wandersman and M. van Hecke, Europhys. Lett., 2014, 105, 24002 CrossRef .
  8. O. Pouliquen, Phys. Fluids, 1999, 11, 542–548 CrossRef CAS PubMed .
  9. L. E. Silbert, J. W. Landry and G. S. Grest, Phys. Fluids, 2003, 15, 1–10 CrossRef CAS PubMed .
  10. Y. Forterre and O. Pouliquen, J. Fluid Mech., 2003, 486, 21–50 CrossRef .
  11. G. D. R. MiDi, Eur. Phys. J. E, 2004, 14, 341–365 CrossRef CAS PubMed .
  12. I. S. Aranson and L. S. Tsimring, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 061303 CrossRef .
  13. L. Bocquet, J. Errami and T. C. Lubensky, Phys. Rev. Lett., 2002, 89, 184301 CrossRef .
  14. J. T. Jenkins, Phys. Fluids, 2006, 18, 103307 CrossRef PubMed .
  15. K. Kamrin and G. Koval, Phys. Rev. Lett., 2012, 108, 178301 CrossRef .
  16. D. L. Henann and K. Kamrin, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 6730–6735 CrossRef CAS PubMed .
  17. K. Kamrin and G. Koval, Comput. Part. Mech., 2014, 1, 169–176 CrossRef PubMed .
  18. D. L. Henann and K. Kamrin, Int. J. Plast., 2014, 60, 145–162 CrossRef PubMed .
  19. D. L. Henann and K. Kamrin, Phys. Rev. Lett., 2014, 113, 178001 CrossRef .
  20. P. Jop, Y. Forterre and O. Pouliquen, J. Fluid Mech., 2005, 541, 21–50 CrossRef .
  21. P. Jop, Y. Forterre and O. Pouliquen, Nature, 2006, 441, 727–730 CrossRef CAS PubMed .
  22. L. Bocquet, A. Colin and A. Ajdari, Phys. Rev. Lett., 2009, 103, 036001 CrossRef .
  23. P. Chaudhuri, V. Mansard, A. Colin and L. Bocquet, Phys. Rev. Lett., 2012, 109, 036001 CrossRef .
  24. J. A. Dijksman, G. H. Wortel, L. T. van Dellen, O. Dauchot and M. van Hecke, Phys. Rev. Lett., 2011, 107, 108303 CrossRef .
  25. A. Schoefield and P. Wroth, Critical State Soil Mechanics, McGraw-Hill, 1968 Search PubMed .
  26. O. Pouliquen, Phys. Rev. Lett., 2004, 93, 248001 CrossRef CAS .
  27. T. Weinhart, A. R. Thornton, S. Luding and O. Bokhove, Granular Matter, 2012, 14, 531–552 CrossRef CAS .
  28. O. Kuwano, R. Ando and T. Hatano, Geophys. Res. Lett., 2013, 40, 1295–1299 CrossRef .
  29. P.-E. Peyneau and J.-N. Roux, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 011307 CrossRef .
  30. W. Beverloo, H. Leniger and J. Van de Velde, Chem. Eng. Sci., 1961, 15, 260–269 CrossRef CAS .


In the numerics, we allow a small pressure on the top surface corresponding to the weight of a layer of (1/2)d thickness, i.e., Ptop = (1/2)ϕρsGd[thin space (1/6-em)]cos[thin space (1/6-em)]θ, to ensure that the coefficient of the g2 term in eqn (7) remains bounded.
In fact, the steady-state-only system corresponding to eqn (17) is given through eqn (4) with local rheology eqn (16) and cooperativity length eqn (8). We note that the coefficient of g2 in eqn (17) is for the moment only defined for μ > μs as needed for the Froude number analysis, but any monotonic extension of 1/Hstop for μ < μs may be assumed without affecting the analysis nor the corresponding steady-state-only relation.

This journal is © The Royal Society of Chemistry 2015