Development of a vector‐based 3D grain entrainment model with application to X‐ray computed tomography scanned riverbed sediment

Sediment transport equations typically produce transport rates that are biased by orders of magnitude. A causal component of this inaccuracy is the inability to represent complex grain‐scale interactions controlling entrainment. Grain‐scale incipient motion has long been modelled using geometric relationships based on simplified particle geometry and two‐dimensional (2D) force or moment balances. However, this approach neglects many complexities of real grains, including grain shape, cohesion and the angle of entrainment relative to flow direction. To better represent this complexity, we develop the first vector‐based, fully three‐dimensional (3D) grain rotation entrainment model that can be used to resolve any entrainment formulation in 3D, and which also includes the effect of matrix cohesion. To apply this model we use X‐ray computed tomography to quantify the 3D structure of water‐worked river grains. We compare our 3D model results with those derived from application of a 2D entrainment model. We find that the 2D approach produces estimates of dimensionless critical shear stress ( τcr∗ ) that are an order of magnitude lower than our 3D model. We demonstrate that it is more appropriate to use the c‐axis when calculating 2D projections, which increases values of τcr∗ to more closely match our 3D estimates. The 3D model reveals that the main controls on critical shear stress in our samples are projection of grains, cohesive effects from a fine‐grained matrix, and bearing angle for the plane of rotation (the lateral angle of departure from downstream flow that, in part, defines the grain's direction of pivot about an axis formed by two contact points in 3D). The structural precision of our 3D model demonstrates sources of geometric error inherent in 2D models. By improving flow properties to better replicate local hydraulics in our 3D model, entrainment modelling of scanned riverbed grains has the potential for benchmarking 2D model enhancements. © 2019 The Authors. Earth Surface Processes and Landforms Published by John Wiley & Sons Ltd.


Introduction and Problem Statement
The initiation of transport of grains in a stream is dependent upon the driving bed shear stress (τ b ) crossing some critical threshold for entrainment: the critical shear stress, τ cr , the dimensionless form of which is known as the Shields parameter, τ * cr (Shields, 1936). Although Shields (1936) defines τ * cr as the ratio of fluid force on the grain to the submerged weight of the grain, the comprehensive review by Buffington and Montgomery (1997) highlights that in natural systems there are many other factors, which affect entrainment. These factors include grain protrusion (defined fully below, but a function of projection from a mean bed elevation into the flow and grain exposure), intergranular geometry (a function of grain shape, sorting and packing), and fine sediment infiltration. These recognised characteristics of natural systems led to a number of definitions linking the threshold of motion to simple geometric properties and grain arrangements, whereupon entrainment models were developed using the literature linked to a combination of a lack of observational data and variation between derived models (Buffington and Montgomery, 1997).
Despite the acknowledged geometrical complexity of natural river systems, the approaches cited above make a number of significant simplifications in that they assume a given grain (i) is spherical or ellipsoidal (or indeed, circular or elliptical in 2D), (ii) pivots about a single contact point in the downstream direction, (iii) is unaffected by mortaring by fine sediment (Barzilai et al. 2013;Hodge et al. 2013), and (iv) is not partially buried beneath other grains. A few 3D entrainment models have been developed that are capable of more realistic geometric arrangements addressing the reliance on a single contact point noted above (Ali and Dey, 2016;Dey, 1999;Nabi et al. 2013); nevertheless, they still rely on spherical approximations of grains with relatively simple particle arrangements when compared to the natural conditions of riverbeds. Furthermore, the simplistic representation of grains as spheres can lead to further assumptions about simple ratios between grain diameter and projection being suitable to characterise the exposure of the grain. A number of studies have highlighted the fact that existing geometric grain entrainment models are developed for idealised laboratory or theoretical conditions and so do not reproduce the conditions found in natural riverbeds. Miller et al. (1977) showed that despite using a carefully selected set of entrainment data from non-cohesive, uniform and wellrounded or spherical grains, considerable scatter still existed in the threshold curves derived from their measurements. Similarly, Buffington et al. (1992) noted the relatively poor performance of the Wiberg and Smith (1987) model of grain entrainment, and advocated extending it to include the effects of packing, grain arrangement (including partial burial) and protrusion-induced flow acceleration as seen in natural systems.
Several studies address entrainment effects from grain exposure and projection through theoretical development (Kirchner et al. 1990;Wiberg and Smith, 1987) or flume experiments using natural river gravels (Fenton and Abbott, 1977), and from grain sheltering through remote sensing of bed topography (Measures and Tait, 2008). Measures and Tait (2008) identified two mechanisms of grain sheltering: direct sheltering, where adjacent particles prevent flow forces from directly acting upon the grain, and remote sheltering, where distant upstream particles modify local flow conditions in the vicinity of the grain. Empirical methods have been employed to identify the effects of friction angles and sorting (Johnston et al. 1998;Komar and Carling, 1991) as well as grain size and shape (Komar and Li, 1986) on sediment entrainment. And whilst this research identified important relationships between entrainment shear stress and granular properties of coarse river sediment and their distributions, these studies are still based on relatively simple geometries. Simplified models, such as those described above, do not take into account more complex grain characteristics such as particle arrangement and orientation as well as additional spatial angles that, along with pivot angle, better describe the full range of motion during grain entrainment.
To improve the performance of geometric grain entrainment models there is a need to better account for the complexities of natural systems, including addressing the fundamental point that riverbed structures are three dimensional in nature. Existing models, such as those identified above, treat the entrainment and resistance forces used in 2D terms and/or use simple grain geometry and arrangements, giving rise to an oversimplification of the problem. Similarly, current entrainment models do not account for the cohesion forces that fine sediment can exert on grains in natural systems (Barzilai et al. 2013;Hodge et al. 2013). Until recently, 2D models have been a necessary simplification because of the difficulty in measuring the 3D geometry of a water-worked bed of grains. However, techniques such as high-throughput X-ray computed tomography (XCT) provide a new opportunity for measuring 3D bed structure, providing a motivation for the development of new 3D models.
Application of any 3D entrainment model requires high-resolution data of a grain arrangement from which to calculate grain parameters including size and geometry, spatial orientation, exposed area, elevation and the proximity of grain-tograin contact points. Non-intrusive imaging methods such as high-throughput XCT have been used in geosciences and geomechanics for a number of years (Wildenschild and Sheppard, 2013;Ahmed et al. 2016b;Callow et al. 2018), for example to explore pore networks (Pierret et al. 2002), the breakthrough of solutes in porous media (Clausnitzer and Hopmans, 2000) or hyporheic zone sediment structure (Chen et al., 2009(Chen et al., , 2010Packman et al. 2006). Despite various hydrological applications (Wildenschild et al. 2002), and some preliminary explorations using magnetic resonance imaging (Haynes et al. 2009;Kleinhans et al. 2008), there are no published examples of using XCT to quantify river bed sediment geometry in a way pertinent to grain-scale entrainment models.
This paper contributes two advances. The first is a novel 3D entrainment model that employs vector mechanics to account for the 3D geometric structure of grains and their surroundings as well as incorporating the effects of cohesive resisting forces. We chose to use the 2D Kirchner et al.(1990) geometric grain entrainment model as a basis of comparison with our 3D entrainment model. The second advance is developing a method for collecting and processing high-resolution XCT data from sediment beds. We use the new 3D model and XCT data to establish the τ * cr for individual grains within two river bed samples taken from a prototype scale flume experiment: one of coarse grains and one with additional cohesive sediments. We compare our results to τ * cr values calculated using the Kirchner et al. (1990) model, and also evaluate the impact of the new cohesion term.

Grain Rotation Threshold of Entrainment Models
Here vector algebra is applied to a 3D moment balance to derive a physically-based threshold of motion entrainment model, which greatly simplifies the 3D modelling process when compared with the oft-used trigonometric method. Our focus is exclusively on grain rotation as it is the primary mechanism of entrainment (Agudo et al. 2014), although it should be recognised that all grains can undergo a mixture of entrainment modes throughout the sediment transport process (Ancey et al. 2006). Since our 3D model implements XCT scanned images of bed sediment, it utilises spatial proximity of surrounding grains and their contact points. Therefore, we choose an existing 2D grain protrusion-based entrainment model to facilitate comparison with our 3D vector-based entrainment model. model is presented, although it is noted that any other model based on grain protrusion could be used to compare with our vector-based 3D entrainment model developed below. The key components of the 2D entrainment model are defined below, but for full details see Kirchner et al. (1990). Kirchner et al. (1990) used the 2D framework derived by Wiberg and Smith (1987), where they account for bed slope in a static force balance used to derive their entrainment model. Kirchner et al. (1990) presumed negligible bed slope for their model by setting bed slope to zero. Here we use a trigonometric method to derive the same 2D entrainment model that Kirchner et al. (1990) used, except based on a static moment balance of forces about a point of contact between two grains, which forms a rigid body centre of rotation in the entrainment model (head of vector C in Figure 1a). The driving forces of entrainment, namely drag, F D , and lift, F L , act upon a circular particle and are functions of a logarithmic velocity flow profile, u(z), (Equation (1a)). These forces oppose a resistant force, namely submerged grain weight, F W . A moment of force is the perpendicular distance from the centre of rotation for a rigid body to the line of action of the force times the magnitude of that force. Summing all of the moments (positive for clockwise rotation; negative for anticlockwise rotation) and setting their sum equal to zero yields the static moment balance whereupon entrainment occurs (Equation (1b)). The entrainment forces of drag and lift depend on grain protrusion, which is defined by two components: projection (p) of the top of the grain above mean bed elevation; and exposure (e) of the top of the grain above local upstream maximum bed elevation (Kirchner et al., 1990, Appendix). It is theoretically assumed for the latter that the grain of interest, or test grain of diameter, D, sits on a bed of uniform grains of median grain size, K 50 , and is a function of the mean friction angle, α, such that their geometric relationship is linked (Figure 1b, Equation (1c)). Thus the three key equations that define the 2D Kirchner model are where u(z) is flow velocity relative to height, z, above mean bed height at z=0, ρ is the density of water, τ b is the boundary shear stress, κ=0.407 is van Karman's constant, f(z) is a logarithmic profile, and z 0 = K 84 /10 is a length scale. The drag force, F D , is expressed as an integral of the product of dynamic pressure, q ¼ 1 2 ρuðzÞ 2 , and the cross-sectional area of the grain, dA=w(z)dz, integrated over the protruding portion of the grain where C D is the drag coefficient, and w(z) is the width of the grain cross-section as a function of elevation z. The lift force, F L , is the product of the dynamic pressure differential across the top and bottom of the grain and the planar cross-sectional area of the grain where C L is a lift coefficient and A is the plan view crosssectional area of the grain. In Equations (2) and (3), Kirchner et al. (1990) used coefficient values of 0.40 and 0.20 for C D and C L , respectively. The submerged weight force, F W , is defined as where ρ s and ρ are the density of sediment and water, respectively, and g is gravitational acceleration of 9.81ms −2 . Substituting Equations (1a), (1c) and (2)-(4) into the 2D moment balance, Equation (1b), Kirchner et al. (1990) expressed the τ cr as Developing a vector-based 3D grain rotation model To account for the geometrical complexity of grains found in natural river systems, we use a vector-based 3D moment balance approach for resolving τ cr . The approach essentially adapts the components of fluid mechanics (Pritchard and Mitchell, 2015) used in the 2D critical shear model (Kirchner et al. 1990) for use within a 3D rigid body mechanics framework (Beer et al. 2013) and also accounts for an additional resisting force associated with cohesion. It is important to point out that the model being developed is fairly simple, assuming a logarithmic velocity flow profile where all forces are configured in a simple orientation relative to the system reference frame. Consequently, any additional driving and resisting forces or fluvial mechanics modifications Figure 1. (a) Diagram of forces acting on a grain: F L , F D and F W are the lift, drag and submerged weight forces, respectively, acting on the grain. U is the mean flow velocity at height Z above the mean bed elevation. (b) Idealized geometry for calculating grain projection (p) and exposure (e) as a function of the diameters of the test grain (D) and the bed particles (K) and the mean friction angle represented by the pivot angle (α). Redrawn after Kirchner et al. (1990).
3D ENTRAINMENT MODEL WITH APPLICATION TO CT SCANNED RIVERBED SEDIMENT could subsequently be included as components within the 3D vector-based framework (e.g. forces related to flow turbulence).
A 3D vector is written as F ¼ F x î þ F y ĵ þ F zk , where F x , F y , and F z are the scalar components, and î; ĵ, andk are the orthonormal unit vectors corresponding to the (x,y,z) axes of a Cartesian reference frame. To simplify the 3D entrainment model, the reference frame is chosen such that force vectors are oriented along single vector components. Here the stream flow velocity is defined in the direction of the x-axis, uî (ms −1 ), and gravity in the direction of the z-axis, gk (9.81ms −2 ). Unless otherwise stated, all vectors in this section are given as 3 × 1 column vectors to ensure correctly applied vector calculations. Typically, in the 2D case, grains are assumed to rotate in the direction of flow, and a pivot angle is formed between the gravity vector and the vector extending from the grain's centre of mass to a single point of contact with an adjacent downstream grain. However, in reality, grains do not pivot about a single contact point but rather they do so about an axis of rotation (AOR), A, formed by the vector difference of a left contact vector, CV L , and a right contact vector, CV R , between two contact points of adjacent grains ( Figure 2a). Each contact vector extends from the grain's centre of mass to a contact point with the adjacent grain. For illustration, think of a grain rotating within the saddle formed by two downstream grains. All of the action of a pivoting grain occurs within the plane of rotation (POR), which is orthogonal to the AOR (x ′ z ′ -plane in Figure 2a, b). The 3D orientation of the POR is induced by the rotation axis, where all of the 3D vector mechanics affecting rotation about an axis essentially collapse onto the 2D subspace we refer to as the POR.
Within any 3D axis rotational framework, an arbitrary vector, V, in 3D vector space, may be resolved into two orthogonal component vectors: one vector parallel to the AOR, V A , and another vector that spans the POR, V P (note vector subscripts in Figure 2a,b). For example, the drag force, F D , shown in Figure 2a, is parallel to the downstream flow direction along the x-axis. Here drag force is resolved into two orthogonal component vectors: one that is projected parallel to the AOR, F DA , the other that is projected onto the POR, F DP . Of these two orthogonal component vectors, only the POR force projection affects entrainment of the grain since the AOR force projection is parallel to the rotation axis. As an analogy, an open door does not swing shut due to gravity, which is a vector parallel to the door's AOR along its hinges; you must apply an external, nonparallel force to push the door shut. Figure 2 shows that all forces involved in grain entrainment act along multiple angles, only one of which is the pivot angle. The overhead view of the grain (Figure 2a) shows that the arrangement of the two contact points and the grain's centre of mass determines a bearing angle (β, Figure 2a) and a tilt angle (γ, not shown) for the POR as well as a pivot angle (α, Figure 2b) that lies within the POR. We define the bearing angle as the angle formed between the POR and the direction of stream flow, and the tilt angle as the angle formed between the AOR vector and the horizontal xy-plane. These two departure angles, one from the downstream flow direction (bearing) and the other from gravity (tilt), define the 3D spatial orientation of the POR as it begins its forward rotation (pivot) at the threshold of entrainment.
To obtain the pivot angle, either contact vector plus the gravity vector must be projected onto the POR before calculating the angle between these projected vectors. The 3 × 3 matrix that projects any vector onto the AOR is λ(λ T λ) −1 λ T , where λ is the unit vector for the AOR ( Figure 2a) and superscript T is the transpose of a vector. The projection that is orthogonal to the AOR projection, which projects any vector onto the POR, is I 3 −λðλ T λÞ −1 λ T , where I 3 is the 3 × 3 identity matrix. Hence, without loss of generality, the projection of CV L onto the POR is given by CV LP ¼ ðI 3 −λðλ T λÞ −1 λ T ÞCV L (Figure 2a). Similarly, the bearing angle is found by projecting the x-axis unit vector î onto an 'un-tilted' vertical POR. Let λ XY be λ, where thek scalar component is set to zero and normalized. Then I 3 −λ T XY ðλ T XY λ XY Þ −1 λ XY is the matrix that projects any vector onto the vertical POR, and Figure 2. (a) Overhead image of a 3D ellipsoidal grain on a sediment bed is overlaid by a vector diagram shown in the xy-plane, with z-axis going into the page. The AOR vector (A, blue arrow) formed by the vector difference of left and right contact vectors (CV L and CV R , red arrows) that extend from the centre of mass of the grain to where the grain makes contact with two neighbouring surface grains (green stones). Vector A forms a pivot for grain rotation when acted upon by a drag force (F D , violet arrow). The drag force resolves as two component vectors: one that is parallel to the AOR (F DA ) and one that resides within the POR (F DP ). Vectors that reside within the POR, denoted by the P subscript, contribute to pivoting the grain (e.g. F DP and CV LP ). (b) The same grain viewed within the 3D POR reference frame (x ′ ,z ′ ) indicates the y ′ -axis is coming out of the page. The 3D POR view shows component vectors of resultant forces (F EP and F RP , violet arrows) and their positions (R EP and R RP , red arrows), which are the two moment arms for applied forces acting about the rigid body centre of rotation (A, circular blue arrow). Angles for pivot, α, and bearing, β, are indicated. Tilt angle, γ, which is the departure angle of the rotation axis from the horizontal plane, is zero in the diagram and not shown. [Colour figure can be viewed at wileyonlinelibrary.com] H. VOEPEL ET AL.
î XY ¼ ðI 3 −λ T XY ðλ T XY λ XY Þ −1 λ XY Þî is the projection of the unit vector î onto the vertical POR. Finally, we use the definition of the dot product to find the angles for pivot, α, bearing, β, and tilt, γ, which are given by where G P is the projection of gravity gk onto the POR determined by the same projection used for contact vector projections, and |A| is the norm of arbitrary vector A. It is important to realize that the angles defined by Equations (6) are unnecessary for a vector-based 3D entrainment model. Nevertheless, since they are important to describe entrainment behaviour and to evaluate the 3D model, these angles are calculated once threshold of entrainment is determined. Our 3D entrainment model uses four forces acting upon a grain: the drag force, F D , the lift force, F L , the submerged weight of the grain, F W , and a cohesive force, F C , due to fine-grained matrix making contact with the particle (Figure 2b). The resultant forces for entrainment F E =F D +F L and for resistance F R =F W +F C are each applied at positions R E and R R , respectively. Each position vector extends from anywhere along the AOR (either contact point is most convenient) to the line of action of the applied force (Figure 2b). The static moment balance about an AOR is given as the scalar triple product The dot product in Equation (7) is the scalar resolute of the resultant moment vector, ∑RÂF, in the direction of the AOR unit vector, λ. Recall that vector R×F is orthogonal to the 2D plane spanned by vectors R and F. Thus the scalar resolute in Equation (7) is equivalent to collapsing all of the R and F vectors onto the POR as vector components (Figure 2b), calculating the resultant of all their moments, and then measuring its magnitude, yielding a signed scalar value (see Appendix A for a geometric proof). Subsequent subsections give details for the entrainment and resistance force vectors and their corresponding position vectors given in Equation (7).

Entrainment force components for the 3D model
We assume non-slip conditions at the bed and a logarithmic velocity profile with increasing elevation, z, up the water column. Flow velocity is zero at the mean bed elevation. Both drag and lift forces are functions of the velocity profile, which is given by where κ=0.407 is van Karman's constant, τ b is the boundary shear stress in Nm −2 , ρ is water density assumed to be 1000 kgm −3 , z 0 = D 84 / 10 is a bed roughness length scale, reference elevation, z=0, starts at the local mean bed elevation, which is taken over a distance of D 84 upstream and downstream of the grain, and z∈Z are discrete non-negative elevations above the mean bed. Drag force in vector format is given by where C D is the drag coefficient (assumed to be 0.91 for natural sediment after Schmeeckle et al. (2007) and will be used in both 2D and 3D models for model comparison), u(z) is the logarithmic velocity profile defined in Equation (8), A ⊥ (z) is the cross-sectional area of the grain perpendicular to flow as a function of discrete elevation height increments z∈Z above the mean bed elevation, F u is the fraction of grain cross-sectional area A ⊥ , where u>0, EF is an exposure factor to account for the sheltering effects of upstream grains, and EF adj is the necessary adjustment for the portion of EF already accounted for by F u .
To estimate exposure factor, EF, for any given grain we update the concept of calculating sphere protrusion areas (Yager et al. 2007) to include a weighted average of multiple area ratios over several viewing angles of the bed. To better account for the turbulent nature of streamflow, we assume that flow velocity (an instantaneous vector we term V) influencing the drag force could interact with the grain at any angle between horizontal and vertical relative to the bed, but that only the horizontal component of V is relevant to the resulting drag force (V sin θ, where θ is the angle from the overhead view).The exposure factor, EF, is the weighted average of exposure area ratios estimated by where ER θ is the exposure area ratio at viewing angle θ, which is defined as the exposed area of the grain viewed from angle θ divided by the total area of the grain viewed from the same angle but with no other grains present, and Θ is the collection of equidistant angles from 0°at the overhead view of the bed to 90°at the downstream horizontal view of the bed (see Appendix B for physical derivation of EF and EF adj ). There are four noteworthy observations about Equation (10): (i) normalizing the exposed-view area by total-view area yields a dimensionless metric; (ii) the metric is physically consistent with the product of dynamic pressure and the cross-sectional grain area, 1 2 ρu 2 A, used in drag force calculations; (iii) we do not need to know velocities or grain areas, as they cancel across numerator and denominator; and (iv) the weights sum to unity as required for any weighted mean. Lift force in vector format is given by where C L is the lift coefficient assumed to be 0.20 for natural grains (Schmeeckle et al. 2007), A ‖ is the overhead cross-sectional area of the grain parallel to flow, and u(z T ) and u(z B ) are the logarithmic velocity profile values at the top and bottom of the grain, respectively.
The position vector for entrainment, R E ¼ R Ex î þ R Ey ĵ þ R Ezk , extends from the AOR to the intersection of perpendicular lines of action for lift and drag forces (see Figure 2b). The line of action for the lift force runs vertically through the grain centre of mass. Due to the location and distribution of the velocity profile, however, the line of action for the drag force is horizontally situated somewhere above the grain's centre of mass. To find the elevation along the velocity profile for the horizontal line of action for drag force, we need to find the midpoint 3D ENTRAINMENT MODEL WITH APPLICATION TO CT SCANNED RIVERBED SEDIMENT of distributed drag force across the face of the grain. By using the drag force distribution as a weighting mechanism, the elevation for the horizontal line of action for drag force is found by calculating where z E is the elevation above the mean bed located along the velocity profile.
Resistant force components for the 3D entrainment model There are two resisting forces: a submerged grain weight and a cohesive force due to contact with any fine-grain matrix. Submerged weight force in vector format is given by where ρ s −ρ is the submerged density of sediment grain assumed to be 1650kgm −3 , g is the gravitational constant of 9.81ms −2 , and V S is the grain volume (m 3 ). Cohesive force in vector format is given by where η is the cohesive force per unit area (N mm −2 ), A S is the surface area of the grain (mm 2 ), and f S is the fraction of the grain surface area in contact with matrix fines. An empirical power law model was derived from a cohesive tensile force experiment where a 25mm diameter marble was pulled with a handheld force gauge from several binary mixtures of sand and clay (p <.001, R2=.54 (log); see Appendix C for details). This model is given by where η ϕ is the force per area power law (N mm −2 ) for grains with median axes within matrix sand size range ϕ, D ϕ is the mean for matrix sand grain size range ϕ, F is the fraction of clay in the sand-clay matrix, and the log-linearized model regression coefficients are a=0.002579, b=−0.5332, and c=0.2328. Sieving analysis rarely yields a uniform distribution of mass across all sand size classifications. As such, we use Equation (14b) to estimate a mass-weighted mean of cohesive force across grain sizes. Metric grain size ranges for sand classifications are Φ={(0.0625,0.125],(0.125,0.25],(0.25,0.5],(0.5,1], (1,2]} mm, and the set of means for each grain size range is D Φ ={0.0938,0.188,0.375,0.75,1.5} mm. The mass-weighted mean force per area η used in Equation (14a) is given by =where η ϕ is the force per area from Equation (14b) and m ϕ is the sediment mass for size range ϕ∈Φ determined from sieving analysis of the sediment used in XCT scanning. Equation (14c) is necessary since the cohesive tensile force experiments were performed on binary mixtures of a narrow range of grain sizes for sand mixed with clay whereas sand in the cohesive matrix from the flume experiments consists of a broad range of grain sizes for sand mixed with clay. Each force contribution from the narrow grain size range, η ϕ , is weighted by the mass of grains, m ϕ , that were sieved within that size range.
The position vector for resistance, R R ¼ R Rx î þ R Ry ĵ þ R Rzk , extends from either contact point on AOR to the centre of mass of the grain (see Figure 2b). Observe that the projection of R R onto the POR is just −CV LP (or the equivalent vector −CV RP ). Hence, either contact point works.
Calculation of τ cr for the 3D entrainment model For any given grain, the centre of mass and coordinates of all nearby contact points are used to construct all vectors and metrics. A set of critical shear stress values, {τ cr }, is calculated for each stone using every viable pair of contact vectors, (CV L , CV R ). Each viable pair of contact vectors induces a unique set of two position vectors, R E and R R , used in the entrainment model (Equation (7)), where their positions extend from either contact point to the grain's centre of mass for the resistance vector and to the horizontal line of action (Equation (12)) for the entrainment vector. For each grain, entrainment computation starts by calculating the resistance force vector, F R =F W +F C , for the entrainment model. A list of all viable pairs of contact vectors, (CV L ,CV R ), and their corresponding position vectors is then constructed. The viability of each pair of contact vectors is established using four tests: (i) the POR passing through the grain's centre of mass must intersect the AOR between the pair of contacts; (ii) the POR must not tilt so severely that both contact vectors are on the same side of a POR-aligned vertical plane containing the gravity vector; (iii) the POR bearing angle is bounded by |β|< 90°; and (iv) any contact pair must be eliminated from the list if a third contact point is positioned such that it prevents the grain from forward rotation. From all viable contact pairs, the smallest critical shear value in the set, τ cr =min{τ cr }, and its corresponding contact pair, (CV L ,CV R ), are denoted as the threshold of motion solution for that grain.
The most efficient way to calculate the scalar triple product, λ·R×F, is to express it as a 3 × 3 matrix with the AOR unit vector, λ ¼ λ x î þ λ y ĵ þ λ zk in the first row, the position vector, R ¼ R x î þ R y ĵ þ R zk , in the second row, and the resultant force vector, F ¼ F x î þ F y ĵ þ F zk in the third row, and then evaluate its determinant Recognizing that we need one determinant matrix for each position vector in the equation for a static moment balance about an AOR, λ·∑RÂF ¼ 0, our 3D moment balance (Equation (7)) yields the following scalar equation Algorithmic solutions, using iterations of boundary shear stress (see Appendix D), should be calculated using the scalar triple product form, Equation (7) or (15a), until the scalar moment reaches zero whereas analytical solutions should use the scalar equation form given by Equation (15b). Factoring by force components with all resistant force terms placed on the right-hand side of Equation (15b) and all entrainment force terms on the left-hand side, we solve for the boundary shear stress yielding the critical shear equation For clarity, we partition the critical shear stress formula (Equation (15c)) into two components: the first term is the 3D equivalent to the 2D entrainment model of Kirchner et al. (1990) and the second term is the cohesive force component due to particle contact with the fine-grain matrix. The Λs are the linear combinations of scalar components, in units of length (m), consisting of products of components from the AOR unit vector and position vectors; they represent the 3D geometry of grain-to-grain rotation for the rigid body mechanics of the entrainment model. Scalar force component functions, F D and F L , are attributed to the continuum mechanics for drag and lift forces, respectively, and are shown here as functions of cross-sectional grain areas (m 2 ) perpendicular and parallel to streamflow, respectively. Terms in the denominator are in units of cubic length (m 3 ) whereas those in the numerator are in units of force-length (N m), making the τ cr in units of force per area (N m −2 ). Dimensionless τ cr is calculated as τ * cr ¼ τ cr =½ðρ s −ρÞgD for grain diameter D.

Application of the Model to XCT Scanned River Grains
We used XCT to image an extracted section of an analogue riverbed, which was water-worked in a prototype scale flume. The 3D image was then processed to derive grain metrics, which were used to apply the new 3D entrainment model.

Development and extraction of analogue riverbed sediments
A prototype scale riffle-pool sequence was constructed in a large (60m long ×2.1m wide×0.7m deep) flume at the University of Southampton's Chilworth Hydraulics facility (Figure 3a), using river gravels with a grain size distribution that was carefully matched through sieving to those found in a natural riffle-pool system at Bury Green Brook, UK (mean D 50 range is 25-34 mm; Hodge et al., 2013). Prior to water-working, a series of five 250 mm diameter metal mesh baskets were buried with the rim 1.5×D 50 below the surface of the bed, allowing the surface grains of the bed to move without being impeded. Waterworking of the gravels was achieved by slowly raising the flow rate at the start of each flume run to a critical flow level (when several painted D 50 stones placed on the riffle mobilized), where it was sustained for 15 minutes prior to reducing the flow F F F to a half critical level for 6 hours. Powell et al. (2016) show that the majority of bed adjustment occurs over this timescale. The mean ± standard deviation of half critical flow pooled over all flume runs was 6.26×10 2 ±5.97×10 −3 m 3 s −1 , with half critical flow for individual flume runs ranging from 5.44×10 2 ±1.87×10 2 m 3 s −1 to 7.03×10 2 ±2.02×10 2 m 3 s −1 . For this work, we examined two end-member extremes in grain size distribution from the pool tail: (R1) coarse gravels (minimum grain size =4mm) with no fines and (R2) the same coarse sediment with the additional input of sands and clay (1-4 μm).
The clay was continuously added to the inlet flow to achieve a steady concentration of 500 mg L −1 throughout the run, measured using a conductivity meter with a predetermined calibration. In both runs the final sediment was framework supported. At the end of each run, the baskets were carefully extracted and briefly submerged in liquid wax to fix the in situ spatial arrangement of bed grains prior to transport to the scanner (Figure 3b).

XCT imaging, reconstruction and registration
The two baskets of sediment were imaged using a micro-focus Nikon Metrology μCT scanner at the μ-VIS X-Ray Imaging Centre, University of Southampton, UK. The image scanner has a 450 kVp X-ray source and a 2048-pixel curved linear detector array (CLDA). Data were acquired using an electron accelerating potential of 440 kVp and a tube current of 922 μA. 951 equiangular projections were acquired through 360°with 4 frames per projection taken to reduce noise. An exposure time of 60 ms was used. A detailed discussion of XCT artefacts and corrections is presented in Clausnitzer and Hopmans (2000). 3D reconstruction of the projections was performed using NI-TRO, a quasi-Newton and trust region differentiable optimization algorithm found in Digi XCT (Digisens, 2014;More and Sorensen, 1983;Nocedal and Wright, 1999). The resulting isometric voxel (the 3D extension of a 2D pixel) size is 600μm. Interactive volume registration was performed using VGStudio Max 2.1 (Volume Graphics, 2011). Three reference points were used to rotate and align the 3D datasets orthogonal to the original directions of flow and gravity, which is an important stage for later calculations of grain geometry and exposure area in relation to downstream flow direction ( Figure 3c).

Image processing and model application
After reconstruction and registration, the XCT data consists of a 3D array of voxels whose greyscale intensity is representative of the material density at that point. Image segmentation was undertaken to split the image into coarse grains (larger than 4 mm; Figure 3d, in red), and finer grains of sands and clay (termed 'matrix' and applicable only in R2) along with air/wax ( Figure 3d, in green), and was performed in Fiji ImageJ (Schindelin et al. 2012). Segmentation is performed using a semi-automated classification process in which the user sequentially defines known regions of air and wax (background), grains and matrix to train the classification algorithm, resulting in binary masks of each type of material (Figure 3e,f). Despite the larger sand grains in the matrix being greater than the image resolution, it was not possible to individually identify them as it takes 2-4 voxel widths at a 600 μm resolution to describe any grain with even rudimentary detail; consequently they were included in the general matrix class. Following classification, an initial separation algorithm was applied to the grain dataset to isolate individual grains. An adverse effect of image segmentation and separation is the inherent surface erosion over all of the separated particles in the image. The extent of erosion is determined by a number of factors including, but not limited to, scan quality, material density differences, segmentation steps and the separation algorithm. No contact actually occurs between separated particles in the images. Two sets of 2D image slices, each set composes a 3D image stack: a binary image stack of coarse grains ( Figure 3e) and a greyscale image stack of fine-grain matrix (similar to Figure 3 f), were imported into MATLAB (MathWorks, 2018). Four categories of grain metrics were extracted prior to entrainment calculations: (i) grain characteristics consisting of unique grain identification, centroid, volume, surface area, maximum grain elevation and its coordinates; (ii) principal axis lengths and their spatial orientation; (iii) grain-to-grain relationships consisting of coordinates of contact points with surrounding grains, their minimum separation distances and identification numbers of all contact grains; and (iv) matrix image processing to remove air pores and to calculate the fractional areas of matrix fines in contact with coarse sample grains. Imaging software packages typically have functions to extract characteristics such as labelling particles and obtaining their centroids, areas and bounding boxes. Herein we used MATLAB with the Image Processing Toolbox and the Parallel Processing Toolbox to compute metrics. Processed image data and developed code are available online at https://github.com/NERCPATCheS/ VectorEntrainment3D.
XCT scan data imported into MATLAB do not intrinsically have any information about scaling required to make real spatial measurements. To obtain spatial dimensions from array calculations for distance, area or volume, the calculated results must be multiplied by the voxel resolution, squared resolution, or cubed resolution, respectively. Dividing any of these converted spatial dimension by the appropriate voxel resolution converts the spatial dimension for each voxel back into an index value; this facilitates coding of entrainment model formulas where indices are required (Equations (8), (9a), (11) and (12)). For our MATLAB coding, we left coordinates, such as centroids, as index values to avoid rounding errors in moment calculations using Equation (7), where we used a binary search algorithm rather than the analytical solution (Equation (15c)) to calculate our results (see Appendix D). Note that the Λs are dimensionless when Equation (15c) is used with indices rather than distances to find analytical solutions of grain entrainment for XCT scanned images.
Individual grain characteristics must first be established before particle group relationships can be derived. The surface area of each grain was calculated by eroding the 3D binary image of grains with a 3D ball structural element to create labelled perimeters for each grain of one voxel thickness (Russ, 2011). Tabulating the surface and internal labelled voxels for each grain allowed the surface area and volume to be readily calculated. The maximum elevation of each labelled grain was also recorded. Axis length and orientation of image particles are usually calculated using ellipsoidal fitting (Ahmed et al. 2016a;Ge et al. 2005), but this method performs poorly when grains are more angular, as they often are in natural systems. Instead, we performed a principal component analysis (PCA) on all voxel coordinates of individual grains to produce axis lengths and grain orientation metrics. Contact points were found using an efficient search algorithm that calculates the minimum sums of squared distances between all surface voxels of an object grain and nearby surface voxels of adjacent grains. A threshold distance is then used during entrainment calculations to identify contact points in close enough proximity for consideration as a potential contact stones. Entrainment results were sensitive to small increments in threshold values as the number of stones having viable entrainment calculations increased from very few to nearly all surface stones. An example of contact points is shown in a transparent mock-up image (Figure 3g), which illustrates a group of grains (grey) shown with partial contact vectors (green) that H. VOEPEL ET AL. extend from the centre of mass of an object particle (centre grain) to the contact points (red patches) of its neighbouring grains. Greyscale image thresholding was used on the matrix images to segment pore space (i.e. darker greyscale values) from the fine-grain matrix, which converts the remaining matrix into a binary image. The labelled image of surface grains was then used with the binary matrix image to determine the amount of surface area in contact with matrix voxels, and a contact area fraction, f S , was calculated for each grain.
Once grain metrics were computed and tabulated, we used the scalar triple product form of our 3D entrainment model (Equation (7)) along with a binary search algorithm (see Appendix D) to calculate τ cr for each grain. The corresponding pivot angle and two POR angles for bearing and tilt were then calculated (Equation (6)). Finally, τ * cr was calculated from the iterated solution of τ cr for each grain that experienced sufficient shear stress for entrainment. The 3D model was run both with and without the cohesion term on the R2 sample to demonstrate how the cohesion term affects estimates of τ * cr when fines are present.
To provide a comparison of our vector-based 3D entrainment model with the 2D model of Kirchner et al.(1990), 2D metrics were also computed for all surface grains, denoted as any grain visible when looking directly down onto the R1 sample from overhead. Vertical 2D image slices were selected for each surface grain whereby each slice was parallel to downstream flow and contained the (x,y,z) coordinates for the 3D grain's centre of mass (see Figure 3e). Each slice was used to obtain the (x,z) coordinates for the contact point and the grain tops for each object grain and its upstream neighbouring grain. The 2D moment arm vector R ¼ Δxî þ Δzk was calculated, where Δx and Δz are the distances from the contact point to the grain's centre of mass. Grain diameters (b-axis and c-axis) were determined by the 3D PCA routine outlined above. D 50 and D 84 were obtained from sieving analysis of the R1 sample, and are used to calculate the empirical protrusion parameters and the local bed area, respectively. Empirical protrusion parameters were calculated by direct measurements from the vertical 2D image slices: projection, p, was defined as the vertical distance from the top of the object grain to the mean bed elevation that was calculated for the 3D results; and exposure, e, was defined as the vertical distance from the top of the object grain to the top of the upstream neighbouring grain. Theoretical protrusion parameters for projection, p, and exposure, e, were calculated for 2D entrainment based on both b-and c-axis lengths (Equation (1c)). The local bed area was defined as the D 84 distance upstream and downstream of the 3D bounding box for the grain, the latter of which was found using image processing functions. The mean bed elevation is thus the average elevation of all surface elevations within this bounded area.

Results
Our new vector-based 3D grain entrainment model and the protrusion-based 2D grain entrainment model of Kirchner et al. (1990) were used to calculate pivot angle, α, projection, p, and dimensionless critical shear stress, τ * cr , for 28 viable surface grains from the R1 sample, allowing comparisons to be drawn between 2D and 3D models. To ensure that 2D versus 3D treatment of geometry and protrusion were compared rather than the input data alone, we ran the 2D model based on empirical protrusion in addition to those based on theoretical protrusion (Equation (1c)). Empirical 2D parameters were measured p and e values taken directly from processed XCT images. We further assessed 2D-to-3D ratios of model parameters for α, p and τ * cr where the 2D models were based on b-axis, c-axis and empirical measurements for p and e. Calculations of τ * cr , both with and without the cohesive force given by the second term in Equation (15c), were made for 40 grains in the R2 sample, allowing direct assessment of cohesive effects. Specific grains from the R1 sample were selected to highlight entrainment metric differences under various complicated grain arrangements and orientations. We used both samples to examine whether water-worked grains lay flat. Finally, we explored the relationship between τ * cr and p both without and with cohesive force using R1 and R2, respectively.
Comparison of the 2D and 3D grain entrainment models R1 sample 2D entrainment metrics calculated using theoretical protrusion parameters based on b-and c-axis grain lengths, 2D entrainment metrics using empirical protrusion measurements, and our 3D entrainment metrics are shown in Table I.
Distributions for α, p and τ * cr were lognormal for all 2D and 3D entrainment metrics. Statistics associated with lognormal distributions are the geometric mean,μ G , and standard deviation,σ G , which are given in Figure 4 (Limpert et al. 2001;McAlister, 1879). Quantile intervals (QIs) for lognormal values,μ G ⋇σ G ¼ ½μ G Âσ G ;μ G ÷σ G , are the geometric counterpart toμ±σ used for normal values (Limpert et al. 2001) and are given for comparison of dispersion aboutμ G . A random variable X is said to have a lognormal distribution if logðXÞ has a normal distribution, denoted by logðXÞ∼Nðμ; σ 2 Þ, with mean μ and variance σ 2 (Casella and Berger, 2001). Since the 2D-to-3D ratios were also lognormal, we can test equality of our 3D metrics with each set of 2D metrics through 2D/3D=1, and take its log-transform, logð2DÞ−logð3DÞ ¼ 0 . Since logð2DÞ−logð3DÞ∼Nðμ d ; σ 2 d Þ with mean μ d ¼ E½logð2DÞ−logð3DÞ and variance σ 2 d ¼ var½logð2DÞ−logð3DÞ, we can use a paired sample t-test to check differences between various 2D models and our 3D model by H 0 :μ d =0 versus H A :μ d ≠0. A paired difference is a blocking method that removes correlated variation due to particle-to-particle variability (e.g. size, shape) allowing only model differences to be tested (Ott and Longnecker, 2001).
Here we examine parameter values for α, p and τ * cr (Figure 4 a-c) using entrainment calculations from the 2D b-axis theoretical method (2Db, green), the 2D empirical method (2De, orange), and our 3D vector entrainment method (3D, violet). Comparison p-value results from paired sample t-tests of logtransformed 2D-to-3D parameter ratios for α, p, and τ * cr (Figure 4d-f) are from the 2D b-axis theoretical method (green diamonds), the 2D c-axis theoretical method (fuchsia squares), and the 2D empirical method (orange circles), each compared with our 3D vector method. Four R1 sample grains were excluded from Figure 4 due to infinite 2D τ * cr values, which occur whenever pivot angles are high enough to generate negative projections (using Equation (1c)).
The pivot angles, α, have aggregate similarity as inferred by a geometric mean ofμ G ¼ 1:07 for the α 2D -to-α 3D ratio (Figure 4  d). There is moderate scatter around the geometric mean,σ G ¼ 2:02, because whereas the 2D model assumes that the grain pivots around a single contact point in the downstream direction, the 3D model allows for pairs of contact points in a wide  (Figure 4a). This implies that 3D entrained grains can pivot about a wide range of pivot angles within the POR, as the POR itself diverges from flow direction (the bearing angle) and gravity (the tilt angle). For any given grain, there is no systematic pattern as to whether an α 2D or α 3D value is higher.
Projections calculated for 2Db are greater than those for 3D, with a geometric mean for p 2D -to-p 3D ratios ofμ G ¼ 2:15 and only five grains having p 2D /p 3D <1. This difference is amplified in the calculation of τ * cr , with a relatively high 2Db projection leading to relatively lower τ * cr values; hence τ * 2D -to-τ * 3D ratios have a geometric mean ofμ G ¼ 0:32 with considerable scatter relative to the mean,σ G ¼ 5:05 (Figure 4f). The 2Db and 3D QIs for p are [8.5mm,36.8 mm] and [3.2mm,21.0mm], respectively (Figure 4b), differing by +163% and +75% at the lower and upper bounds, respectively. This discrepancy at the endpoints is reflected in the corresponding 2Db and 3D QIs for τ * cr at [0.033,0.750] and [0.105,2.328], respectively, differing by an average of about −112% at the endpoints.
The difference in 2Db and 3D projection values reflects the fact that the 2Db projection is calculated under the assumption that the grain is sitting on a uniform bed of D 50 grains (Equation (1); Kirchner et al., 1990), and does not include the increased sheltering effect of larger grainssomething that is explicitly accounted for with our 3D exposure factor. In the 2Db model, the entrained grain is assumed to be circular, with a diameter equal to the b-axis of the original grain. We also calculated the 2D p and τ * cr values using the c-axis, 2Dc, which reflects the fact that water-worked grains are typically arranged with their c-axis vertical (Komar and Li, 1986), and thus c-axis may be more relevant for calculating grain protrusion used to obtain entrainment forces of drag and lift. Using the c-axis reduces the p 2D -to-p 3D ratio mean fromμ G ¼ 2:15mm toμ G ¼ 1:35mm (Figure 4e), and almost triples the τ * 2D -to-τ * 3D ratio mean fromμ G ¼ 0:32 toμ G ¼ 0:86 (Figure 4f). Model 2Dc had a geometric mean and standard deviation ofμ G ¼ 11:1 mm andσ G ¼ 2:44mm for p, andμ G ¼ 0:425 toσ G ¼ 4:54 for τ * cr (not shown in Figure 4). The 2Dc model QIs for p and τ * cr are [4.6mm,27.1mm] and [0.094,1.932], respectively. The p endpoint discrepancy between 2D and 3D QIs decreases from +163% and +75% for 2Db values to +41% and +29% for 2Dc values at the lower and upper bounds, respectively. Similarly the τ * cr endpoint discrepancy between 2D and 3D QIs increases from −214% and −210% for 2Db values to −12% and −21% for 2Dc values. This reduction in p and τ * cr QI endpoint discrepancy from 2Db to 2Dc models, each relative to the vector-based 3D model, suggests 2D calculations using the c-axis may be more appropriate than using the b-axis. Note: Values in bold font indicate grains illustrated in Figure 6. Blank grain nos. are not listed in Figure 4 due to infinite 2D τ * cr values. Grain sizes were determined by the PCA method. Metrics labelled '2Db', '2Dc' and '2De' are calculated from b-axis, c-axis and empirically derived lengths (from XCT images), respectively. Grain no. is used in Figure 4; grain ID is used in Figure 6.
As an attempt to place 2D calculations on an equal footing with our 3D entrainment calculations, we used direct measurements of projection, p, and exposure, e, to derive an empirical entrainment model, 2De. p 2D -to-p 3D ratios for 2De are 27% smaller than those for 3D ( Figure 4e) with geometric metrics ofμ G ¼ 0:79mm andσ G ¼ 1:40mm, respectively. This small difference nearly doubles τ * cr values for 2De relative to our 3D model (Figure 4f) with geometric metrics for τ * 2D -to-τ * 3D ratios ofμ G ¼ 1:95 andσ G ¼ 2:38, respectively. The 2De model QIs for p and τ * cr are [2.3mm,18.5mm] and [0.187,4.982], respectively (Figure 4bc). The p endpoint discrepancy between 2D and 3D QIs increased from −40% and −13% for 2De values to +41% and +29% for 2Dc values at the lower and upper bounds, respectively. Similarly, the τ * cr discrepancy between 2D and 3D QIs decreases from +78% and +114% for 2De values to −12% and −21% for 2Dc values. The 2D-to-3D discrepancy shift from 2De to 2Dc is a much larger swing than that from 2Db to 2Dc, the former having a respective total percent change of 81% and 42% for lower and upper bound of p and 90% and 135% for that of τ * cr . Paired sample t-test of log-transformed 2D-to-3D parameter ratios were performed to evaluate mean ratio differences from unity of various 2D methods and our 3D method (Figure 4df). There were insignificant differences (at the 5% level) in α ratios across all methods. p 2D -to-p 3D and τ * 2D -to-τ * 3D ratios for 2Dc are insignificantly different from unity, whereas those for 2Db and 2De have significant differences. We infer from paired ratio t-test that 2D estimates of τ * cr using the c-axis (p=.623) is a more suitable proxy for 3D grain entrainment calculations compared to the commonly used b-axis (p=.001) or using empirical measurements (p<.001).

Effect of the cohesion term on critical shear stress
We use the data from the R2 sample, which contained additional fine sediment, to investigate the effect of our cohesion term in the 3D entrainment model (Equation (15c)), by calculating τ * cr with and without the cohesion term included. On average, the cohesion term approximately doubles the aggregate τ * cr ( Figure 5), but there is a considerable spread in the data. Eight out of 40 surface grains are unaffected by cohesion, whereas for one grain (stone no. 35) τ * cr increases nearly tenfold when adding the cohesion term. This variation is because different grains in the sample have different matrix contact areas, with some grains on the surface of our sample having no matrix contact, whereas those that are buried deeper may have more stone surface area contact with matrix fines. The scatter plot of τ * Fc ¼0 =τ * Fc ≠0 ratios suggests this to be the case as the grain numbers shown are ordered from highest to lowest surface elevations (Figure 5b).

Analysis of specific exemplar grains from the XCT image
To demonstrate the veracity of the new 3D model we look in detail at six exemplar grains from the R1 sample. Table I lists (d-f) Parameter ratios 2D-to-3D of individual surface grains using the 2D Kirchner model for theoretical b-axis (green diamonds), c-axis (fuchsia squares) and empirical p and e values (orange circles). Ratios are shown for: (d) pivot angle, α 2D -to-α 3D , (e) projection, p 2D -to-p 3D , and (f) dimensionless critical shear stress, τ * 2D -to-τ * 3D . Plots (d-f) show statistics for ratios (geometric mean and standard deviation, and p-value for paired ratio t-test) along with vertical lines for geometric means (coloured dot-dash) relative to unity (solid black). Solid symbols highlight grain IDs 3, 6, 7, 12, 14 and 19, which are analysed further below and in Figure 6. Four surface grains were excluded as obtuse 2D pivot angles generates infinite τ * cr values. [Colour figure can be viewed at wileyonlinelibrary. com] 3D ENTRAINMENT MODEL WITH APPLICATION TO CT SCANNED RIVERBED SEDIMENT 2D and 3D entrainment calculations for these exemplar stones (shown in bold font). Figure 6 shows the 3D representation of these grains from the XCT image, each grouped with their left and right contact grains, identified by their Grain ID (see Table I). Qualitative characteristics of each exemplar stone, relative to theoretical protrusion values based on b-axis calculations, are given in Table II and briefly discussed below.
Grains 3 and 7 have 3D τ * cr values that are considerably less than the 2D equivalents, but for different reasons. Grain 3 shares a downstream contact point between the 2D and 3D models, but the additional contact point for the 3D model (the right contact) lies below the grain, resulting in a low pivot angle of 11.6°in the 3D model relative to high pivot angle of 82.2°in the 2D model. Alternatively, the lower 3D value of τ * cr for grain 7 is attributed to the larger projection of 15.0 mm for 3D compared to 9.14 mm for 2D, as the pivot angles used for the 2D and 3D models are similar. Grain 12 requires a higher τ * cr using the 3D model compared to the 2D implementation (over 38 times higher than 2D). The combined effect of much lower projection (∼74% relative difference) and much larger pivot angle (about 3 times) for the 3D model relative to the 2D model contributes to the large τ * cr difference. The two downstream grains visibly making contact with grain 12 are not accounted for in the 2D model since neither is included in the 2D vertical plane passing through the grain 12 centre of mass. Grains 6, 14 and 19 have similar τ * cr values for the 2D and 3D entrainment models; however, we note that this is not because the 2D model accurately accounts for entrainment geometry; rather, interacting metrics counterbalance one another to coincidentally arrive at similar estimates. For example, grain 14 has a wide angle between contact points, reducing the distance Critical shear, τ * No X X X Note: Median angle values for R1 are used to partition high and low bearing and tilt angle values. Comparisons are based on b-axis calculations (2D-b in Table I).
between the centre of mass and the AOR in the 3D model relative to the 2D, but this reduction in length of both position vectors is offset by a moderately high bearing angle. Grains 6 and 19 have lower pivot angles in the 3D model implementation, but these are offset by relatively high bearing angles. These exemplar grains illustrate how 3D grain orientation combined with contact grain arrangement affects adjustments in τ * cr values in ways that are far too complex for 2D modelling.
Do water-worked grains lay flat?
There has been some research into whether grains are preferentially deposited with either the b-or c-axis aligned with gravity (Komar and Li, 1986). Our XCT-derived grain geometries from the R1 and R2 samples suggest that there is considerable variation in our water-worked grains (Figure 7) but a general preference for them to be positioned with their c-axis more closely aligned to vertical than either of the other two principal axes. The box plot in Figure 7a indicates a c-axis alignment with gravity is within about 45°for about half of all grains, whilst the a-axis and b-axis fall within about 65°for about half the grains. Similarly, results in Figure 7b show the c-axis alignment with the downstream direction are up to as much as 75°for about half of the grains. This result further supports our earlier conjecture that c-axis lengths would be a better proxy upon which to base 2D grain entrainment calculations than b-axis lengths.
Grain projection relationship to τ * cr Grain protrusion (and likely projection, p) is known to have a fairly strong correlation with τ * cr (Fenton and Abbott, 1977). Data from both the Coarse R1 and Coarse + Fines R2 samples show a strong to moderate power law relationship between τ * cr and p with 85% and 54% explained variation in logðτ * cr Þ, respectively ( Figure 8a). Dummy variable regression of linearized power law models combines R1 and R2 data for testing linear model differences (Draper and Smith, 1998). To perform this test, let logðτ * cr Þ ¼ α þ βlogðpÞ þ ϵ be the linearized power law for (i) Coarse and (ii) Coarse + Fines models. Let dummy variable Z=1 for the Coarse model and let Z=0 otherwise. Then the combined model that allows us to simultaneously test for 3D ENTRAINMENT MODEL WITH APPLICATION TO CT SCANNED RIVERBED SEDIMENT equal slope, α≡logðα ′ Þ, and scaling, β, parameters is logðτ * cr Þ ¼ α 1 þ α 2 Z þ β 1 logðpÞ þ β 2 Z logðpÞ þ ϵ , or, written in parsed form logðτ * cr Þ ¼ ðα 1 þ α 2 Þ þ ðβ 1 þ β 2 ÞlogðpÞ þ ϵ Coarse; where the error term assumption is ϵ∼N(0,σ 2 ). Using a partial F-test to test the hypothesis H 0 :α 2 =0andβ 2 =0 versus H A :α 2 ≠0orβ 2 ≠0, model parameters are significantly different (p<.001). The power law models suggest τ * cr is higher for grains in contact with matrix fines, and yet τ * cr differences between models decreases with increased projection. τ * cr for Coarse + Fines scales faster (β=−1.78) with increasing projection than Coarse grains alone (β=−1.52). This is likely the result of less surface contact with matrix fines as the grains project further from the river bed. Cohesion comparison ( Figure 5) suggests this might be the case as the lowest grain numbers in the plot are stones positioned at the highest elevations. Plotting the combined R1 and R2 samples (Figure 8b) suggests matrix fines might play a dominant role as the scaling parameter (β=−1.75) is similar to R2.

Sources of geometric error inherent in 2D models
The wide-ranging differences in entrainment parameter values between the 2D and 3D models are primarily attributed to a number of geometric factors: (i) the 2D model is a special case of the 3D model, with a range of motion restricted to a 2D flowparallel vertical plane, and where the AOR coincides with a single contact point; (ii) in the 3D model, the AOR (and therefore the POR) orientation is induced by variations in contact vector pair arrangements affecting entrainment results; (iii) variations of position vector lengths induced by contact pair arrangements affect the magnitude of the position vector and its contribution to the moment balance; and (iv) changes in force vector magnitude and direction induced by variations in POR orientation also affect moment balance. These factors affect how forces are geometrically transferred through the POR into the rotational mechanics of entrainment, and each is discussed below by considering a 2D entrainment model with fixed force and position vectors in static equilibrium about a single contact point (Figure 1a). We first show that the 2D model is a special case of the 3D model, case (i) above, and then we proceed to use it in a geometric argument to show the remaining cases. When combined, these factors show the geometric biases inherent in a 2D model. Given the careful consideration in our attempt to match the continuum mechanics for 2D and 3D entrainment models, we suggest these 2D geometric biases as the primary difference between 2D and 3D entrainment models.
For this discussion, consider a vector-based 2D entrainment model with fixed vectors in static equilibrium about a single contact point (Figure 1a). Consider the xz-plane as the POR for the 2D model (shown in Figure 1a). Now, by the right-hand rule for reference frames, the AOR has unit vector λ ¼ − ĵ. Reversing the direction of the single contact vector, C, in Figure 1 a, the position vector for the 2D model is R ¼ − 1 2 D sin αî þ 1 2 D cos αk. The force vectors for drag, lift and submerged grain Using the matrix determinant method to calculate the scalar triple product, we arrive at the same 2D moment balance as Equation (1b), thereby showing that the 2D entrainment model is just a special case of the 3D model restricted to a 2D subspace in the form of a flow-parallel vertical plane. Now the POR of our fixed 2D model contains all of the force vectors involved in the 2D moment balance. By parsing 2D force vectors into those with horizontal lines of action (drag force) and those with vertical lines of action (lift and submerged weight), we can examine the effects of bearing and tilt angles separately. By increasing the bearing angle, β, from the 2D position, thus reflecting a changing pair of contact vectors in 3D, the increased POR deflection from downstream flow (x ′ -axis in Figure 2a) has no effect on the contribution from vertical forces since those vectors remain within the POR. However, the drag force is effectively reduced by a factor of cosβ as the POR swings away from the horizontal line of action (vector F DP in Figure 2a). For the moment to remain in balance, the boundary shear stress must be increased to compensate for the loss in driving force being projected onto the POR, which is caused by an increased bearing angle. Similarly, increasing tilt angle, γ, from the 2D position for the POR has no effect on contribution of the horizontal force on entrainment, but the tilted POR does affect how vertical force contributes to overall entrainment, as vertical forces are effectively reduced by a factor of cosγ. Again, this requires an adjustment in the boundary shear stress for the moment to remain in balance. However, since the resultant of vertical forces, ðF L −F W Þk, is largely cancelled out by their difference, any increase in bed shear stress contributes more towards the magnitude of drag force in the total entrainment, F E =F D +F L , which is also partly attributed to the fact that drag and lift forces scale proportionally to their coefficients, C D =0.91 and C L =0.20, respectively.
By the definition of cross-product of position and force vectors, RÂ F ¼ ‖R‖ ‖F‖sin θn, the magnitude of the moment, ‖R × F‖, is reduced whenever the magnitude of either position vector or force vector is reduced, or whenever the angle between their vectors diverges from an orthogonal orientation. Consider the location of the AOR shown in Figure 2b (blue spiral arrow), which shows the AOR passing through the body of the stone rather than rotating at its edge. Acute angles between contact vectors places the AOR near the edge of the stone (e. g. grain 7 in Figure 6), whereas very obtuse angles place the AOR much closer to the stone's centre of mass (e.g. grain 14). Increasing the angle separating contact vectors on the same grain reduces the length of both position vectors (Figure 2b) as the AOR shifts towards the centre of mass, thus affecting the magnitude of the moment balance for resisting and driving forces.
Consider again increasing the bearing angle of the POR from a 2D position, which reduces the magnitude of the drag force vector projection on the POR. The reduced drag force projection (Figure 2a,b) results in an entrainment vector, F EP ¼ F DP þ F LP , that changes in both magnitude and direction within the POR (i.e. F DP decreases as F LP remains fixed). As bearing angle increases, the angle between the entrainment vector, F EP , and its position vector, R EP (shown in Figure 2b), increases from slightly acute to slightly obtuse as the drag force projection decreases in magnitude. This result corresponds to increasing angle θ in the cross-product ‖R EP ‖ ‖F EP ‖sin θn from an acute angle to the maximum ‖R EP ÂF EP ‖ value when vectors become orthogonal only to decreasing again with increasing obtuse angle. Any of these changes in vector magnitude, direction, or angle between vectors requires adjustments in the boundary shear stress to bring the moment back into static equilibrium.
When combined, these geometric factors affecting moment balance form the bulk of entrainment parameter differences between 2D and 3D entrainment models (see Figure 4). Moreover, using a vector-based 3D framework to calculate static moments from scanned images of the actual grain structure yields precise results with respect to the rigid body mechanics of the 3D entrainment model. Thus any 2D-to-3D differences in entrainment values on the same grain are primarily attributed to a measurable, geometric bias, which is inherent in any 2D model. The preceding arguments also demonstrate that friction angles obtained in the field are largely an amalgamation of angles for pivot, bearing and tilt.

Sources of error in the 3D entrainment model
Inaccuracies in the 3D model are primarily due to: (i) simplifying assumptions for the fluid mechanics of the 3D model; (ii) errors related to XCT scanning and image processing; (iii) missing grains and contact points from the images; and (iv) the rule set used to define viable pairs of contact vectors. We attempted to minimize 2D-to-3D model differences with respect to fluid mechanics used in the 3D model so that model differences would primarily be attributed to the rigid body mechanics. Indeed, since fluid mechanics developed for our 3D entrainment model are a discrete analogue of the continuous functions used in the 2D Kirchner et al. (1990) model, they are mechanistically comparable. In examining both models using the coarse sample of R1, without cohesive influence, and using the same drag and lift coefficients, we are effectively able to compare the influence of a 3D framework relative to a 2D framework by controlling for continuum influences.
Errors attributed to XCT scanning and image processing includes quality and resolution of the scan, segmentation algorithm and steps taken to segment grains from matrix, and separation methods and the unavoidable grain erosion that ensues. Although their cumulative error is relatively minor, when compared with geometric errors, they can lead to missing grains and contact points, which can only be quantified through a detailed error analysis that is beyond the scope of this work. Smaller grains are usually the culprits in this regard since larger stones tend to retain their shape and proximity to neighbouring grains post separation.
We developed a rule set for viable contact pairs based on fairly relaxed assumptions with respect to the limits of POR orientation. More investigations into these limits using carefully controlled experiments should establish a physically rigorous set of rules that are likely to be more restrictive. Differences between the new and current rule sets would reveal any error that exists with the current rule set.

Grain orientation and representative axis
Our results show that, for our samples, the Kirchner 2D model typically predicts a lower τ * cr than our vector-based 3D model. The lower 2D prediction can be explained by two simplifications that are made by the Kirchner et al. (1990) model. First, grain geometry is simplified down to a case of a circular grain sitting on a bed of uniform circular grains. Second, grains are assumed to pivot around a single point in the downstream direction. The Kirchner model, and many similar 2D models, represents grains as a circle with a diameter equal to the grain's b-axis. However, natural grains are almost never spherical, instead having a variety of complex forms that are usually represented as an ellipsoid with principal axes: a (longest), b (intermediate) and c (shortest). The variety in grain size and shape means that the protrusion of any grain is determined by the grain's orientation and position relative to those of its neighbouring grains.
Numerous studies have evaluated the preferential orientation of water-worked grains in both experimental and field scenarios. Aberle and Nikora (2006) found that for flume-derived armour surfaces the majority of grains were oriented with their long axis aligned with the flow, which was attributed to interaction dynamics between the bed and constant discharge in the experiment. In contrast, the largest grains, which did not move during water-working and armouring, were arbitrarily oriented (Aberle and Nikora, 2006). Other studies have also found flow parallel orientation of the long axes of grains, although a range of causal mechanisms have been proposed, including low transport rates or the presence of static armour layers (Hodge et al. 2013;Powell et al.2016).
Our results cast doubt on the suitability of using the b-axis as a proxy for protrusion when estimating 2D entrainment thresholds, or in roughness calculations more generally. In these situations we recommend use of the c-axis instead of the b-axis, since doing so reduces the discrepancy between τ * cr predicted by the 2D and 3D models by 62.8% as relative difference with respect to c-axis calculations (Figure 4f). Conversion to c-axis could be as simple as estimating a mean for c-to-b axis ratios obtained from a Wolman pebble count (Wolman, 1954) to use as a size multiplier on b-axis grain size distributions.

Grain entrainment angles
The pivot (and other) angles that we calculate from our 3D XCT data and model are not limited to a defined range of values as the model does not rely on trigonometric functions, instead using a 3D vector space for calculations. In contrast, the Kirchner et al. (1990) model limits pivot angles to between 0°and 90°(e.g. see Table I). Consequently, we sometimes record pivot angles that are greater than 90°, for example for a grain where the contact points are above the grain's centre of mass. Such a grain will sit mostly below the mean bed surface, but it is still subject to a component of the driving forces. Although such grains are unlikely to be very mobile, the 3D model enables these grains to be assigned a τ * cr . Therefore, when calculating the range of τ * cr values over a gravel surface, we can incorporate all grains for which a proportion of their surface can be seen from directly overhead. The inclusion of what might be considered outliers in a 2D grain entrainment implementation accounts for some of the extreme values and scatter in the data reported for the 3D results (Figures 4 and 5). Being able to calculate τ * cr for these grains shows that the 3D model can be applied beyond the realm of idealized situations.
To assess the impact of cohesion we have used a simple empirical formulation, which shows the potential importance of the term and how easily it can be incorporated into the 3D model. Cohesion effects will be a function of the matrix GSD, and so further work is needed to develop and test the cohesion term including a test for field site dependence.

Fluvial response within the roughness layer
We used a logarithmic velocity profile assuming zero velocity at the local mean bed elevation to facilitate model comparison; however, this assumption ignores interfacial flows that occur within the roughness layer. A more appropriate profile model to describe flow over a coarse gravel bed assumes either a 3D ENTRAINMENT MODEL WITH APPLICATION TO CT SCANNED RIVERBED SEDIMENT linear (Nikora et al. 2004;Mignot et al. 2009) or a power law (Sarkar and Dey, 2010;Dey and Das, 2012) profile within the roughness layer that increases from zero near the bed troughs coupled with a logarithmic profile upon entering the surface layer. By extracting roughness layer elevation boundaries from images at a distance, say D 84 , upstream of the grain, a profile may be constructed to accommodate a roughness layer that is unique to each surface grain. These variations in local hydraulic properties affecting the profile should generate more appropriate entrainment metrics across all surface grains in the sample.
The potential for grain mobilization is spatially dependent on surrounding sediment; upstream particles within close enough proximity are found to strongly affect grain entrainment (Measures and Tait, 2008). For our 3D model we developed an exposure factor that is a physically-based weighted mean of exposed fraction of grain area. XCT scanned images are rotated to provide multiple viewing angles of surface grains to mimic turbulence fluctuation of flow velocity over the bed, where the horizontal component of the vector is assumed to contribute to drag force during grain entrainment. Measures and Tait (2008) used surface images of bed sediment to isolate individual grains and to extract upstream projected areas and their midpoint elevations, and then used these metrics for individual grain calculations. They considered two mechanisms of sheltering from upstream grains: direct sheltering, which reduces exposed grain area; and remote sheltering, which modifies flow due to particle protrusion (Measures and Tait, 2008). Both mechanisms affect EF in our 3D entrainment model, since a sheltering grain increases EF values with increased distance from the object grain.

Applications of 3D models
Although vector-based 3D entrainment models yield more realistic rolling motion physics than their 2D counterparts, their downside is the impracticality of collecting the 3D scanned images that are necessary for model parametrization. The time and costs involved in obtaining and processing sediment scans of riverbed samples are currently prohibitive for routine monitoring applications. The real value for developing vector-based 3D models to interpret scanned images of grains is their potential for improving existing 2D models (e.g. switching from baxis to c-axis entrainment calculations). Furthermore, our finding that grain entrainment is primarily controlled by projection suggests that high-resolution 2.5D measurement techniques such as terrestrial laser scanning (TLS) offer exciting opportunities for in situ, field-based estimates of entrainment using exposure as a proxy for projection.
Here we developed a 3D model that utilized simple constructs typically used in fluid mechanics of 2D entrainment models to facilitate comparison with a simple 2D model. Future development should focus on improving vector-based 3D modelling with respect to the fluid mechanics; three such areas of enhancement are computational fluid dynamics simulation over bed topography (Hardy, 2006;Hardy et al. 2005;Lane et al. 2004), coupling interfacial velocity profiles for the roughness layer with logarithmic profiles for the surface layer (Nikora et al. 2004;Mignot et al. 2009;Sarkar and Dey, 2010;Dey and Das, 2012;Blois et al. 2014;Cooper et al. 2018), and entrainment effects due to turbulent pressure fluctuations (Amir et al. 2014;Cooper et al. 2018;Vollmer and Kleinhans, 2007). In doing so, a vector-based entrainment model representative of 3D rolling motion physics could serve as a benchmark from which enhancements to existing 2D models can be developed and tested. Kirchner et al. (1990) advocated a better understanding of the variability within gravel bed surfaces, not simply between different locations. We have developed a vector-based 3D grain entrainment model that is parametrized using high-resolution XCT data of individual grains within water-worked gravel beds. Our new approach explicitly accounts for variability by making no simplifying assumptions about the grain geometry and location of the contact point(s) for pivoting, instead relying upon accurately measured grain geometries and the spatial position all contact points. Our 3D entrainment model does not rely on pivot angles to calculate critical entrainment shear stress; however, in addition to pivot angles, bearing and tilt angles for the plane of rotation are easily calculated once grain entrainment is established.

Conclusion
The 3D model presented here has made geometric entrainment formulations relevant and applicable to real-world data, removing the need for idealized situations and gross assumptions. In the event that one cannot fully resolve 3D subsurface grain geometry, we show that use of the grain c-axis, rather than the oft-used b-axis, for 2D grain entrainment calculations results in critical shear stress estimates that more closely resemble those calculated using the full 3D implementation. Furthermore, the new model ultimately reveals that entrainment is predominantly controlled by (i) grain projection, (ii) plane of rotation bearing angle, and (iii) cohesion force due to grain contact with a matrix of fines. Whilst protrusion (and as we have shown here, projection) has long been recognized as a dominant control on critical shear stress, the impact of bearing angle and cohesion has not yet been accounted for. We demonstrate major sources of geometric error inherent in 2D entrainment models to show where improvements can be made in 2D model development; the magnitude of this error varies with changes to orientation for the plane of rotation.
Grain entrainment is a function of both the overall structure of the gravel bed and the overlying flow. Our model has significantly improved the representation of the former, but still uses a simplified logarithmic flow profile to represent the latter. Lamb et al. (2008) noted that local average velocity is not the only relevant velocity scale of interest and that the role of turbulent fluctuations on grain entrainment should be considered, with turbulent flow and pressure differentials aiding entrainment of grains from a bed (Vollmer and Kleinhans, 2007). The structure of our vector-based 3D grain entrainment model is such that it could be developed in the future to address the impact of flow properties on grain entrainment. By improving flow properties to better replicate local hydraulics on vector-based 3D models, entrainment modelling of 3D scanned riverbed grains has the potential for use as a benchmark for developing and testing enhancements in existing 2D models.
shown to be equivalent to projecting both vectors, R and F, onto the plane of rotation (POR) before calculating the scalar moment. The scalar triple product effectively produces a signed magnitude of the resulting moment of the system. First, resolve arbitrary vectors for resultant force, F, and position of applied force, R, into two components: one parallel to the AOR, F 1 and R 1 , and one orthogonal to the AOR, F 2 and R 2 , which lie within the POR ( Figure A1). Then the scalar triple product is given by The first term of distributed cross-products, R 1 × F 1 , in the second equality is zero because R 1 and F 1 are parallel vectors. The second and third terms, R 1 × F 2 and R 2 × F 1 , in the third equality are zero because R 1 and F 1 are parallel to unit vector λ; hence the cross-product of either of these two vectors with any other vector is obviously orthogonal to λ. The remaining nonzero term is a scalar, equivalent to projecting force and position vectors onto the POR, R 2 and F 2 , evaluating their vector cross-product, R 2 ×F 2 , and then calculating the signed magnitude of the result, sgn{R 2 ×F 2 }‖R 2 ×F 2 ‖, the latter of which is the consequence of taking the dot product of the unit vector λ with the resulting moment vector R×F.
Appendix B: Exposure Factor (EF) Model and its Adjusted Value (EF adj )

Exposure factor
The exposure factor, EF, for our 3D model is the functional equivalent of grain protrusion used in a 2D model, but is necessary since the elevation of sheltering effects from upstream grains is not laterally continuous across the face of the downstream grain. The exposure factor, EF∈[0,1], is defined as the average proportion of the face of the grain that is exposed to the flow when the grain is viewed over multiple angles. We use the turbulence of flow over the bed as a physical basis for developing a 3D exposure factor. One way of picturing flow turbulence is to imagine a camera attached to the flow velocity vector, V, that points in the same direction as the vector. Turbulence is a departure from the dominant downstream direction of flow, and so provides different viewing angles, θ, of the bed surface (see colour bands in Figure B1a). Viewed from overhead at θ=0°, the entire stone is seen (blue band), whereas looking downstream at θ=90°, only a fraction of the stone is exposed (red band). Under turbulent flow regimes, a sphere experiences viscous shear and pressure forces up to the moment of flow separation at about 115°beyond the front of the sphere (Celik et al. 2014;Constantinescu and Squires, 2004); thus we can use overhead views of exposed areas in EF calculations. Let ER θ be the exposure ratio defined as the visible area of exposed stone (solid colour bands) divided by the total area (A θ , thick dashed colour lines) of the stone from the same viewing angle θ. Now consider the magnitude of the horizontal component of the velocity vector, ‖V sinθ‖, as the speed of river current at an arbitrary point along the velocity profile, u, which is used in calculating drag force. Observing that ðV sin θÞ 2 ¼ V·V sin 2 θ ¼ ‖V‖ 2 sin 2 θ ¼ ‖V sin θ‖ 2 ; the drag force on the stone is proportional to the product of dynamic pressure and stone area at the same viewing angle Figure A1. Image of a plane of rotation (blue plane) and its orthogonal axis of rotation (dotted line). Arbitrary force and position vectors, F and R, that induce a moment R×F in 3D space are resolved into components parallel to the axis of rotation, F 1 and R 1 , and perpendicular to the axis of rotation, F 2 and R 2 , spanning the plane of rotation. The unit vector for the axis of rotation is λ. Redrawn after Beer et al. (2013).
[Colour figure can be viewed at wileyonlinelibrary.com] Figure B1. (a) Development of the exposure factor, EF. For three different view angles, solid coloured bands show the areas of viewable stone, whereas the dashed coloured lines indicate total area of the stone that could be viewed in the case of no upstream sheltering grains. Viewing angle range from 0°(directly overhead) to 90°(in the direction of flow). The exposure ratio, ER θ , is the ratio of these two area values. (b) The adjusted exposure factor, EF adj , is necessary because we are only interested in the relative exposure of the areas of the grain that are above the elevation at which the flow velocity is zero. To calculate EF adj , we partition the face of the stone into two initial regions: above the mean bed elevation where the velocity profile is nonzero (white and light-grey areas) and below mean bed elevation where flow velocity is zero (dark-grey area). The upper region is further partitioned into the hidden fraction (light grey) and the exposed fraction (white) of the nonzero velocity profile. [Colour figure can be viewed at wileyonlinelibrary.com] 3D ENTRAINMENT MODEL WITH APPLICATION TO CT SCANNED RIVERBED SEDIMENT 1 2 ρðV sin θÞ 2 A θ (equivalent to 1 2 ρu 2 A from Equation (9a)). Hence we can use the magnitude of the horizontal component of the velocity vector (illustrated by colour-matched vectors each from the same velocity V; Figure B1a) as a physically-based weighting mechanism in estimating a weighted mean of ER θ , which is the exposure factor EF used in Equation (9b). Thus EF is calculated by Under idealized conditions of spherical grains, 1 2 ρ‖V‖ 2 A θ is not dependent on viewing angle, so it factors out of both summations and cancels across the fraction. However, for natural stones, which are typically not spherical, this weighting mechanism provides a good approximation as an exposure factor. The summed denominator is necessary since the sum of all the weightings used in a weighted mean must be unity. To minimize edge effects when calculating EF on scanned images, we extended the edges of basket of grains by replicating copies of a binary image of sample particles in a circular pattern around the entire basket with a slight buffered overlap around the outer edge.

Adjusted exposure factor, EF adj
The EF is the proportion of the total face of the grain that is exposed to the flow; however, for the entrainment calculation we are only interested in the proportion of the grain face that is exposed to nonzero flow velocities. Since the nonzero portion of the logarithmic velocity profile is defined as starting at the mean local bed elevation, which is generally above the base of a grain, it is necessary to rescale EF so that portions of the grain where the profile is zero are not included in the drag force calculations. Without this step, drag force calculations used in the moment balance will be underestimated, resulting in increased critical shear values. The adjustment to EF depends upon whether the fraction of stone experiencing no flow (i.e. zero velocity profile) is greater or less than the fraction that is sheltered by upstream grains (as expressed by EF). Let F u be the fraction of a stone above the mean local bed elevation where u>0, and let F u =F 1 +F 2 be the partition of this fraction such that F 2 =EF is the exposure factor (see Figure B1b). Then F 1 =F u −F 2 =F u −EF. Since EF represents the exposed fraction for the entire stone, the hidden fraction of F u is F 1 /F u (light-grey section) and the exposed fraction of F u is the adjusted exposure factor EF adj =F 2 /F u (white section). EF adj further reduces drag force on the stone due to the additional sheltering effects of EF, which are only partially accounted for by F u (dark-grey section). Then, from the defined partition: Now whenever 0<F u <EF≤1, EF already covers F u in its entirety, so there is no need to reduce drag force any further than already accounted for by F u . Hence, in this case EF adj =1. For computational efficiency, we combine both conditional domains for EF adj into a single, unconditional domain; thus the EF adj becomes ; Fu; EF ∈ 0; 1 ð :

Appendix C: Empirical Cohesive Force (F C ) Experimental Design and Model
Cohesive force experiments were designed to determine the effect of burial depth, sand size and clay fraction on the vertical tensile force required to extricate a 25 mm diameter spherical glass marble from a sand-clay mixture. The chosen marble size represents the lower end of the D 50 size range from field data (Hodge et al. 2013). The marble was also pulled from pure clay to establish a baseline from which treatment effects could be modelled relative to the effects of pure clay. Treatments used in the experiment were burial depth (4 mm, 8 mm and 12 mm), size range of sand ((0.1,0.3] mm, (0.5,1.0] mm and (1.0,2.0] mm), and clay fraction of the sand-clay mixture (0, 0.25, 0.50 and 0.75) as a 3×3×4 main effects factorial design. For the experimental setup, a three parts sand-clay mixture was thoroughly mixed with one part water, and then poured into a small tray that could accommodate 15 buried marbles (three burial depths × five replicates) with sufficient spacing as to avoid surface deformation effects of neighbouring grains during force pulls. The sand-clay mixture was spread over the tray, and the surface was smoothed without compaction to ensure consistency between pours. String was glued to each marble prior to vertical pulls. Force pulls were made of each freestanding marble prior to burial to adjust their extraction force by the marble's weight. Marbles were then buried to their depth markings, and the surface was allowed to set undisturbed, yet still damp, until force measurements, which were obtained using a handheld force gauge.
Force pulls for pure clay were found to be independent of burial depth, as the simple linear regression was found to be insignificant (p=.526). Tensile force per buried surface area of the marble for all sand-clay mixtures had substantial nonlinear relationships with sand size and clay fraction. Therefore, a single power law model was developed where the response was defined as the force per surface area for the mixtures divided by the overall mean of force per surface area for pure clay (i.e. the comparison baseline). Then the zero value clay fraction was adjusted to 0.001 before log-linearizing the power law model for linear regressing over mean sand sizes (0.20 mm, 0.75 mm and 1.50 mm) and clay fractions (0.001, 0.25, 0.50 and 0.75). The significant regression of the linearized power law (p<.001,R 2 =.54 (log)) was exponentiated and multiplied through by the baseline value for pure clay force per surface area, yielding the final model ( Figure C1). Other empirical definitions could be derived for use within our vector-based 3D model. Data and model development are located on our Github page at https://github. com/NERCPATCheS/VectorEntrainment3D.

Appendix D: Algorithmic Solution for 3D Vector-Based Grain Entrainment
An algorithmic solution is provided for calculating τ * cr so that more complexity can be added into the 3D framework (e. g. adding turbulence) without deriving an analytical solution. To arrive at a value of τ cr requires finding the boundary shear stress, τ b , used in Equation (8) that satisfies the static moment balance. Calculate entrainment for each stone in the sample from top to bottom as only surface stones experience shear stress. For each stone with entrainment potential, where initial force calculations are nonzero, get a list of N contact points of neighbouring stones that are within some specified distance. Then generate a list of N 2 À Á contact pairs to construct vector pairs (CV L ,CV R ), using the rule set to test all vector pairs for viability, keeping only the pairs that are viable for entrainment calculations. Next, calculate a set of τ cr values, where each value in the set {τ cr } corresponds to each viable pair of contacts. The smallest τ cr value in the set {τ cr } and its corresponding pair of contacts (CV L ,CV R ) are the entrainment solution for that particular stone. Calculate the corresponding τ * cr value, all entrainment angles using Equations (6), and any other useful information. What follows is the algorithm used to obtain a single τ cr value for a given contact pair (CV L ,CV R ) that will be added to the set {τ cr } from which the entrainment solution for a given grain is τ cr =min{τ cr }.
Once we reach an interval [τ 1 ,τ 2 ] that generates m 1 ×m 2 <0, we then use a binary search algorithm as follows: (i) split the interval [τ 1 ,τ 2 ] in half at the mean,¯τ ¼ ðτ 1 þ τ 2 Þ=2, and generate moment intervals ½m 1 ; m¯τ and ½m¯τ ; m 2 corresponding to ½τ 1 ;¯τ and ½¯τ ; τ 2 , respectively; (ii) test each half interval ½m 1 ; m¯τ and ½m¯τ ; m 2 to find which half contains zero; (iii) choose a new [τ 1 ,τ 2 ] that corresponds to the moment interval that contains zero going back to step (i) and repeat until maxfjm 1 j; jm 2 jg≈0; (iv) calculate¯τ from the final [τ 1 ,τ 2 ] saving its value to the set {τ cr } for the current (CV L ,CV R ) pair. Get the next (CV L ,CV R ) pair for the current grain, setting [τ 1 ,τ 2 ] =[0,Δτ b ] at n=1, and repeat the algorithm for all valid (CV L , CV R ) pairs of the current stone. Repeat the algorithm for all remaining stones in the sample to obtain a τ * cr value plus angles and other metrics for each stone with entrainment potential. A flow chart detailing the steps taken in the algorithmic solution is provided ( Figure D1).
Binary search algorithms are computationally efficient as they have a geometric rate of convergence in narrowing [τ 1 ,τ 2 ] intervals to obtain a τ cr value. By initializing the search for zero in [m 1 ,m 2 ] over n=k,k+1,…, the [τ 1 ,τ 2 ] interval used in the binary search is obtained in just a few iterations. Figure D1. Flow chart of algorithmic solution for 3D vector-based grain entrainment.
3D ENTRAINMENT MODEL WITH APPLICATION TO CT SCANNED RIVERBED SEDIMENT Combining both of the these enhancements with parallel processing the for-loops over contact pairs (CV L ,CV R ), this algorithm is highly efficient. The Λ values in Equation (15c) indicated the potential for τ cr <0, which is possible due to the relaxed β<|90°| assumption of the current rule set for viable (CV L ,CV R ) pairs. For our analysis, we chose to ignore any negative τ cr values as extraneous information.