Local Rheology Relation with Variable Yield Stress Ratio across Dry, Wet, Dense, and Dilute Granular Flows

TL;DR

This study uses extensive simulations to unify the rheology of various granular flows, revealing a common scaling law that breaks down near flow thresholds, leading to a variable yield stress ratio.

Contribution

It introduces a generalized rheology relation that applies across dry, wet, dense, and dilute granular flows, accounting for variable yield stress ratios.

Findings

Mohr-Coulomb friction coefficient scales with the square root of Péclet number in most conditions.

Scaling breaks down at low Péclet number and high temperature gradients, indicating a variable yield stress ratio.

The results unify diverse granular flow behaviors under a common rheological framework.

Abstract

Dry, wet, dense, and dilute granular flows have been previously considered fundamentally different and thus described by distinct, and in many cases incompatible, rheologies. We carry out extensive simulations of granular flows, including wet and dry conditions, various geometries and driving mechanisms (boundary driven, fluid driven, and gravity driven), many of which are not captured by standard rheology models. For all simulated conditions, except for fluid-driven and gravity-driven flows close to the flow threshold, we find that the Mohr-Coulomb friction coefficient $\mu$ scales with the square root of the local P\'eclet number $\mathrm{Pe}$ provided that the particle diameter exceeds the particle mean free path. With decreasing $\mathrm{Pe}$ and granular temperature gradient $M$, this general scaling breaks down, leading to a yield condition with a variable yield stress ratio characterized by $M$.

Figures

Tables

Table 1. Table 1: Summary of simulated granular flows.

Flow geometry	Driven by	Contact model ($e$, ${\mu }_c$)
Sediment transport (2D)	Fluid	Linear (0.9, 0.5)
Rapid gravity flows (2D)	Gravity	Linear (0.9, 0.5)
Sheared suspensions (2D)	Boundary	Linear (0.1, 0.4)
Dry shear flows (2D)	Boundary	Linear (0.1, 0.4)
Lubricated drum flows (3D)	Boundary	Hertz (0.5, 0.5)
Dry drum flow (3D)	Boundary	Hertz (0.5, 0.5)

Table 2. Table 2: Sediment transport regimes  Pähtz and Durán ( 2018a ) .

Sediment transport regime	Condition
Viscous bedload transport
Turbulent bedload transport
Bedload-saltation transition
Viscous saltation transport
Turbulent saltation transport	$\sqrt{s}\mathrm{Ga}\ge 80\wedge s\ge 10\wedge \sqrt[4]{s}\mathrm{Ga}\ge 32$

Table 3. Table S1: Summary of the simulated sediment transport and gravity flow conditions included in Figs. 2, S2, and S4.

Category	Density ratio $s$	Galileo number $\mathrm{Ga}$	Range of $\sqrt{\Theta }$ or $\alpha$	# of simulations
Viscous bedload transport	$2.65$	$2$	$\sqrt{\Theta }\in [0.79,1.39]$	$11$
Viscous bedload transport	$2.65$	$5$	$\sqrt{\Theta }\in [0.73,1.18]$	$10$
Viscous bedload transport	$2.65$	$10$	$\sqrt{\Theta }\in [0.77,1.12]$	$8$
Viscous bedload transport	$100$	$0.5$	$\sqrt{\Theta }\in [0.84,1.02]$	$4$
Viscous bedload transport	$2000$	$0.1$	$\sqrt{\Theta }\in [0.84,1.05]$	$4$
Turbulent bedload transport	$2.65$	$20$	$\sqrt{\Theta }\in [0.67,0.91]$	$7$
Turbulent bedload transport	$2.65$	$50$	$\sqrt{\Theta }\in [0.61,1.45]$	$19$
Turbulent bedload transport	$2.65$	$100$	$\sqrt{\Theta }\in [0.6,1.1]$	$12$
Bedload-saltation transition	$100$	$2$	$\sqrt{\Theta }\in [0.52,0.7]$	$8$
Bedload-saltation transition	$100$	$5$	$\sqrt{\Theta }\in [0.4,0.6]$	$5$
Bedload-saltation transition	$2000$	$0.5$	$\sqrt{\Theta }\in [0.73,0.9]$	$3$
Viscous saltation transport	$100$	$10$	$\sqrt{\Theta }\in [0.34,0.64]$	$11$
Viscous saltation transport	$2000$	$2$	$\sqrt{\Theta }\in [0.34,0.72]$	$7$
Turbulent saltation transport	$100$	$20$	$\sqrt{\Theta }\in [0.35,0.65]$	$4$
Turbulent saltation transport	$100$	$50$	$\sqrt{\Theta }\in [0.28,0.31]$	$2$
Turbulent saltation transport	$100$	$100$	$\sqrt{\Theta }\in [0.28,0.31]$	$2$
Turbulent saltation transport	$2000$	$5$	$\sqrt{\Theta }\in [0.3,0.6]$	$4$
Turbulent saltation transport	$2000$	$10$	$\sqrt{\Theta }\in [0.35,0.45]$	$2$
Gravity flow in ambient air	$2000$	$2$	$\alpha \in [{51}^{\circ },{60}^{\circ }]$	$4$
Gravity flow in ambient air	$2000$	$5$	$\alpha \in [{48}^{\circ },{60}^{\circ }]$	$5$
Gravity flow in ambient air	$2000$	$10$	$\alpha \in [{42}^{\circ },{60}^{\circ }]$	$7$
Gravity flow in ambient air	$2000$	$20$	$\alpha \in [{36}^{\circ },{60}^{\circ }]$	$9$
Gravity flow in ambient air	$2000$	$50$	$\alpha \in [{33}^{\circ },{60}^{\circ }]$	$10$
Gravity flow in ambient air	$2000$	$100$	$\alpha \in [{30}^{\circ },{60}^{\circ }]$	$11$

Equations

$$\mu =c_{\mu }\sqrt{\mathrm{Pe}},$$

Full text

Reliable large-scale simulations and thus predictions of geophysical and industrial processes require a deep understanding of the continuum properties of granular flows. However, existing theories of the granular flow rheology are limited to small subsets of the physical conditions under which such processes can occur. For example, although geophysical granular flows are often wet (i.e., significantly affected or driven by ambient fluid) Courrech du Pont et al. (2003); Houssais and Jerolmack (2017); Delannay et al. (2017) and consist of coexisting dense (liquidlike) and dilute (gaslike) flow layers Börzsönyi et al. (2009); Holyoake and McElwaine (2012); Brodu et al. (2013, 2015); Delannay et al. (2017), even understanding comparably simple dry, dense-only or dilute-only flows has remained a major challenge MiDi (2004); Andreotti et al. (2013); Jop (2015); Kumaran (2015).

Existing rheologies for noncohesive and nonquasistatic flows of sufficiently hard granular particles can be classified in terms of the particle volume fraction $\phi$ (the fraction of space covered by particles), particle-fluid-density ratio $s\equiv\rho_{p}/\rho_{f}$, and Stokes number $\mathrm{St}\equiv\rho_{p}\dot{\gamma}d^{2}/\eta_{f}$, where $d$ is the particle diameter, $\dot{\gamma}$ the granular shear rate, and $\eta_{f}$ the fluid viscosity. Dilute, dry flows ($\phi\lesssim 0.5$, $s\rightarrow\infty$, $\mathrm{St}\rightarrow\infty$) have been described by the kinetic theory of dry granular gases Sela and Goldhirsch (1998); Saha and Alam (2014, 2016); Garzó and Dufty (1999), dense, dry flows ($\phi\gtrsim 0.5$, $s\rightarrow\infty$, $\mathrm{St}\rightarrow\infty$) by the local viscoplastic rheology Jop et al. (2006) and its nonlocal extensions Kamrin and Koval (2012); Bouzid et al. (2013, 2015); Zhang and Kamrin (2017), dense solid-liquid suspensions ($\phi\gtrsim 0.5$, $s\simeq 1$, variable $\mathrm{St}$) by different (partially incompatible Guazzelli and Pouliquen (2018)) viscoinertial rheologies Boyer et al. (2011); Trulsson et al. (2012); Ness and Sun (2015, 2016); Amarsid et al. (2017); DeGiuli et al. (2015), and sediment transport driven by liquids (variable $\phi$, $1<s<3$, variable $\mathrm{St}$) by modified viscoplastic or viscoinertial rheologies Houssais et al. (2015, 2016); Maurin et al. (2016). Furthermore, different rather complex and controversial approaches exist to extend kinetic theory to solid-gas suspensions Garzó et al. (2012); Chamorro et al. (2015); Saha and Alam (2017); Alam et al. (2019) or the dense regime Chialvo and Sundaresan (2013); Vescovi et al. (2014); Berzi and Vescovi (2015).

Here we show that, despite their fundamental differences, granular flows from the entire phase space $(\phi,s,\mathrm{St})$ actually obey a common scaling law for the Mohr-Coulomb friction coefficient $\mu$, the knowledge of which is essential for any rheological description.

We carry out discrete element method-based simulations of granular flows for a variety of geometries and driving mechanisms (Table 1 and Fig. 1), which cover the entire phase space: (i) two-dimensional sediment transport driven by a large variety of Newtonian fluids; (ii) two-dimensional rapid gravity-driven flows in ambient static air of varying viscosity, many of which are highly convective (e.g., they can exhibit a strong kinetic heat transfer normal to the flow direction) and/or “supported” Brodu et al. (2015); (iii) two-dimensional uniformly sheared viscous suspensions in density-matched fluid of varying viscosity; (iv) two-dimensional dry uniform shear flows; (v) three-dimensional rotating drum flows lubricated by a density-matched fluid; and (vi) a three-dimensional dry rotating drum flow. Among these flows, rapid gravity flows and rotating drum flows are known to elude the description by standard rheology models Cortet et al. (2009); Börzsönyi et al. (2009); Holyoake and McElwaine (2012); Brodu et al. (2013, 2015); Govender (2016). In all simulations, contacting particles interact via normal repulsion (restitution coefficient $e$, modeled through viscous damping), governed either by a linear or Hertzian law, and tangential friction (contact friction coefficient $\mu_{c}$, Table 1). Details are described below.

*Sediment transport and gravity flows.—*The numerical model couples a discrete element method for the particle motion (stiffness $k=5000mg/d$) under gravity, buoyancy, and fluid drag with a continuum Reynolds-averaged description of hydrodynamics (described in detail and/or validated in Refs. Durán et al. (2011, 2012, 2014); Pähtz and Durán (2017, 2018a, 2018b)). Spherical particles ($\sim 10^{4}$) with mild polydispersity are confined in a quasi-two-dimensional, vertically infinite domain of length $\sim 10^{3}d$ with periodic boundary conditions in the flow direction. For gravity flows, the ambient fluid is kept static.

Simulations are carried out for varying density ratio $s$, Galileo number $\mathrm{Ga}\equiv\rho_{f}\sqrt{(s-1)gd^{3}}/\eta_{f}$, Shields number $\Theta\equiv\tau_{f}/[(\rho_{p}-\rho_{f})gd]$, and inclination angle $\alpha$, where $g$ is the gravitational constant and $\tau_{f}$ is the bed fluid shear stress. For gravity flows, we simulate conditions with $s=2000$, $\mathrm{Ga}\in[2,100]$, $\Theta=0$, and $\alpha$ between the flow threshold and $60^{\circ}$. For sediment transport, we simulate conditions with $s\in[2.65,2000]$, $\mathrm{Ga}\in[0.1,100]$, $\alpha=0$, and $\Theta$ above the flow threshold, which correspond to five different transport regimes (Table 2) Pähtz and Durán (2018a). Following the symmetry along the flow direction, simulation data are averaged over horizontal layers of variable thickness depending on the particle volume fraction Durán et al. (2012).

*Uniformly sheared particle and suspension flows.—*The numerical model couples a discrete element method for the particle motion ($k=2000P_{zz}d$) under viscous fluid drag and torque with the Stokes equations for laminar flow (described in detail in Refs. Trulsson et al. (2012, 2017)). Two-dimensional disks ($\sim 10^{3}$) with moderate polydispersity are confined within a shear cell composed by two rough walls, created by gluing together two dense layers of grains, with periodic boundary conditions along the flow direction parallel to the walls. The position of the walls is controlled to ensure constant confining pressure $P_{zz}$ and mean shear rate.

Simulations are carried out for varying volume fraction (in the range $\phi>0.24$, where $\phi$ is calculated as $2/3$ of the disk area fraction, like for spheres confined in two dimensions) and two general cases: no ambient fluid (dry condition) and an ambient density-matched liquid with varying dimensionless viscosity ($s=1$, $\eta_{f}/\sqrt{\rho_{p}P_{zz}d^{2}}=[10^{-3},10^{-2},10^{-1},10^{0},\infty]$).

*Rotating drum flows.—*The numerical model uses a discrete element method for the particle motion ($k=17000mg/d^{3/2}$) under lubrication forces Ness and Sun (2015) and gravity. The contact model employs the LIGGGHTS implementation of Hertzian contacts, which ensures a constant value of $e$ Di Renzo and Di Maio (2004, 2005). Spherical monodisperse particles ($\sim 10^{4}$) are confined within a closed horizontal cylinder (drum) of radius $20d$ and width $20d$ rotating at a constant rate $\omega$.

Simulations are carried out for dry conditions and an ambient density-matched liquid with varying dimensionless viscosity ($s=1$, $\eta_{f}/(\rho_{f}\omega d^{2})=[1/160,1/16,3/16]$). Simulation data are averaged using an anisotropic Gaussian smoothing function of dimension $3d\times 3d\times 20d$.

*General rheology relation.—*Contact ($\mathbf{P}^{c}$) and kinetic ($\mathbf{P}^{k}$) granular stresses are calculated from the simulation data using the method given in Ref. Artoni and Richard (2015), which ensures that the granular temperature $T=\sum_{i}P^{k}_{ii}/(\rho_{p}\phi D)$ (the root-mean-square of the particle fluctuation velocity), where $D$ is the number of space dimensions, is insensitive to the coarse-graining width. Furthermore, the shear rate $\dot{\gamma}$ is calculated as the norm of the deviatoric component of the strain rate tensor $\epsilon_{ij}\equiv\partial_{i}\langle v_{j}\rangle+\partial_{j}\langle v_{i}\rangle$, which reads $\dot{\gamma}\equiv\lVert\boldsymbol{\epsilon}^{d}\rVert=\sqrt{\sum_{ij}\epsilon^{d}_{ij}\epsilon^{d}_{ij}/2}$, where $\epsilon^{d}_{ij}=\epsilon_{ij}-\sum_{k}\epsilon_{kk}\delta_{ij}/D$. Finally, we calculate the Mohr-Coulomb friction coefficient $\mu\equiv\max_{i,j}|P^{c}_{i}-P^{c}_{j}|/|P^{c}_{i}+P^{c}_{j}|$ from the principal components $P^{c}_{i}$ of $\mathbf{P}^{c}$ (a Drucker-Prager definition of $\mu$ and/or a definition that includes kinetic stresses yield slightly but statistically significantly worse results).

For both dilute and dense flow conditions (defined shortly), we find that the friction coefficient $\mu$ is only a function of the Péclet number $\mathrm{Pe}\equiv\dot{\gamma}d/\sqrt{T}$ Chialvo and Sundaresan (2013) and scales as

[TABLE]

almost everywhere within all simulated conditions [Fig. 2(a), closed symbols], except for sediment transport and gravity flows too close to the flow threshold (which are excluded from Fig. 2) because of nonlocal effects (see Supplementary Material for details Sup ). The scaling parameter $c_{\mu}$ slightly depends on the driving conditions, with the smallest value found for viscous bedload transport. For uniformly sheared flows, $c_{\mu}$ varies significantly with $\mu_{c}$ ($c_{\mu}\simeq 0.3$ for $\mu_{c}=0$, $c_{\mu}\simeq 0.55$ for $\mu_{c}=100$) but not with $e$, even in the extreme cases $\mu_{c}=100$ (no sliding) and $\mu_{c}=0$ (always sliding) Sup . We therefore propose that Eq. (1) relates the structure anisotropy to $\mu$, where $c_{\mu}$ partly encodes the anisotropy of the tangential contact force network (which increases with $\mu_{c}$) and $\mathrm{Pe}$ encodes the anisotropy of the particle assembly (i.e., the unit normal contact vector), relating the diffusion forces toward isotropic configurations to that of the anisotropic compression and extension of the shear Chialvo and Sundaresan (2013). This proposition is consistent with dry uniform shear flows, for which both anisotropies correlate with $\mu$ Azéma and Radjaï (2014). Alternatively, $\mathrm{Pe}$ may also be interpreted as the ratio between the rates of the macroscopic shearing motion ($\dot{\gamma}$) and microscopic kinetic rearrangements ($\propto\sqrt{T}/d$), which is similar to the original interpretations of the inertial number $I\equiv\dot{\gamma}d/\sqrt{P^{c}/\rho_{p}}$ (where $P^{c}\equiv\sum_{i}P^{c}_{ii}/D$) and viscous number $I_{\nu}\equiv\eta_{f}\dot{\gamma}/P^{c}$ Cassar et al. (2005). The difference is that, for $I$ (and $I_{\nu}$), the kinetic rearrangement rate has been obtained from assuming a particle fall driven by pressure (and opposed by viscous drag), whereas for $\mathrm{Pe}$, the kinetic rearrangement rate is obtained from the actual relative motion between neighboring particles. The latter rate should be more general, which is supported by the fact that, in contrast to $\mathrm{Pe}$, neither $I$ (even when limited to dry flows), nor $I_{\nu}$, nor a combination of the two collapse the $\mu$ data [Fig. 2(b) and Supplementary Material Sup ], except the viscoinertial number $\sqrt{I^{2}+5I_{\nu}}$ for uniformly sheared flows [Fig. 2(b)]. Note that a standard nonlocal rheology model Kamrin and Koval (2012); Zhang and Kamrin (2017) and extended kinetic theory also fail to describe our flows Sup . In particular, in Navier-Stokes order, the latter predicts $\mu\simeq c^{\mathrm{kin}}_{\mu}\mathrm{Pe}$ for dense flows Chialvo and Sundaresan (2013); Berzi and Vescovi (2015) (with a proportionality constant $c^{\mathrm{kin}}_{\mu}$ that depends on $e$ but not on $\mu_{c}$), inconsistent with Eq. (1) and our dry and wet flow data (for dry flows, adding higher-order terms may remedy this discrepancy Berzi and Jenkins (2018)).

We define dilute and dense conditions – as opposed to rarefied ones – in terms of the mean free path $\lambda$ through the condition $\lambda<d$, where $\lambda(\phi)=\sqrt{2}d/(12\phi)$ for spherical particles and $\lambda(\phi)=\pi d/(12\sqrt{2}\phi)$ for spheres confined in two dimensions. In fact, we hypothesize that at least some of the few deviations from the scaling in Eq. (1) at large $\mathrm{Pe}$ [Fig. 2(a)] are related to a transition from dilute to rarefied flows at large shear rates, where $\mu$ is limited by the geometrical constraints of high energy collisions Pähtz and Durán (2018b).

*Variable yield stress ratio.—*Interestingly, deviations from the scaling in Eq. (1) at small $\mathrm{Pe}$ (larger-than-predicted values of $\mu$) are well characterized by the dimensionless granular temperature gradient $M\equiv d|\nabla T|/T$ and seem to occur whenever $M\lesssim 0.7$ and $\mathrm{Pe}\lesssim 0.5$ [Fig. 2(a), open symbols]. For uniformly sheared flows (squares in Fig. 2), where temperature gradients are negligible ($M\simeq 0$), these deviations are owed to the fact that $\mu$ converges to the yield stress ratio $\mu_{s}$ ($\mu^{\mathrm{2D}}_{s}=0.277$ Trulsson et al. (2012), $\mu^{\mathrm{3D}}_{s}=0.382$ Chialvo and Sundaresan (2013)) in the limit of vanishing shear rate. From Eq. (1), we find that this yield transition in homogeneous flows starts at $\mathrm{Pe}=(\mu_{s}/c_{\mu})^{2}\approx 0.5$.

For inhomogeneous flows, $\mu$ can be substantially smaller than $\mu_{s}$ when $\mathrm{Pe}\lesssim 0.5$ [Figs. 2(a) and 3(a)], at which point deviations from the scaling in Eq. (1) are controlled by the condition $M\lesssim 0.7$. These deviations have several elements in common with a yield transition, as illustrated in Fig. 3 for turbulent saltation transport, which is a nearly dry granular flow because of a large density ratio ($s=2000$) and large Stokes numbers [$\mathrm{St}\in(10,200)$]. First, $\mu$ seems to converge to a finite value $\mu_{s}^{\ast}$ in the limit of vanishing shear rate [Fig. 3(a)]. Second, the dimensionless characteristic length $M^{-1}$, associated with spatial changes in the granular temperature, collapses as a function of $\mu/\mu^{\ast}_{s}$, peaks with a finite value at $\mu/\mu^{\ast}_{s}=1$, and once $\mu/\mu^{\ast}_{s}\lesssim 1$, the data scatter [Fig. 3(b)]. This peak is similar to the divergence of the relaxation length associated with spatial changes of the shear rate ($\dot{\gamma}$) and granular stresses ($\mathbf{P}^{c}$) when approaching the yield condition ($\mu\rightarrow\mu_{s}$) in existing nonlocal rheology models Bouzid et al. (2015). Finally, $\phi$ approaches the packing fraction as $\mu\rightarrow\mu^{\ast}_{s}$ [Fig. 3(c)]. We thus conclude that $\mu^{\ast}_{s}$ is the analog of $\mu_{s}$ for inhomogeneous flows.

The onset of the yielding transition is thus controlled by the local values of either $\mathrm{Pe}$ for relatively uniform flows (i.e., relatively small $M$) or $M$ for relatively inhomogeneous flows (i.e., relatively large $M$). In the latter case, the yielding transition can expand over a range of $M$ [e.g., $0.37\lesssim M\lesssim 0.7$ for turbulent saltation transport, see Fig. 3(b)] and coincides with nonlocality in the relation between $\mu$ and $\mathrm{Pe}$. These behaviors are consistent with $\mathrm{Pe}$ playing the role of a granular fluidity. Fluidity inhomogeneities, which are associated with nonlocality Bouzid et al. (2015), then would decrease with increasing $M$ as shear rate inhomogeneities get compensated by temperature gradients, $(\nabla\mathrm{Pe})/\mathrm{Pe}=(\nabla\dot{\gamma})/\dot{\gamma}-(\nabla T)/(2T)$, rendering the rheology local. They are also consistent with vanishing normal stress differences (e.g., $P_{xx}/P_{zz}-1=0$) at $\mu=\mu^{(\ast)}_{s}$ in the large-system limit Clark et al. (2018) [Fig. 3(d)] because large velocity and temperature gradients generate normal stress differences Sela and Goldhirsch (1998); Saha and Alam (2014, 2016) and thus prevent the yielding transition. In particular, it seems that, for $M\gtrsim 0.7$, the generation of normal stress differences is sufficient to prevent the yielding transition, even for comparably small $\mathrm{Pe}$ and thus $\mu$ [Eq. (1)].

*Conclusions.—*In this study, we have shown that, under certain relatively weak constraints, the Mohr-Coulomb friction coefficient $\mu$ obeys the general scaling $\mu=c_{\mu}\sqrt{\mathrm{Pe}}$, with the Péclet number defined as $\mathrm{Pe}\equiv\dot{\gamma}d/\sqrt{T}$ (but, in general, disobeys scaling laws from viscoinertial rheology models and extended kinetic theory). This calls for the development of hydrodynamic models for dense granular flows involving granular temperature. Apart from extended kinetic theory, several such models have already been proposed to reproduce several aspects such as hysteresis Lee and Huang (2012); DeGiuli and Wyart (2017) and drag on an object Seguin et al. (2011). The scaling parameter $c_{\mu}$ varies with the tangential friction coefficient $\mu_{c}$ but not with the normal coefficient of restitution $e$, which led us to propose that $\mu=c_{\mu}\sqrt{\mathrm{Pe}}$ encodes the effects of the anisotropies of the particle assembly ($\mathrm{Pe}$) and tangential contact force network ($c_{\mu}$) on $\mu$.

The yield stress ratio of granular media, below which granular flows either stop or fundamentally change Bouzid et al. (2015), is often found to be independent of the flow geometry Clark et al. (2018). However, for wind-driven sediment transport, the scaling $\mu\propto\sqrt{\mathrm{Pe}}$ even holds for friction coefficients as low as $\mu\simeq 0.08$ and is accompanied by very low yield stress ratios (as low as $\mu^{\ast}_{s}\simeq 0.04$, reminiscent of vibrated granular flows Gaudel and De Richter (2019)), which seem to be caused by relatively large values of the dimensionless temperature gradient $M\equiv d|\nabla T|/T$. Future studies should investigate this link between $\mu^{\ast}_{s}$ and $T$ because it may play a role in explaining long-standing open problems, such as the reduction of friction in long-runout landslides Legros (2002); Lucas et al. (2014); Johnson et al. (2016).

I Supplementary Material

I.1 Nonlocality in sediment transport and gravity flows

When the dimensionless fluid shear stress (“Shields number”) $\Theta$ of the fluid driving sediment transport is too close to its value at the flow threshold or when the inclination angle $\alpha$ of the gravity flows in ambient air is too close to its value at the flow threshold, the simulation data do not obey the scaling $\mu\propto\sqrt{\mathrm{Pe}}$ (Fig. S1), which relates the local friction coefficient $\mu$ to the local Pèclet number $\mathrm{Pe}$ (i.e., the granular flow rheology is nonlocal). This finding is not surprising because nonlocal effects are already known to be crucial for dry gravity flows (e.g., they are responsible for the stopping angle dependency on the flow thickness Bouzid et al. (2015)) and because it is already known that sediment transport tends to creep below the surface of the granular bed Houssais et al. (2015), which is also associated with nonlocality Bouzid et al. (2015). In particular, for turbulent saltation transport, the granular bed does not flow liquidlike but creeps even relatively far from the flow threshold because the fluid shear stress at the bed surface is insufficient to mobilize particles Durán et al. (2011). Only for sufficiently intense conditions (i.e., those in Table S1), when collisions between particles of the rarefied transport layer and granular bed are so frequent that the bed no longer recovers between collisions, does the granular bed flow like a liquid and thus the rheology become local Pähtz and Durán (2018b).

I.2 Alternative manner to evaluate the Péclet number scaling

Figure S2 shows the scaling $\mu=c_{\mu}\sqrt{\mathrm{Pe}}$ for the same data as in Fig. 2(a) of the paper, but including the standard deviation for all our simulated flows (for uniform flows, it is usually smaller than the symbol size and therefore not shown).

I.3 Influence of contact parameters on scaling law

Figure S3 shows the influence of contact parameters on the scaling $\mu=c_{\mu}\sqrt{\mathrm{Pe}}$ for uniformly sheared particle and suspension flows. It can be seen that the normal restitution coefficient $e$ does not affect the scaling parameter $c_{\mu}$. In contrast, $c_{\mu}$ increases with the contact friction coefficient $\mu_{c}$.

I.4 Failure of standard rheology models

I.4.1 Failure of viscoplastic rheology and fluidity scaling

The viscoplastic $\mu(I)$ rheology predicts that there is a general relationship between the friction coefficient $\mu$ and the inertial number $I$ across dense, dry granular flows Jop et al. (2006). However, Figs. S4(a) and S5(a) show that this prediction fails for our simulated nearly dry flows. This failure cannot be remedied by nonlocal extensions of the $\mu(I)$ rheology. For example, it has been demonstrated that the nonlocal model of Ref. Kamrin and Koval (2012) is based on a general scaling of the rescaled fluidity $\dot{\gamma}d/(\mu\sqrt{T})$ with the particle volume fraction $\phi$ Zhang and Kamrin (2017), which is disobeyed across our dry flows [Fig. S5(b)].

I.4.2 Failure of viscoinertial rheology models

Figures S4(b) and S6(b) show that the failure of the $\mu(I)$ [Figs. S4(a), S5(a), and S6(a)] also cannot be remedied by the viscous number rheology $\mu(I_{\nu})$ Boyer et al. (2011). Interestingly, turbulent saltation transport roughly obeys the same power law with the viscoinertial number $\sqrt{I^{2}+5I_{\nu}}$ as uniformly sheared particle and suspension flows [$\mu\propto(I^{2}+5I_{\nu})^{0.15}$], but the proportionality constant is much smaller [Fig. S6(c)].

I.4.3 Failure of extended kinetic theory

Extended kinetic theory predicts that the particle shear stress $\lVert\tau^{c}\rVert\equiv|P^{c}_{1}-P^{c}_{2}|/2$ (for two-dimensional flows), rescaled by $\rho_{p}d\sqrt{T}\dot{\gamma}$, and the rescaled particle pressure $P^{c}/(\rho_{p}T)$ are functions of only the particle volume fraction $\phi$. However, it is well known that this prediction, in general, does not hold for wet granular flows Garzó et al. (2012); Chamorro et al. (2015); Saha and Alam (2017); Alam et al. (2019), which is shown in Fig. S7 for uniformly sheared particle and suspension flows.