## Appendix A. Ill-posedness of the one-dimensional
$\unicode[STIX]{x1D707}(I),\unicode[STIX]{x1D6F7}(I)$
-rheology

Although there are some technical complications, effectively the dynamic equations of
$\unicode[STIX]{x1D707}(I),\unicode[STIX]{x1D6F7}(I)$
-rheology for two-dimensional flow are always ill-posed (Heyman *et al.*
2017). As in Barker *et al.* (2015), the ill-posedness has a directional character: plane-wave solutions (in all space) in certain directions suffer uncontrolled growth while those in other directions decay smoothly to uniform flow. If, as in § 3, solutions depending on only one spatial variable are sought, these one-dimensional equations may be either well-posed or ill-posed, depending on whether the one independent direction retained corresponds to one of the stable or unstable directions. In this appendix, conditions for the ill-posedness of the one-dimensional system (3.1)–(3.8) are derived.

The equations (3.1)–(3.3) together with (3.6) and (3.8) are linearized around a base flow
$(\unicode[STIX]{x1D719}^{0},u^{0},w^{0})$
in the dependent variables
$(\unicode[STIX]{x1D719},u,w)$
by introducing perturbations
$(\breve{\unicode[STIX]{x1D719}},\breve{u},\breve{w})$
of the form

and then freezing the base flow coefficients. In the linearization the highest-order derivatives of the perturbed variables are retained in each equation, together with the convective derivatives. This leads to linearized equations of the form

where
$a_{i}$
and
$b_{ij}$
are constant coefficients derived from the pressure (3.6) and the deviatoric stresses (3.8).

The linearized system (A 2) admits normal-mode solutions of the form

in which
$\unicode[STIX]{x1D709}$
is the wavenumber,
$\unicode[STIX]{x1D706}$
is the temporal growth rate and
$\tilde{\unicode[STIX]{x1D719}}$
,
$\tilde{u}$
and
$\tilde{w}$
are constants. Substituting these into (A 2) reveals that
$\unicode[STIX]{x1D706}$
must be an eigenvalue of the matrix

The imaginary terms
$-\text{i}w^{0}\unicode[STIX]{x1D709}$
on the diagonal shift the eigenvalues
$\unicode[STIX]{x1D706}$
, but do not affect stability or well-posedness, which are governed by the real part of
$\unicode[STIX]{x1D706}$
. The vertical component of the base state velocity
$w^{0}$
is therefore set to zero in what follows.

Denoting the bottom right
$2\times 2$
matrix of terms multiplying
$-\unicode[STIX]{x1D709}^{2}$
as

the eigenvalues are calculated by solving

where the coefficients are

To test for ill-posedness, the large-wavenumber limit
$\unicode[STIX]{x1D709}\rightarrow \infty$
is taken. Balancing the order of terms in (A 6), it is clear that
$\unicode[STIX]{x1D706}(\unicode[STIX]{x1D709})\sim \unicode[STIX]{x1D709}^{2}$
, so the substitution
$\unicode[STIX]{x1D706}=\unicode[STIX]{x1D6EC}\unicode[STIX]{x1D709}^{2}$
is made. Then the terms of maximal order in (A 6), i.e.
$O(\unicode[STIX]{x1D709}^{6})$
, have the coefficient
$\unicode[STIX]{x1D6EC}^{3}+(\text{tr}\,\unicode[STIX]{x1D71E})\unicode[STIX]{x1D6EC}^{2}+(\det \unicode[STIX]{x1D71E})\unicode[STIX]{x1D6EC}$
. Thus, to leading order, either
$\unicode[STIX]{x1D6EC}=0$
, or
$\unicode[STIX]{x1D6EC}$
is a solution of

i.e. one of the two roots

To determine the signs of these roots, the coefficients in (A 2) are now evaluated. If the superscript 0 notation is dropped, the coefficients in (A 2) are then

where all values are evaluated with the base-state fields

From (A 10) and (A 5) it follows that

and

The two roots in (A 9) are real since

Since
$\unicode[STIX]{x1D706}_{+}\sim \unicode[STIX]{x1D6EC}_{+}\unicode[STIX]{x1D709}^{2}$
as
$\unicode[STIX]{x1D709}\rightarrow \infty$
, the system is linearly ill-posed if the larger root
$\unicode[STIX]{x1D6EC}_{+}$
is positive, which occurs if either

It may be seen from (A 12) and (A 13) that if
$\text{tr}\,\unicode[STIX]{x1D71E}<0$
then
$\det \unicode[STIX]{x1D71E}<0$
, so it suffices to consider only case (b). Thus,

is a sufficient condition for ill-posedness.

Conversely, suppose

Then
$\text{tr}\,\unicode[STIX]{x1D71E}$
and
$\det \unicode[STIX]{x1D71E}$
are both positive, so
$\unicode[STIX]{x1D6EC}_{\pm }$
are both negative, and hence
$\unicode[STIX]{x1D706}_{\pm }\sim \unicode[STIX]{x1D6EC}_{\pm }\unicode[STIX]{x1D709}^{2}$
are bounded above. From (A 4), the third eigenvalue
$\unicode[STIX]{x1D706}_{0}$
asymptotes to
$-C/B$
as
$\unicode[STIX]{x1D709}\rightarrow \infty$
. But
$C/B$
is the ratio of two quartics and
$B\neq 0$
since
$\det \unicode[STIX]{x1D71E}\neq 0$
. Thus
$\unicode[STIX]{x1D706}_{0}$
also remains bounded as
$\unicode[STIX]{x1D709}\rightarrow \infty$
. In other words, (A 17) is a sufficient condition for well-posedness.

To rewrite (A 16), (A 17) in a way that is useful for the computations in § 3.2, let us define

The one-dimensional computations with
$\unicode[STIX]{x1D707}(I),\unicode[STIX]{x1D6F7}(I)$
-rheology will be

## Appendix B. Positive stress power and Onsager symmetry

Goddard & Lee (2018) examined equations for granular flow in light of Edelen’s (1972) theory of irreversibility. They showed that the model of Heyman *et al.* (2017), which includes
$\unicode[STIX]{x1D707}(I),\unicode[STIX]{x1D6F7}(I)$
-rheology, is not Onsager symmetric and not even properly dissipative. In this appendix the iCIDR equations are shown to have positive stress power and Onsager symmetry.

Positive-definite stress power means that all deformations with non-zero stress dissipate energy; i.e. in symbols, the dissipation rate
${\mathcal{D}}$
satisfies

Using the decompositions of the stress (2.3) and the strain rate (2.5)

and it follows that

Substituting the alignment condition (2.6) and the dilatancy law (4.2) and recalling that
$\Vert \unicode[STIX]{x1D64E}\Vert ^{2}=\text{tr}(\unicode[STIX]{x1D64E}^{2})/2$
, implies that

The yield condition (4.1), the yield function (4.7) and the dilatancy function (4.8) then imply that the dissipation rate is

where
$M=\unicode[STIX]{x1D707}(\unicode[STIX]{x1D6F9}(\unicode[STIX]{x1D719}))$
and
$\unicode[STIX]{x1D6F9}=\unicode[STIX]{x1D6F9}(\unicode[STIX]{x1D719})$
. Thus, the dissipation rate is non-negative and in fact positive if
$\Vert \unicode[STIX]{x1D64E}\Vert >0$
(and hence, by (4.12),
$p>0$
). To determine what happens when
$\Vert \unicode[STIX]{x1D64E}\Vert =0$
, the definition of the inertial number (2.9) is substituted into (B 5) to show that

As before, of course,
${\mathcal{D}}>0$
if
$\Vert \unicode[STIX]{x1D64E}\Vert >0$
. If
$\Vert \unicode[STIX]{x1D64E}\Vert =0$
but
$\text{div}\,\boldsymbol{u}<0$
, then by (4.12) the second term in (B 6) is positive. If
$\Vert \unicode[STIX]{x1D64E}\Vert =0$
but
$\text{div}\,\boldsymbol{u}>0$
, then indeed
${\mathcal{D}}=0$
, but
$p=0$
and by (4.14)
$\Vert \unicode[STIX]{x1D749}\Vert =0$
, so
$\Vert \unicode[STIX]{x1D748}\Vert =0$
. This completes the verification of (B 1), i.e. the iCIDR equations always dissipate energy provided there is some deformation and the stresses are non-zero.

Since the constitutive laws have positive stress power, it follows from Edelen (1972) that (i) all deformations are irreversible and increase entropy, (ii) the iCIDR constitutive laws satisfy Onsager symmetry and (iii) each component
$\unicode[STIX]{x1D70E}_{ij}$
of the stress equals the derivative of a convex dissipation potential
$\unicode[STIX]{x1D713}(\unicode[STIX]{x1D63F})$
with respect to
$\unicode[STIX]{x1D60B}_{ij}$
. Regarding point (iii), it is convenient to derive
$\unicode[STIX]{x1D713}$
through its Legendre transform,

where
$\unicode[STIX]{x1D6FF}$
is a shorter notation for
$\text{div}\,\boldsymbol{u}$
when it appears as a dummy variable. Consider the function

The claim is that the iCIDR constitutive relations can be expressed as the derivatives

The first relation implies

which, in light of (2.6), (4.1) and (4.7), verifies (B 9*a*
). Whilst the the second implies

which, in light of (4.2) and (4.8), verifies (B 9*b*
). The dissipation potential
$\unicode[STIX]{x1D713}(\unicode[STIX]{x1D63F})$
is the inverse Legendre transform of
$\unicode[STIX]{x1D711}$
, which may be written

where
$P(\unicode[STIX]{x1D64E},\text{div}\,\boldsymbol{u})$
is the pressure given by (4.12). Mixing stress and strain-rate variables,
$\unicode[STIX]{x1D713}$
can be calculated as a function of
$\unicode[STIX]{x1D64E}$
and
$p$
by replacing
$P$
in (B 12) by
$p$
and substituting (B 11) for
$\text{div}\,\boldsymbol{u}$
, which yields

Observe that

a result that follows from the formula (Goddard & Lee 2018)

and the observation that
${\mathcal{D}}$
is homogeneous in
$\unicode[STIX]{x1D63F}$
of degree 3. It is not proved directly that
$\unicode[STIX]{x1D713}$
is convex, but rather the fact that the iCIDR constitutive relations are well-posed (Barker *et al.*
2017) is used to deduce this (Goddard & Lee 2018).

## Appendix C. Details of the DEM calculations

Two-dimensional DEM simulations (Cundall & Strack 1979) were performed in a shear-box geometry, both to confirm the well-established steady behaviour and to explore transient flows. The domain is a rectangle in
$(x,z)$
space, with all units non-dimensionalized with the scalings (3.12), such that
$0\leqslant \hat{x}<\hat{L}$
and
$0\leqslant \hat{z}<{\hat{H}}$
. Boundary conditions and the system size are then designed to suppress confinement effects so that the volume can be taken to be a representative bulk element. Periodicity is enforced in the
$\hat{x}$
direction and Lees–Edwards boundary conditions (Lees & Edwards 1972) are applied at the bottom
$\hat{z}=0$
and top
$\hat{z}={\hat{H}}$
of the domain. There is no gravity applied to the system and the only driving is provided by the difference in horizontal velocity between the top and bottom
$V_{0}$
, as prescribed by the Lees–Edwards algorithm.

Details of the precise DEM simulation algorithm can be found in Otsuki & Hayakawa (2011) as an identical method is employed here. Normal forces
$\boldsymbol{f}^{(n)}$
between particles are calculated from a linear spring–dashpot arrangement with an associated spring constant
$k^{(n)}$
and viscous dissipation constant
$\unicode[STIX]{x1D702}^{(n)}$
. Tangential forces
$\boldsymbol{f}^{(t)}$
may either stick or slip depending on whether a Coulomb friction criterion, with particle friction constant
$\unicode[STIX]{x1D707}_{p}$
, is satisfied. Stick interactions are defined as those with
$|\boldsymbol{f}^{(t)}|<\unicode[STIX]{x1D707}_{p}|\boldsymbol{f}^{(n)}|$
and, like the normal force, are calculated from a linear spring–dashpot with parameters
$k^{(t)}$
and
$\unicode[STIX]{x1D702}^{(t)}$
. Interactions with greater computed tangential forces are labelled as slip events and the tangential force is truncated to
$|\boldsymbol{f}^{(t)}|=\unicode[STIX]{x1D707}_{p}|\boldsymbol{f}^{(n)}|$
. In this paper parameters are chosen which give
$k^{(n)}/p>10^{4}$
so that calculations are in the rigid particle regime of da Cruz *et al.* (2005). Unless stated otherwise the values
$\unicode[STIX]{x1D707}_{p}=0.4$
,
$k^{(n)}=10^{4}$
and
$k^{(t)}=0.5k^{(n)}$
are used. Particle interactions with these values are very short-lived so that results are invariant of the viscous dissipation. The tangential dashpot is therefore neglected (
$\unicode[STIX]{x1D702}^{(t)}=0$
) and
$\unicode[STIX]{x1D702}^{(n)}=4.2$
is chosen, as in Silbert *et al.* (2001), so that the particles have a restitution coefficient of
$e=0.9$
.

To select a mean packing fraction
$\unicode[STIX]{x1D719}_{0}$
, the system size along the
$\hat{x}$
direction is determined as
$L=\bar{A}N/(H\unicode[STIX]{x1D719}_{0})$
, where
$\bar{A}$
is the average grain area. The shear cell is then populated with
$N$
particles with density
$\unicode[STIX]{x1D70C}_{\ast }$
and mean diameter
$d$
and
$V_{0}$
set to unity in order to match the non-dimensional iCIDR equations. In order to avoid crystallization effects, the individual particle diameters are chosen randomly from a discretized distribution. Here an even spread from
$0.8d$
,
$0.9d$
,
$d$
,
$1.1d$
and
$1.2d$
is taken so that the number of particles of each diameter is
$N/5$
. These particles, which are initially elastic but not frictional, are randomly distributed in the domain. This results in overlaps which would normally cause very large elastic forces, so firstly there is a period during which the total kinetic energy of all particles is scaled to a small constant value so that the system can reach an equilibrium state. The arrangement which results from this procedure has almost uniform packing density and very small overlaps. Then, the interaction algorithm is altered so that the particles are approximately rigid (very large
$k^{(n)}$
) and frictional. The true simulation begins after the velocity fields are prescribed.

The steady-state
$\unicode[STIX]{x1D707}(I)$
and
$\unicode[STIX]{x1D6F7}(I)$
relations found in previous works (e.g. da Cruz *et al.*
2005) are first confirmed in order to obtain macroscopic rheological parameters. Both curves are derived from the same set of experiments in which the solids volume fractions
$\unicode[STIX]{x1D719}_{0}$
take values in the range
$0.76$
–
$0.8$
. This range is expected to lie in the inertial regime and, due to the
$\unicode[STIX]{x1D6F7}(I)$
relation, each packing corresponds to a unique inertial number. Once the system has reached a steady state, the flow fields are coarse-grained. This is achieved by averaging in the
$\hat{x}$
direction and then averaging within bins which discretize
$\hat{z}$
into boxes of height
$2d$
. Each run is then repeated 10 times in order to calculate error estimates. The solids volume fraction
$\unicode[STIX]{x1D719}$
and the two velocity components
$\hat{u}$
and
${\hat{w}}$
are clearly defined in the DEM data and allow the inertial number to be readily calculated from its definition (1.1). Calculation of the bulk friction coefficient
$\unicode[STIX]{x1D707}$
requires the macroscopic stress components to be defined. As in Silbert *et al.* (2001) this involves a sum over all particles
$\unicode[STIX]{x1D6FC}$
in the sampling volume:

where
$\boldsymbol{r}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$
is the centre-to-centre vector,
$\boldsymbol{f}^{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$
is the total force between particle pairs and
$\unicode[STIX]{x1D6FF}\hat{\boldsymbol{u}}^{\unicode[STIX]{x1D6FC}}=(\hat{u} ^{\unicode[STIX]{x1D6FC}}-\hat{V}_{0}\hat{z}/{\hat{H}},{\hat{w}}^{\unicode[STIX]{x1D6FC}})$
is the velocity fluctuation of particle
$\unicode[STIX]{x1D6FC}$
. From the stress decomposition (2.3) these stress components are used to calculate the pressure
$p$
, deviatoric stress
$\Vert \unicode[STIX]{x1D749}\Vert$
and hence
$\unicode[STIX]{x1D707}$
across the domain. As these quantities are all found to be invariant of
$\hat{z}$
, the mean value is taken for the
$\unicode[STIX]{x1D707}(I)$
data presented in figure 1(*a*). The volume fraction is also found to be constant at steady state so that the
$\unicode[STIX]{x1D6F7}(I)$
relation, plotted in figure 1(*b*), is simply verified using the mean packing density
$\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D719}_{0}$
.

## References

Bagnold, R. A.
1954
Experiments on a gravity-free dispersion of large solid spheres in a Newtonian fluid under shear. Proc. R. Soc. Lond. A
225 (1160), 49–63.

Barker, T. & Gray, J.
2017
Partial regularisation of the incompressible 𝜇(*I*)-rheology for granular flow. J. Fluid Mech.
828, 5–32.

Barker, T., Schaeffer, D. G., Bohorquez, P. & Gray, J. M. N. T.
2015
Well-posed and ill-posed behaviour of the 𝜇(*I*)-rheology for granular flow. J. Fluid Mech.
779, 794–818.

Barker, T., Schaeffer, D. G., Shearer, M. & Gray, J. M. N. T.
2017
Well-posed continuum equations for granular flow with compressibility and 𝜇(*I*)-rheology. Proc. R. Soc. Lond. A
473 (2201), 20160846.

Bouzid, M., Trulsson, M., Claudin, P., Clément, E. & Andreotti, B.
2013
Nonlocal rheology of granular flows across yield conditions. Phys. Rev. Lett.
111 (23), 238301.

Canuto, C., Hussaini, M. Y., Quarteroni, A. & Zang, T. A.
1988
Spectral Methods in Fluid Dynamics, Springer Series in Computational Physics. Springer.

Chialvo, S., Sun, J. & Sundaresan, S.
2012
Bridging the rheology of granular flows in three regimes. Phys. Rev. E
85 (2), 021305.

da Cruz, F., Emam, S., Prochnow, M., Roux, J. & Chevoir, F.
2005
Rheophysics of dense granular materials: discrete simulation of plane shear flows. Phys. Rev. E
72, 021309.

Cundall, P. A. & Strack, O. D.
1979
A discrete numerical model for granular assemblies. Geotechnique
29 (1), 47–65.

Daerr, A. & Douady, S.
1999
Two types of avalanche behaviour in granular media. Nature
399, 241–243.

DeGiuli, E. & Wyart, M.
2017
Friction law and hysteresis in granular materials. Proc. Natl Acad. Sci. USA
114, 9284–9289.

Edelen, D. G.
1972
A nonlinear Onsager theory of irreversibility. Intl J. Engng Sci.
10 (6), 481–490.

Edwards, A. N. & Gray, J. M. N. T.
2015
Erosion–deposition waves in shallow granular free-surface flows. J. Fluid Mech.
762, 35–67.

Edwards, A. N., Viroulet, S., Kokelaar, B. P. & Gray, J. M. N. T.
2017
Formation of levees, troughs and elevated channels by avalanches on erodible slopes. J. Fluid Mech.
823, 278–315.

Forterre, Y. & Pouliquen, O.
2008
Flows of dense granular media. Annu. Rev. Fluid Mech.
40, 1–24.

GDR MiDi
2004
On dense granular flows. Eur. Phys. J. E
14 (4), 341–365.

Goddard, J. D. & Lee, J.
2018
Regularization by compressibility of the 𝜇(*I*) model of dense granular flow. Phys. Fluids
30, 073302.

Gray, J. M. N. T. & Edwards, A. N.
2014
A depth-averaged 𝜇(*I*)-rheology for shallow granular free-surface flows. J. Fluid Mech.
755, 503–534.

Heyman, J., Delannay, R., Tabuteau, H. & Valance, A.
2017
Compressibility regularizes the 𝜇(*I*)-rheology for dense granular flows. J. Fluid Mech.
830, 553–568.

Jackson, R.
1983
Some mathematical and physical aspects of continuum models for the motion of the granular materials. In Theory of Dispersed Multiphase Flow (ed. Meyer, R.). Academic Press.

Jop, P., Forterre, Y. & Pouliquen, O.
2006
A constitutive relation for dense granular flows. Nature
44, 727–730.

Kamrin, K. & Koval, G.
2012
Nonlocal constitutive relation for steady granular flow. Phys. Rev. Lett.
108 (17), 178301.

Lagrée, P.-Y., Staron, L. & Popinet, S.
2011
The granular column collapse as a continuum: validity of a two-dimensional Navier–Stokes model with a 𝜇(*I*)-rheology. J. Fluid Mech.
686, 378–408.

Lees, A. & Edwards, S.
1972
The computer study of transport processes under extreme conditions. J. Phys. C: Solid State Phys.
5 (15), 1921–1928.

Liu, A. & Nagel, S.
1998
Jamming is not just cool any more. Nature
396, 21–22.

Martin, N., Ionescu, I. R., Mangeney, A., Bouchut, F. & Farin, M.
2017
Continuum viscoplastic simulation of a granular column collapse on large slopes: 𝜇(*I*) rheology and lateral wall effects. Phys. Fluids
29, 013301.

Otsuki, M. & Hayakawa, H.
2009
Critical behaviors of sheared frictionless granular materials near the jamming transition. Phys. Rev. E
80 (1), 011308.

Otsuki, M. & Hayakawa, H.
2011
Critical scaling near jamming transition for frictional granular particles. Phys. Rev. E
83 (5), 051301.

Pouliquen, O., Cassar, C., Jop, P., Forterre, Y. & Nicolas, N.
2006
Flow of dense granular material: towards simple constitutive laws. J. Stat. Mech.
2006, P07020.

Pouliquen, O. & Forterre, Y.
2002
Friction law for dense granular flows: application to the motion of a mass down a rough inclined plane. J. Fluid Mech.
453, 133–151.

Pouliquen, O. & Forterre, Y.
2009
A non-local rheology for dense granular flows. Phil. Trans. R. Soc. Lond. A
367, 5091–5107.

Russell, A. S., Johnson, C. G., Edwards, A. N., Viroulet, S., Rocha, F. M. & Gray, J. M. N. T.
2019
Retrogressive failure of a static granular layer on an inclined plane. J. Fluid Mech.
869, 313–340.

Schiesser, W. E.
2012
The Numerical Method of Lines: Integration of Partial Differential Equations. Elsevier.

Schofield, A. & Wroth, P.
1968
Critical State Soil Mechanics, vol. 310. McGraw-Hill.

Silbert, L. E., Ertaş, D., Grest, G. S., Halsey, T. C., Levine, D. & Plimpton, S. J.
2001
Granular flow down an inclined plane: Bagnold scaling and rheology. Phys. Rev. E
64 (5), 051302.

Singh, A., Magnanimo, V., Saitoh, K. & Luding, S.
2015
The role of gravity or pressure and contact stiffness in granular rheology. New J. Phys.
17 (4), 043028.

Staron, L., Lagrée, P.-Y. & Popinet, S.
2012
The granular silo as a continuum plastic flow: the hour-glass versus the clepsydra. Phys. Fluids
24, 103301.

Trefethen, L. N.
2000
Spectral Methods in MATLAB, Software, Environments, and Tools, vol. 10. Society for Industrial and Applied Mathematics (SIAM).

Trulsson, M., Bouzid, M., Claudin, P. & Andreotti, B.
2013
Dynamic compressibility of dense granular shear flows. Eur. Phys. Lett.
103, 38002.