The critical shear stress (τc) for grain entrainment is a poorly constrained control on bedload transport rates in rivers. Direct calculations of τc have been hindered by the inability to measure the geometry of in situ grains; i.e., the shape and location of each grain relative to surrounding grains and the bed surface. We present the first complete suite of three-dimensional (3-D) grain geometry parameters for 1055 water-worked grains, and use these to parameterize a new 3-D grain entrainment model and hence estimate τc. The 3-D data were collected using X-ray computed tomography scanning of sediment samples extracted from a prototype scale flume experiment. We find that (1) parameters including pivot angle and proportional grain exposure do not vary systematically with relative grain size; (2) τc is primarily controlled by grain protrusion, not pivot angle; and (3) larger grains experience larger forces as a result of projecting higher into the flow profile, producing equal mobility. We suggest that grain protrusion is a suitable proxy for assessing gravel-bed stability.


Predictions of bedload transport rates are important for predicting sediment sorting and channel evolution and for river management and restoration, but current approaches are accurate to an order of magnitude at best. Bedload transport models are typically produced from analysis of bulk average hydraulic and sediment transport data, and their poor performance may be improved by taking into account the grain-scale complexity of the process. Grain entrainment occurs when the local entraining forces of the fluid (applied shear stress, τ) exceed the resisting forces of the individual grain (critical shear stress, τc). Although the controls on τ are increasingly well understood (Bottacin-Busolin et al., 2008; Schmeeckle, 2014; Lamb et al., 2017), the controls on τc are still poorly constrained.

For each grain, τc is determined by its geometry within the bed; i.e., the grain’s size, shape, and position relative to surrounding grains and the velocity profile (Fig. 1). τc has been estimated from physically based grain-entrainment models that calculate the balance of forces or moments acting on a grain (e.g., White, 1940; Li and Komar, 1986; Wiberg and Smith, 1987; Kirchner et al., 1990; Vollmer and Kleinhans, 2007; Yager et al., 2018). These models use parameters to describe the grain geometry, including exposure (the area of the grain surface that is exposed to the flow), projection (the elevation of the grain with respect to the velocity profile), the pivot angle (ϕ) through which the grain rotates to exit its pocket, and the weight of any overlying grains. Another common parameter, protrusion, is the combination of exposure and projection. However, the difficulty of measuring the three-dimensional (3-D) geometry of sediment beds means that most models incorporate simplifications, e.g., spherical grains, assumed relationships between parameter values and relative grain size (grain diameter relative to median grain size, D/D50), movement in the downstream direction only, and 1-D simplification of 2-D properties like grain exposure.

Quantifying the different contributions to τc is necessary to develop predictive models, e.g., of the impact of cohesive sediment in a bed; but, they are difficult to disentangle when using traditional methods to collect field or laboratory data. Direct field-based measurements of the forces required to move an in situ grain (e.g., Johnston et al., 1998; Prancevic and Lamb, 2015) do not allow the individual parameter values to be resolved, because the measured force is the sum of the forces needed to remove the grain from its pocket, to displace any overlying grains, and to overcome any cohesive material around the grain (Jain and Kothyari, 2009; Barzilai et al., 2012; Hodge et al., 2013; Perret et al., 2018). Force measurements also do not account for hydraulic interactions with the bed; under the same applied shear stress, grains with different exposures would experience different forces depending on their exposed area and position within the velocity profile.

The limitations of models and field data mean that for grains within water-worked beds, we still cannot predict how grain geometry parameters vary with D/D50, nor which parameter exerts the most influence on τc. 3-D imaging techniques have been used to gain insights into granular processes such as infiltration (Kleinhans et al., 2008; Haynes et al., 2009). We have applied this approach to grain entrainment by using X-ray computed tomography (XCT) scans of water-worked sediment beds to produce the first complete suite of grain-geometry parameters from 1055 in situ grains. Measuring parameters directly removes the need for the previous model simplifications, and for each grain, we use our measurements to calculate τc in a vector-based 3-D moment-balance model for grain entrainment (Voepel et al., 2019). We demonstrate how grain geometry varies with grain size and controls τc, and compare the τc values to those derived from standard bedload-transport models.


The sediment assemblages used for this study were extracted from prototype scale flume (60 × 2.1 × 0.7 m) experiments that replicated the riffle-pool sequence morphology (Fig. DR1 in the GSA Data Repository1) and poorly sorted grain-size distribution (GSD) (D50 = 23 mm, D0 = 4 mm, D100 = 64 mm) of Bury Green Brook (Bishop’s Stortford, UK; presented in Hodge et al., 2013). In six flume runs, 0.25-m-diameter baskets (Figs. DR2A and DR2B) were buried in the riffle, deep pool, and pool tail, with their rims 1.5 D50 below the surface of the bed. In two runs, baskets were also buried in the shallow pool, giving a total of 20 baskets. Data from two pool-tail baskets were analyzed by Voepel et al. (2019). The baskets were covered with the bulk GSD, leaving surface grains free to move under flow. Six hours of flow at half the discharge at which D50 grains started to move created a water-worked bed that had a structure comparable to that of a natural river bed. Ockelford and Haynes (2012) showed that application of half the critical shear stress that mobilized D50 caused significant grain-scale rearrangement in gravel beds, and Powell et al. (2016) showed that the majority of bed restructuring occurred over flow durations of much less than 6 h. The baskets were then excavated, stabilized by dipping in molten wax, and XCT scanned. The bed was manually turned over prior to resetting the morphology and burying the baskets for the next run. Four flume runs (10 baskets) had additional fine sand and clay added to the bed.

The samples were scanned using a micro-focus Nikon Metrology μCT scanner. Following image reconstruction, the 3-D images, with a voxel size of 0.600 mm, were segmented into coarse grains (>4 mm), matrix (sand and clay), and background (air and wax) (Figs. DR2C–DR2E). All individual coarse grains were then isolated, and contact areas between them were identified (Fig. DR2F). After scanning, we measured the surface layer GSD of eight baskets. Here we analyze the combined data set from all 20 baskets, in order to identify trends that are common to water-worked sediment, rather than identify results that are specific to particular morphological units or the presence or absence of fine sediment. Before describing the grain-geometry parameters that we measured, we introduce the 3-D entrainment model.


The entrainment model (Voepel et al., 2019) resolves all moments for entrainment and resistance force vectors onto a grain’s axis of rotation, and identifies the point at which τ equals τc and, hence, the grain is on the threshold of motion (Fig. 1). In the model, the entrainment forces are drag and lift (FD and FL, respectively) and the resisting forces are the submerged grain weight and resistance caused by a cohesive matrix (FW and FC, respectively). The orientation of the rotation plane and vector moment arms are determined by the grain center of mass and the two contact points that form the axis about which the grain is rotating, and so, in contrast to previous models, the grain can in theory be entrained in any direction.

The applied τ determines the values of FD and FL, which are calculated using a 3-D version of the approach used by Kirchner et al. (1990), assuming a logarithmic velocity profile with the zero-velocity elevation equal to the mean bed elevation over one D84 upstream and downstream of the grain. For FD and FL, we use the actual projected grain areas perpendicular (A) and parallel (A||) to the flow, rather than a spherical approximation. FD is modified by an exposure factor (EFadj), which accounts for the sheltering effect of upstream grains and the possibility of flow deviation from the horizontal (Voepel et al., 2019). FC was not included in the Kirchner model, and is a linear function of the grain contact area with the cohesive matrix (FC = ηATfM, where η is cohesive force per unit area, AT is total grain surface area, and fM is the proportion of grain surface area in contact with matrix; Voepel et al., 2019). The model does not incorporate the resisting weight of any overlying grains, but we analyze only grains that can move without displacing adjacent grains. To apply the model to a grain, we measured grain volume, local bed elevation, center of mass, contact points with adjacent grains, and the previously defined parameters A, A||, EFadj, AT, and fM. We also calculated a 1-D projection value (maximum grain height above zero-velocity elevation, P), and the proportion of the grain area that is both above the zero-velocity elevation and not sheltered by upstream grains (Aprop = AEFadj). Aprop thus reflects protrusion, in that it incorporates both projection and exposure.

We calculated τc for all grains that were not fully below the zero-velocity elevation and that had at least one viable pair of contact points that allowed the grain to rotate without obstruction by another grain. Most measured grains had multiple viable pairs of contact points, so for each grain, we calculated τc for all viable pairs, and then chose the lowest τc. The contact geometry of this scenario determines the grain pivot angle (ϕ) and the entrainment bearing angle (β, where 0° is downstream and 90° is cross-stream); consequently, φ and β are outputs, rather than inputs as in previous models.


Data from the 1055 grains across six replicate flume runs (e.g., Figs. 2A–2F) show that there is not a strong relationship between any grain-geometry parameter and D/D50 (Figs. 2G–2L). (We use the bulk D50, but analysis using surface D50 for the eight-basket subset produces the same result.) The range of P/D decreases as D/D50 increases. However, Aprop, which combines the projected area of a grain and the effect of upstream sheltering, shows no relationship with D/D50. This is because for small grains with high P/D values, the effects of projection and upstream sheltering effectively cancel each other out. fM and φ have very weak negative and positive relationships with D/D50, respectively, and β has no relationship. Thirty-six percent (36%) of grains have φ > 90°, indicating grains with contact points that are above the center of mass of the grain. The wide range of β shows that it is not always easiest for grains to move in the downstream direction.

Critical shear stress (τc) varies weakly with D/D50. Instead, τc is most strongly a function of P/D and Aprop (Figs. 2N and 2M). For a given P/D, smaller grains tend to have a higher value of τc than larger grains (Fig. 2N). τc has a weak relationship with φ and with fM, and no relationship with β (Figs. 2O, 2P, and 2Q). The strongest identified relationship is between dimensionless critical shear stress (τc*) and P/D50 (Fig. 2R; R2 = 0.71). The data suggest equal mobility; values of the 5th percentile of τcc5), which we use to indicate when grains start to be entrained, are approximately constant across the GSD (Fig. 2L).


We compare our τc values to the Wilcock and Crowe (2003) and Meyer-Peter and Muller (1948) bedload transport models (Fig. 3). Using both our data and bedload models, at a range of dimensionless shear stresses (τ*), we estimate the number of grains entrained from an area equal to all twenty baskets over a suitable time step (see the Data Repository, section 3). The number of grains entrained at each τ* is similar for both the XCT data and the bedload models, indicating that our τc values are realistic. The two bedload models produce relationships with different shapes; our new approach can constrain this relationship and thus improve bedload transport models. At high τ*, we do not account for changes caused by entrainment of adjacent grains or ballistic entrainment by mobile grains, but the similarity between the bedload models and our data suggest that these are not substantial effects.


We find that D/D50 is a poor predictor for φ, β, and Aprop. This is consistent with force-gauge measurements from river beds (Johnston et al., 1998; Sanguinito and Johnson, 2012; Hodge et al., 2013; Prancevic and Lamb, 2015), suggesting that the bed structure we observed is common rather than specific to this GSD. Compared to our data (Fig. 2I), pivot angle measurements of grains placed on gravel surfaces (e.g., Kirchner et al., 1990; Buffington et al., 1992) have a smaller range, and also show a negative relationship between median φ and D/D50. These differences suggest that grains placed on the bed are not representative of natural river beds.

Our results show that P/D and Aprop are important controls on τc. Yager et al. (2018) used different theory and data and also identified the importance of protrusion, suggesting again that our findings are not specific to our particular GSD. Further analysis of our data, and data from other rivers, are required to test whether the specific form of relationships such as that between τc* and P/D50 (Fig. 2R) varies within and between natural river beds, and as a function of GSD, grain shape, and flow history. The importance of projection means that differences in grain mobility between different sediment beds could potentially be estimated from measuring the projection of sediment grains using high-resolution topographic data, and applying a relationship such as that derived from Figure 2R. Further work is required to establish whether differences in projection are adequately captured by bulk bed-roughness properties (as a proxy for projection), or whether the projection of individual grains would need to be measured (e.g., Masteller and Finnegan, 2017). Further analysis of grain projection could also have applications for flow resistance, because a measure of surface roughness has been found to be a more appropriate measure of flow resistance than grain size (e.g., Aberle and Smart, 2003).

Our finding of equal mobility (Fig. 2L) is consistent with previous analyses from Bury Green Brook and other low-gradient channels (Parker, 2008; Hodge et al., 2013), despite our analyses not considering how the bed structure might evolve during a transport event. Equal mobility has been attributed to smaller grains having higher pivot angles than larger grains (Kirchner et al., 1990), but our data do not support this. Instead, the area of grain that is exposed to the flow is a more important control on τc. For a given P/D, τc decreases with increasing D (Fig. 2N). This is because the increased mass is more than compensated for by the projection of larger grains higher into the logarithmic flow profile, such that for the same P/D, a larger grain will be exposed to higher velocities and therefore a larger force per unit exposed area.

Finally, our data show the potential impact of a cohesive matrix. Our median φ is 72°, which is less than the median φ of 85° back-calculated from force-gauge data from Bury Green Brook (Hodge et al., 2013). Additional components (including fM and overlying grains) are therefore contributing to the field force measurements, suggesting that bedload transport calculations should incorporate the effect of a cohesive matrix where present.


This work was funded by Natural Environment Research Council grants NE/K012304/1 and NE/K013386/1. The flume facility was provided by the University of Southampton Chilworth Hydraulics Facility, UK. XCT scanning was undertaken in the University of Southampton μ-VIS X-Ray Imaging Centre. Many thanks to Michael Church and anonymous reviewers for their helpful comments.

1GSA Data Repository item 2020044, details of the flume setup, example XCT data, and comparison with bedload transport models, is available online at http://www.geosociety.org/datarepository/2020/, or on request from editing@geosociety.org. Data in Figure 2 are available from http://doi.org/10.5281/zenodo.3478124.
Gold Open Access: This paper is published under the terms of the CC-BY license.