Nonlocal modeling of granular flows down inclines

Abstract

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;¹ 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,¹² 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 MiDi¹¹) 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 Koval¹⁷). The inertial rheology for steady flow of dense granular media relates the local stress state and the strain-rate^1,20,21 and is expressed via the dimensionless relationship: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 μ₂. A common fit is²⁰

μ_loc (I) = μ_s + Δμ/(I₀/I + 1) (3)

where Δμ ≡ μ₂ − μ_s and I₀ ≡ Δμ/b.

The inertial rheology works well in describing uniform flows (e.g., planar shearing) over a wide range of flow rates.¹ 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.³ 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 simulations¹⁵ and experiments.¹⁶ It is the first continuum model to quantitatively describe all experimentally-obtained flow data in the complex split-bottom family of geometries.¹⁶ 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.¹⁹

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

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 satisfiesroughly in-line with prediction of the kinetic-elasto-plastic (KEP) theory on which other fluidity models are based.²² 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,²³ which is encouraging for our current study.

2.1 Dynamical system for fluidity

In Henann and Kamrin,¹⁸ 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:where t₀ > 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 = g_loc + δ for some small function δ, where g_loc 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 givingWhen 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) isNote that the limit of eqn (8) as μ₂ 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.²⁴ To the latter point, Reynolds dilation during shear initiation induces transient variations in material strength, which are commonly described using critical-state models²⁵ 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.²² The dynamical form of the nonlocal fluidity model shows mathematical similarities to order-parameter-based rheological approaches, which also feature a diffusing state variable.¹²

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 θ, P = ϕρ_sGz cos θ. (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 μ_s, i.e., any layer tilted above θ_r would be predicted to flow.

Fig. 1 (a) Setup of the inclined plane flow geometry. (b) Theoretically predicted H_stop locus (−), eqn (11), as calibrated for glass beads on a fully rough base, compared to experimentally determined values (○) from Pouliquen.⁸

A signature of the cooperativity of granular motion is the fact that this universal repose angle is contradicted in experiments. As shown initially by Pouliquen⁸ 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 H_stop (θ) 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 = g_loc + δ 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 g_loc 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,⁸ where glass beads were used as the media. The local continuum parameters for glass beads have been calibrated in Jop et al.²⁰ and are μ_s = 0.382 = tan(20.9°), μ₂ = 0.643 = tan(32.76°), b = Δμ/I₀ = 0.938, and ρ_s = 2450 kg m⁻³. The nonlocal amplitude for glass beads was calibrated in Henann and Kamrin¹⁶ 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.,⁹ 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 g² term negligible. The remaining PDE is linear and, substituting μ = tan θ, it is solved by

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 ifA similar mathematical technique was utilized in Aranson and Tsimring.¹² The above shows that for tall layers, H ≫ O(d), the material is predicted to have a common repose angle of θ₁ = 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 θ₂ = arctan(μ₂) 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 H_stop data in inclined plane flows in multiple experiments and particle simulations involving different materials and substrates.¹¹ When it comes to specific functional forms for H_stop data, several fit curves have been proposed, a topic we shall return to later, but all invoke critical θ₁ and θ₂ values as described here. It bears noting that the derived H_stop relation, eqn (11), is linearly proportional to the cooperativity length, eqn (8), i.e., H_stop(arctan μ) = (π/2)ξ(μ). This is in line with experimental observations of Pouliquen,²⁶ who reported that velocity correlation lengths in inclined plane flows vary with θ (and hence μ) in a similar manner as H_stop. 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 H_stop(θ), 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.¹⁶ 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 H_start(θ) – 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 H_stop, though it shares the same qualitative features (e.g., θ₁ and θ₂ values where the curve asymptotes and vanishes respectively). Our analysis above is tailored to H_stop 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 H_start curve does not arise theoretically.

4 Velocity profiles

For angles between θ₁ and θ₂, and H exceeding H_stop, the material has a steady flow state down the incline. As reported collectively in GDR MiDi¹¹ and explored directly in Silbert et al.,⁹ 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) ∝ H^3/2 − z^3/2. (12)

Fig. 2(a) shows the discrete particle simulation data of Silbert⁹ for inclination angle θ = 24° and various values of H > H_stop(θ = 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 H_stop/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/H_stop 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 Δx ≪ d and the Bagnold profile as the initial guess.†

Fig. 2 (a) Discrete particle simulation data of Silbert et al.⁹ 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/H_stop ≅ 1.1, 1.9, 2.7, 4.1, 6.4, as well as H/H_stop = 20.

The model shows that for H near H_stop 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 g² 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 H_stop, 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 H_stop proportionally, it is the value of H/H_stop 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,⁸ and others thereafter,^9–11 have observed that the Froude number of the flow, , appears to be a relatively well-defined, one-to-one function of the relative height, H/H_stop, for all angles and heights as long as H is not close to H_stop. In this range, they findwhere, 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.

Fig. 3 Predicted Froude number, , 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 H_stop(θ) as shown in (a) eqn (11) and (b) eqn (14). (c) Theoretical result when the local rheology is chosen consistent with the theoretical H_stop, per eqn (15). In (a)–(c), dashed line is the experimental collapse line, eqn (13), with β = 0.136.

5.1 A more precise H_stop curve

Collapse of the Froude number is sensitively dependent on the fit one uses for H_stop(θ). In the form obtained from our theory, we notice some deviations from the experiment which magnify, unsurprisingly, as one approaches the asymptote at θ₁ (see Fig. 1). A different fit function,¹⁰ empirically matched to the data in Pouliquen,⁸ takes the formwhere L₀/d ≈ 1.65 (see the caption to Fig. 8 of Forterre and Pouliquen¹⁰). If we use the above relation for H_stop, 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 H_stop 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.²⁰). First, a form for H_stop 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 H_stop function, into eqn (13) to solve for μ_loc(I). The result is a local rheology having the formwhich generates eqn (3), upon substituting eqn (14) for H_stop. As was done in Jop et al.,²⁰ 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 H_stop 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 H_stop function, eqn (11). The result is

With this, the nonlocal fluidity dynamics are expressible conveniently asOnly the g² term has changed so the previous stability argument identically produces our H_stop 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̃ = z/H and , writing P = ϕρ_sGz cos θ, recalling the definition of B, and allowing θ₁ ≤ θ ≤ θ₂, eqn (17) produces steady solutions given by

The only varying parameter in the above is the ratio H/H_stop. Hence, all solutions have the form, g̃ = g̃(z̃; H/H_stop). Likewise, the velocity field obeysFinally, we arrive at a relation for the Froude number:We conclude that when the same H_stop 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/H_stop for all H > H_stop. 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/H_stop 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̃ approaches that of a purely local material behavior. Likewise the function Fr(H/H_stop) obeys Fr(H/H_stop) ∼ βH/H_stop for H ≫ H_stop, 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/H_stop = 1, for round grains on a fully rough surface many data trends show the Froude number exceeds zero at this height⁸ and other data suggests it equals zero,²⁷ 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 H_stop function. The new relation, eqn (16), is qualitatively similar to the previous, with μ increasing monotonically from μ_s to μ₂ as I increases. One difference is that the new relation has dμ_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.²⁹ 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 H_stop 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 H_stop and transitioning to the Bagnold profile for large H. While the Froude number does not exhibit a direct collapse against H/H_stop 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 H_stop 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,³⁰ 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.

Acknowledgements

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.

References

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

---

Footnotes

† 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., P_top = (1/2)ϕρ_sGd cos θ, to ensure that the coefficient of the g² 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 g² in eqn (17) is for the moment only defined for μ > μ_s as needed for the Froude number analysis, but any monotonic extension of 1/H_stop for μ < μ_s may be assumed without affecting the analysis nor the corresponding steady-state-only relation.

---