Intense sediment transport: Collisional to turbulent suspension

Steady, uniform transport of sediments in shearing flows, of interest in many engineering applications and geophysical flows, depends on the characteristics of solid particles and fluid, gravity, strength of the shearing flow and boundary conditions. Even in the idealized condition of mono-sized particles transported by a turbulent fluid, a comprehensive theoretical approach based on mechanics and able to describe the flow regimes and their transitions is not yet available.

When the mean distance between the moving particles is so small that inter-particle collisions cannot be neglected, a continuum approach based on kinetic theory of granular gases^1,2 is appropriate to describe the stresses in the particle phase.³ We call collisional suspension the regime in which, on average, inter-particle collisions represent the only physical mechanism able to support the weight of the particles in the direction perpendicular to the flow. Recently, simplified approaches based on an approximate integration scheme of the balance equations and the constitutive relations of kinetic theory have been proposed to obtain full analytical solutions to steady, uniform, collisional suspensions driven either by pressure gradient⁴ or gravity.⁵ These analytical solutions have been successfully tested against laboratory experiments on transport of natural and artificial particles in water.

When the strength of the turbulent, shearing fluid increases, correlations in fluctuations of particle concentration and fluid velocity⁶ induce an additional force in the momentum balances—the turbulent lift—which provides a further mechanism to suspend the particles. A simplified, global approach to deal with the transition to this turbulent-collisional regime has been proposed.⁷ In the collisional and turbulent-collisional suspensions, the large-scale fluid turbulence is suppressed, as experimentally revealed.^8,9

A further increase in the strength of the shearing fluid would cause the weight of the particles to be entirely balanced by the turbulent lift, in a portion of the domain where the turbulence is more intense. This turbulent lift mechanism to suspend the particles is considered important when the ratio between the particle settling velocity and the fluid shear velocity is lower than or about one.^10,11

In turbulent suspensions, the seminal theory of Rouse¹² for the distribution of the concentration, based on a kinematic balance between sediment fluxes due to downward settling velocity and upward velocity of turbulent eddies, is generally accepted,⁸ although reasons to add further flux terms have been devised.^13,14 Attempts for an improvement of the original Rouse approach dealt with turbulence modulation induced by the presence of particles, represented as reduction of either the Karman’s parameter or the mixing length.^15–17 Little attention has been paid so far to the boundary conditions to assign at the bottom of the fully turbulent transport layer: when an erodible bed is present, heuristic assignment of concentration values at some distance from the bed itself has been proposed.^18,19

Some authors, in the past, tried to address the problem in a different way, based on mechanical equilibrium and continuous two-phase approach.^6,20–24 Ni et al.²⁵ proposed a solution of the Boltzmann equation in the dilute and dense cases, where the liquid effects on particle concentration distribution is roughly included through the particle settling velocity and the mixing length. Hsu et al.³ offered a continuous two-phase model, based on momentum and kinetic-energy conservations, able to include both turbulent and collisional suspensions. Recently, Chiodi et al.²⁶ developed a two-phase model where simplified, phenomenological constitutive relations are adopted for both the particles and the fluid. Apart from improvements necessary to model in a proper way the turbulent and collisional effects in regions of high particle concentration, all these approaches do not allow for an analytical solution to the steady, uniform turbulent suspension. Also, they do not deal with the relative importance of the different physical mechanisms acting on the flow.

Here, we extend our previous work on collisional suspensions^4,5 to include the effect of turbulent lift, and obtain analytical solutions to steady, uniform sediment transport over horizontal erodible beds driven by a pressure gradient. To do that, we use extended kinetic theory of granular gases^27,28 for the particles, and treat the flow as a succession of layers: one close to the bed—where we assume that the particles are at yield—in which the concentration exceeds 0.49 and the particle velocities are correlated;²⁹ one above it, in which the interparticle collisions and the turbulent lift cooperate in supporting the weight of the particles; one at the top in which the weight of the particles is balanced by the turbulent lift and collisions are absent. This description is an improvement on previous proposals where there was no distinction between correlated and uncorrelated particle velocities and the role of the turbulent lift was not taken into account.^30,31 We approximately integrate the governing equations in each layer and obtain distributions of particle velocity and concentration, from which global quantities can be easily extracted. We then compare our predictions with available laboratory experiments performed with glass spheres and sand particles in water, over a wide range of strength of the shearing fluid. Finally, we provide criteria to distinguish between the different regimes of pressure-driven transport of mono-sized particles over erodible beds and show how those criteria translate into a regime map in which the inputs are the intensity of the shearing fluid, and the properties of fluid and particles.

The sketch of the flow configuration is depicted in Fig. 1. We focus on the pressure-driven, horizontal flow of identical, inelastic spheres, of diameter d and density ρ_p, over an erodible bed immersed in a turbulent fluid. The fluid has density ρ and molecular viscosity η. We assume that both the particle and the liquid motion are steady and uniform. We take x and y to be the direction parallel and perpendicular to the horizontal bed, located at y = 0, and neglect the variation along the spanwise direction. We assume that at least part of the weight of the particles is supported by the turbulent lift; this is the case when the ratio of the settling velocity of a single particle to the fluid shear velocity at the bed is less than one^7,32—0.8 according to Sumer et al.¹¹ We take H to be the total height of the particles over the erodible bed; above H the turbulence is not enough to suspend the particles. We will explain later how to determine its value. The local particle velocity along x and concentration are u and ν, respectively; U is the fluid velocity in the x-direction. In what follows, all quantities are made dimensionless using the particle density and diameter, and the reduced gravitational acceleration, g(σ − 1)/σ, where g is the gravitational acceleration and σ the ratio of particle to liquid mass density. With this, the inverse of the dimensionless molecular viscosity of the liquid is the fall particle Reynolds number R.

The balance of particle momentum perpendicular to the flow is

where p is the particle pressure, C drag coefficient, S fluid shear stress, and l turbulent mixing length. Here, and in what follows, a prime indicates a derivative along y. The second term on the right hand side of Eq. (1) is the turbulent lift associated with correlated fluctuations in concentration and fluid velocity in the vertical direction as modelled by McTigue.⁶ We can relate the drag coefficient to the settling velocity W of many particles in a resting fluid, C = 1/W; where $W=w1−νn$⁠, with w single particle settling velocity, and n = ln(10⁻³Rσν_c/w)/ln(1 − ν_c).³³ The quantity ν_c is the volume fraction at which rate-independent components of the particle stresses start to develop.^34,35 We take this to be the concentration at the interface with the erodible bed. Then, Eq. (1) reads

We characterize the flow through the different layers depicted in Fig. 1. As already mentioned, we assume that the erodible bed is where the particles experience shear-independent component of the stresses, so that the concentration exceeds the value ν_c;³⁴ the interface between the erodible bed and the flow is where ν = ν_c, and there the ratio of the particle shear stress to the particle pressure has a characteristic yielding value α.³⁵ Close to the erodible bed, there is a layer comprising between y = 0 and y = Δ in which the concentration exceeds 0.49, the value above which the molecular chaos assumption of classic kinetic theory is not valid and extended kinetic theory applies.^{27–29,35,36} Close to the erodible bed there might be a sublayer where the fluid viscous forces become dominant and the mixture behaves as a viscous suspension.⁴ Here, we are not interested in distinguishing between the dense collisional and the dense macroviscous regime, as done in a previous work,⁴ and we generically refer to the region comprised between the bed and y = Δ as the dense layer. Moving away from that layer, between y = Δ and y = h, the concentration is less than 0.49, so that the particle velocities are uncorrelated and classical kinetic theories apply. Here, the turbulence increases, so that the weight of the particles is at least partially supported by the turbulent lift: the mixture behaves as a turbulent-collisional suspension.⁷ We treat the region between the bed and y = h as a boundary layer, i.e., we assume that there the dimensionless total shear stress of the mixture, indicated by θ, is constant, and therefore equal to the dimensionless shear stress exerted on the bed. The latter represents the Shields parameter of the flow. The sum of the fluid, S, and the particle shear stresses, s, are equal to the total shear stress; therefore, S + s = θ. Above y = h, and up to y = H, the turbulent lift balances exactly the weight of the particles. Then, the right hand side of Eq. (2) vanishes and so does the particle pressure p and the particle shear stress s; the inter-particle collisions cease and the mixture behaves as a fully turbulent suspension. We will see that, for strong enough shearing flows, the fully turbulent layer can extend down to y = Δ and beyond. In the fully turbulent layer, the shear stress decreases linearly with the distance from the bed.

The particle shear stress at y = Δ is a fraction of the total shear stress of the mixture there. As in Berzi,⁷ we assume that the ratio of the particle shear stress s to the total shear stress of the mixture at y = Δ is equal to the ratio of the particle pressure to the total submerged weight per unit area of the particles above Δ. The latter is given by the integration of Eq. (2). Making use of the trapezoidal rule to approximate the integral, and taking l = κy, with κ Karman’s constant, and the concentration linearly distributed in the dense layer, we obtain, from S_Δ = θ − s_Δ, where S is the fluid shear stress,

where ξ is a coefficient of order unity to be determined from fitting with experiments which contains all the neglected correlations. Hereafter, the index indicates the position y at which the quantity is evaluated. It is worthwhile mentioning that, in purely collisional suspensions, the turbulent mixing length is local, i.e., it is only a function of the local value of the mean interparticle distance (a fraction of the particle diameter).⁹ Here, given the influence of the turbulence on the particle motion, we take the turbulent mixing length to be nonlocal (proportional to the distance from the bed), as appropriated for turbulent fluids in absence of sediments. Further studies are necessary to understand the transition from local to nonlocal turbulence in suspensions. Equation (3) provides s_Δ as a function of the Shields parameter; s_Δ initially increases and then decreases with the Shields parameter. It vanishes for values of θ larger than $w2/ξ2κ2σ$⁠, when the shearing flow is sufficiently strong to entirely suspend the particles by means of the turbulent lift in the entire region above y = Δ. In this condition, the turbulent-collisional layer of Fig. 1 vanishes. This picture implies a situation in which the turbulence of the interstitial fluid is not suppressed at large concentrations. To our knowledge, this is a novel inference and there is currently no experimental validation or falsification of this result.

We assume that, in the dense suspension layer, the divergence of the fluctuation energy flux can be neglected in the balance of fluctuation energy for the particles.²⁸ In this case, the value, k, of the ratio of particle shear stress to particle pressure at y = Δ is^28,37

where $J=1+eΔ/2+π/43eΔ−11+eΔ2/24−1−eΔ11−eΔ$ and e is the coefficient of collisional restitution—the negative ratio of post to pre-collisional normal, relative velocity between two colliding particles—here taken to change with y. With k and the particle shear stress s_Δ, the particle pressure p_Δ at y = Δ can then be evaluated as p_Δ = s_Δ/k. The dependence of the coefficient of restitution on y is due to the lubrication forces that develop when two particle approach each other in presence of a viscous fluid, providing an additional damping to the collisions. This effect can be modeled by taking the restitution coefficient in collisions to depend on a local Stokes number,³⁸ St ≡ σT^1/2R/9, where T is the local granular temperature (one-third the mean square of the particle velocity fluctuations). Here, as in previous work,^4,5,7 we use the formula of Barnocky and Davis³⁹ 

where ε is the effective coefficient of collisional restitution of the particles in dry condition which also takes into account the role of friction in the energy dissipation.^40,41 The coefficient of restitution cannot be less than zero; hence, from Eq. (5) and the definition of the Stokes number, we calculate the coefficient of restitution at y = Δ as

The granular temperature at y = Δ can be calculated using the constitutive relation for the particle pressure of kinetic theory,² and the radial distribution function at contact of Torquato,⁴² 

Equations (3), (4), (6), and (7), along with p_Δ = s_Δ/k, allow to determine the values of k, e_Δ, T_Δ, and p_Δ through a simple iterative method. When e_Δ = 0, the iteration is not necessary, and we immediately obtain from Eqs. (4) and (7) that k = 0.82, so that p_Δ = 1.22s_Δ, and T_Δ = 2.44p_Δ. The coefficient of restitution at y = Δ is zero when $sΔ≤0.34T⌣$⁠, where, from Eq. (6), $T⌣=62.11+ε/εσR2$ is the minimum granular temperature for having e_Δ greater than zero. If the maximum particle shear stress $w2/4ξ2κ2σ$⁠, from Eq. (3), is less than $0.34T⌣$⁠, the restitution coefficient at y = Δ vanishes for every value of the Shields parameter.

If we approximate Eq. (2) in the dense suspension layer as

and we integrate using the trapezium rule, we obtain, with Eq. (3) and the relation between the particle shear stress and pressure at y = Δ,

where we have taken that the fluid shear stress is zero at the bed, so that p₀ = θ/α.

We adopt the mixing length approach to express the constitutive relation for the fluid shear stress as

Then, using Eq. (10) and the constitutive relation for the particle shear stress of kinetic theory, $s=8Jν2g0T1/2u′/5π1/2$⁠,³⁵ where g₀ is the radial distribution at contact,⁴² and assuming that s + S = θ and u′ ≃ U′ in the dense layer,

By taking the particle shear rate in the dense layer to be constant and equal to the average between y = Δ and y = 0, where it vanishes, we calculate the velocity at the top of the dense layer as

where we have taken, as already mentioned, l_Δ = κΔ. The latter assumption implies that large scale turbulent structures are present even at large concentrations. A direct experimental validation of this assumption is lacking; however, when the turbulence is strong enough to cause the vanishing of the turbulent-collisional layer (when p_Δ = 0 and S_Δ = θ), Eq. (12) indicates that the particle velocity at the top of the dense suspension layer is proportional to the square root of σθ, i.e., to the fluid shear velocity, in accordance with the experiments.¹¹ 

The particle flow rate per unit width at y = Δ is then

with the distributions of concentration and velocity taken to be linear in the dense layer, as shown in experiments,⁴³ when performing the integration.

The total shear stress θ is assumed constant here, like in the dense layer. Integrating Eq. (2) in the turbulent-collisional layer by using the trapezium rule, with the conditions that the particle pressure vanishes at y = h—the turbulent lift supports the weight of the granular material there—gives approximately

Equations (14) and (9) show that h is approximately linear in the Shields parameter as obtained in the experiments.⁴⁴ When s_Δ is null (see Eq. (3)), p_Δ is null too, h equals Δ (Eq. (14)), and the turbulent-collisional layer vanishes.

The concentration at the top of the turbulent-collisional layer can then be obtained by taking the right hand side of Eq. (2) zero there, with $ν′=νh−νΔ/h−Δ$ and s_h = 0, so that S_h = θ,

where we have ignored the concentration-dependence of the settling velocity at y = h, and we have taken l_h = κh.

The integration of Eq. (11) in the turbulent-collisional layer, using the trapezium rule, gives

We then obtain the particle flow rate per unit width at y = h by taking the concentration and the velocity linearly distributed in the turbulent-collisional layer, and integrating

In the fully turbulent layer, the weight of the particles is entirely supported by the turbulent lift. The shear stress decreases linearly with the distance from the bed. Previous studies^7,11,32 indicate that, in a global sense, there is a transition to a turbulent suspension regime when the ratio of the particle settling velocity to the fluid shear velocity at the bed is less than a certain value (0.8 according to Sumer et al.).¹¹ Here, we use this criterion as a local indication of where the particles cease to be suspended by the turbulence, by replacing the bed shear velocity, (σθ)^1/2, with the local shear velocity, (σS)^1/2, at position y. Hence, we determine the upper limit of the fully turbulent layer as the distance H from the bed at which $w/σSH1/2=0.8$⁠. If the fluid shear stress decreases linearly with the distance from the bed and vanishes at y = Y, we obtain

Actually, given that H cannot be less than h, we take the former to be the maximum between the value obtained from Eq. (18) and the value of h given by Eq. (14). It is worth emphasizing that the determination of Y in pressure-driven sediment transport over erodible beds is a hard problem. In the following, for simplicity, we will use the value of Y valid in absence of particles and postpone a more specific analysis to future work. In so doing, we are assuming that the increment of the cross section surface producing shear on the bed⁴⁵ is compensated by the reduction of area due to the sediment deposition.

If the weight of the particles is entirely supported by the turbulent lift, the particle pressure vanishes in the fully turbulent layer. By taking the right hand side of Eq. (2) zero when y exceeds h, we obtain

where, for simplicity, we have taken the fluid shear stress constant and equal to its average value $θ+w2/0.64σ/2$ between y = h and y = H. Integrating, we obtain the distribution of the concentration in the fully turbulent layer,

Although the way to get it and its boundary condition are different from previous works, Eq. (20) is similar to the well-known Rouse solution.¹² 

We assume that the velocity is constant in the fully turbulent layer and equal to u_h. This is corroborated by experimental results.^8,46 Assuming a logarithmic distribution would not substantially alter the results but would prevent us from obtaining an analytical expression for the total particle flow rate per unit width, q = q_H. The latter is

With the integration of Eq. (20), and the already mentioned assumption of linear distribution of the concentration in the dense and turbulent-collisional layers, the depth-averaged concentration in the flow reads

Figure 2 shows the particle flow rate q as a function of the Shields parameter as predicted by the present theory and measured in the experiments of Matoušek et al.⁴⁷ The experiments were performed with glass spheres of diameter equal to 0.18 mm and water (R = 5.8; σ = 2.45; w = 0.55) flowing in a rectangular pressurized conduit of height equal to 284 particle diameters. We take the shear stress to be zero at the centerline of the conduit, so that Y = 142 particle diameters. As appropriated for glass spheres, we use α = 0.38,⁴⁸ ε = 0.73 (corresponding to a normal coefficient of restitution of 0.9,⁴⁹ and a coefficient of surface friction equal to 0.5, which gives exactly 0.38 as yield stress ratio),³⁵ ν_c = 0.587,³⁴ and ξ = 0.01 and 0.2—at the moment the only adjustable parameter of the theory. If ξ is much larger than 0.2, the value of H of Eq. (18) becomes less than h; as already mentioned, if this is the case, we take H = h, so that q_H = q_h. The dashed line in Fig. 2 represents therefore also the prediction of the theory for large ξ. At the smallest values of θ, the theoretical flow rate q_h is closer to the experimental points. Apart from uncertainties in the experimental evaluation of the Shields parameter in pressurized conduits and in presence of erodible beds, this seems to indicate some overestimate of the concentration at y = h. When ξ = 0.2 is used, the present theory, based on mechanical arguments, satisfactorily reproduces the experimental data. A similar agreement can be obtained using the formula proposed by Matoušek,⁵⁰ which is however empirical and contains a number of coefficients and exponents that had to be calibrated by fitting with experimental data.

Figure 3 shows the depth-averaged particle concentration, ν_A, versus the Shields parameter as predicted by the present theory and measured in the experiments performed by Matoušek⁵⁰ with sand of diameter equal to 0.37 mm and water (R = 17.6; σ = 2.65; w = 1.09) flowing in a circular pressurized pipe of radius (and therefore Y) equal to 203 particle diameters. As for glass spheres, we use α = 0.38, ε = 0.73, ν_c = 0.587, and ξ = 0.2 in the theory. The theoretical curve presents a non-monotonic behavior, which is not apparent in the experiments, as more experimental points would be required at the lowest values of the Shields parameter. In absence of the fully turbulent layer—i.e., in collisional sediment transport⁵—and without changing the bed slope, the depth-averaged concentration decreases with the Shields parameter. When a fully turbulent layer is present, the concentration at y = h increases with the Shields parameter and tends to ν_Δ, so that also ν_A eventually increases with the Shields parameter. The non-monotonic behavior of the depth-averaged concentration is therefore a signature of the transition from collisional to turbulent suspension.

Figure 4 shows the prediction of the present theory in terms of particle velocity and concentration profiles at four different values of the Shields parameter for the flow of 0.37 mm sand in water. Also shown in Fig. 4(b) are the experimental results of Matoušek.⁵⁰ For consistency with our analysis, the experimental profiles of concentration have been shifted to have ν = ν_c at y = 0. The qualitative behavior of the experimental concentration profiles is well captured by the present theory. Even the quantitative agreement is satisfactory, considering that the values of ν_c and ν_Δ that we adopt are valid for spheres, not for irregular natural particles, and that we only roughly estimate the total depth H.

Figure 5 shows the theoretical values of the velocity at the top of the turbulent-collisional layer and at the top of the dense layer, along with the experimental measurements performed by Pugh and Wilson⁴⁴ with sand of diameter equal to 0.30 and 0.56 mm and water (R = 12.8 and 32.7; σ = 2.65; w = 0.94 and 1.35) flowing in a circular pressurized pipe— Y = 172 and 92 diameters, respectively. The experimental points correspond to the velocity at the top of the shear layer, defined by Pugh and Wilson as the layer in which the concentration is approximately linearly distributed. As for glass spheres, we use α = 0.38, ε = 0.73, ν_c = 0.587, and ξ = 0.2 in the theory. The velocity is divided by the bed shear velocity; in the experiments, this ratio is constant. Wilson⁵¹ showed that this must be the case, if the bed shear velocity is the only scale velocity in the problem. The present analysis suggests that there is a saturation process, as a consequence of the trade-off between the collisional and the turbulent components of the velocities: at large Shields parameters, the collisional components in Eqs. (12) and (16) vanish, and the fluid shear stress at y = Δ tends to θ. Hence, the proportionality of the velocities with (σθ)^1/2, indicating that the fluid turbulence dominates the particle motion. It is worthwhile emphasizing that this is the case only if the mixing length is taken to be nonlocal, so that the depths of the layers cancel out from Eqs. (12) and (16). The agreement with the experiments is good, considering that the position at which the experimental concentration profile ceases to be approximately linear should be located somewhere below the position h above which the concentration follows the Rouse-like profile of Eq. (20).

The present theory allows to define the boundaries between the different regimes of pressure-driven sediment transport over horizontal, erodible beds. In the dense and turbulent collisional layer, we made use of kinetic theory, which implies a continuous description of the flow. This assumption breaks down when the predicted height h is less than one particle diameter. As in our work on inclined sediment transport,⁵ we refer to that regime as ordinary bedload. The Shields parameter at which there is transition to ordinary bedload can be calculated from Eqs. (9) and (14), with h = 1 and taking s_Δ = θ (the turbulent lift can be ignored in the dense layer at small values of θ)

At Shields parameters larger than $θ⌣$⁠, collisions are able to continuously support the weight of the particles above the bed: we call this regime collisional suspension. The turbulent lift provides an additional mechanism to continuously support the weight of the particles above the bed when the ratio of the single particle settling velocity to the bed shear velocity is equal to 0.8,¹¹ i.e., at $θ =θ⌢=w2/0.64σ$⁠: at larger values of the Shields parameter, the mixture behaves as a turbulent-collisional suspension. Finally, when the particle pressure in Δ vanishes, i.e., as already mentioned, at $θ=θT=w2/ξ2κ2σ$⁠, the turbulent lift becomes the main mechanism to continuously suspend the particles everywhere in the flow, but in a dense sub-layer above the bed: there is a transition to a fully turbulent suspension.

The particle settling velocity w is a function of the density ratio and the fall particle Reynolds number.⁵² Hence, for a given type of particles (e.g., glass spheres), a given fluid (e.g., water) and a nearly horizontal, erodible bed, we can build a regime map which illustrates the range of existence of the different type of steady sediment transport in terms of the only two remaining free parameters of the problem: the Shields parameter and the fall particle Reynolds number. Such regime map is depicted in Fig. 6. Also shown is the critical Shields parameter that represents the threshold for having particle motion over erodible beds.⁵³ It is worth recalling that, at any given particle Reynolds number, a different regime map can be drawn in case of free surface, inclined sediment transport over erodible beds, in which the coordinates are the angle of inclination of the bed and the Shields parameter.⁵ In that case, another regime of sediment transport, i.e., debris flow, shows up.⁵ These regime maps allow for a quantitative classification of granular-fluid flows, therefore improving upon previous qualitative diagrams.⁵⁴ We postpone to future work the extension of the present analysis to inclined, sediment transport in presence of turbulent suspension, for which, however, detailed experimental measurements are still lacking.

In this paper, the steady and fully developed pressure-driven flow of sediments immersed in a turbulent fluid over an erodible bed has been modeled using kinetic theory of granular gases for the particle phase, and a turbulent mixing length approach for the fluid phase. We have distinguished between a dense layer close to the bed, in which the particle velocities are correlated; a turbulent-collisional layer, in which collisions and turbulent lift cooperate in suspending the particles; and a fully turbulent layer, in which the weight of the particles is entirely balanced by the turbulent lift. The approximate integration of the governing differential equations permits to obtain simple analytical expressions for the depths of the layers and the distributions of particle velocity and concentration, as functions of the Shields parameter. The agreement between the theory and laboratory experiments performed with either glass spheres or sand particles in water is notable, with few, measurable particle and fluid properties required as inputs to the model. Only one parameter of order unity necessitates to be tuned against experiments. The inclusion of the turbulent lift permits the theory to be valid for a range of Shields parameters enclosed between those of ordinary bedload—approximately four times the critical Shields parameter for the motion threshold—and infinity.

This research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915 to the Kavli Institute of Theoretical Physics.