Entrainment and suspension of sand and gravel

The entrainment and suspension of sand and gravel are important for the evolution of rivers, deltas, coastal areas, and submarine fans. The prediction of a vertical profile of suspended sediment concentration typically consists of assessing (1) the concentration near the bed using an entrainment relation and (2) the upward vertical distribution of sediment in the water column. Considerable uncertainty exists in regard to both of these steps, especially the near-bed concentration. Most entrainment relations have been tested against limited grainsize-specific data, and no relations have been evaluated for gravel suspension, which can be important in bedrock and mountain rivers. To address these issues, we compiled a database with suspended sediment data from natural rivers and flume experiments, taking advantage of the increasing availability of high-resolution grain size measurements. We evaluated 12 dimensionless parameters that may determine entrainment and suspension relations and applied multivariate regression analysis. A best-fit two-parameter equation (r2 = 0.79) shows that near-bed entrainment, evaluated at 10 % of the flow depth, decreases with the ratio of settling velocity to skin-friction shear velocity (wsi/u∗skin), as in previous relations, and increases with Froude number (Fr), possibly due to its role in determining bedload-layer concentrations. We used the Rouse equation to predict concentration upward from the reference level and evaluated the coefficient βi , which accounts for differences in the turbulent diffusivity of sediment from the parabolic eddy viscosity model used in the Rouse derivation. The best-fit relation for βi (r2 = 0.40) indicates greater relative sediment diffusivities for rivers with greater flow resistance, possibly due to bedform-induced turbulence, and larger wsi/u∗skin; the latter dependence is nonlinear and therefore different from standard Rouse theory. In addition, we used empirical relations for gravel saltation to show that our relation for near-bed concentration also provides good predictions for coarse-grained sediment. The new relations extend the calibrated parameter space over a wider range in sediment sizes and flow conditions compared to previous work and result in 95 % of concentration data throughout the water column predicted within a factor of 9. Published by Copernicus Publications on behalf of the European Geosciences Union. 486 J. de Leeuw et al.: Entrainment and suspension of sand and gravel


Introduction
The suspension of sediment by water plays a critical role in the dynamics of rivers, river deltas, shallow marine environments, and submarine canyons and fans. For example, suspended sediment dominates the load of lowland rivers and builds land in subsiding river deltas and coastal landscapes (Ma et al., 2017;Syvitski et al., 2005;Wright and Parker, 2004a). The transport of sediment on the continental shelf is dominated by the suspension of mud and sand due to wave and current bottom traction (Cacchione et al., 1999;Nittrouer et al., 1986). Suspended sediment provides the negative buoyancy of turbidity currents that move sand and gravel to the deep sea. The suspension of gravel is important in large floods, such as outburst floods (Burr et al., 2009;Larsen and Lamb, 2016), and in steep mountain canyons (Hartshorn et al., 2002), where it can contribute to bedrock erosion (Lamb et al., 2008a). Suspended sediment transport is also important in landscape engineering, such as river restoration (Allison and Meselhe, 2010), fish habitat (Mutsert et al., 2017), and the capacity of dams and reservoirs (Walling, 2006). The balance between entrainment and deposition from suspension determines patterns of deposition and erosion in these environments and therefore controls landform morphology and stratigraphic evolution (Garcia and Parker, 1991;Paola and Voller, 2005). To predict suspended sediment flux across these environments, we need robust theory for the entrainment of sediment from the bed and the vertical distribution of suspended sediment in the water column.
We focus here on understanding the suspended sediment load of cohesionless grains that are entrained from the bed (i.e., suspended bed material) rather than wash load. Most models for sediment suspension are based on the application of Rouse theory (Rouse, 1937;Vanoni, 1946), where C (L 3 L −3 ) is the volumetric sediment concentration at elevation (z; L) above the bed, C a (dimensionless) is the reference near-bed concentration at z = a (unit: L), and H (L) is the flow depth. P denotes the dimensionless Rouse number, in which w s (L T −1 ) is the particle settling velocity, κ is the dimensionless von Karman constant of 0.41, u * (L T −1 ) is the bed shear velocity, and β is a dimensionless factor that accounts for differences between the turbulent diffusivity of sediment and the parabolic eddy viscosity model used to derive the Rouse equation (e.g., Graf and Cellino, 2002). Although more sophisticated models exist, some of which abandon the Rouse theory entirely in favor of a more rigorous turbulence model (Mellor and Yamada, 1982), the Rouse equation remains a useful and tractable approach for modeling and field application (Graf and Cellino, 2002;van Rijn, 1984;Wright and Parker, 2004b). The Rouse equation was derived assuming an equilibrium suspension whereby the upwards volumetric flux of sediment per unit area due to turbulence (F z ; L T −1 ) is balanced by a downwards gravitational settling flux (Cw s ) ( Fig. 1) (Rouse, 1937). It predicts the shape of the concentration-depth profile, with a greater gradient in concentration for larger P (Fig. 1). By mass balance, the difference between entrainment, F za , and settling, C a w s , fluxes per unit area near the bed (at z = a) controls the bed deposition rate, dz b /dt (L T −1 ) (Parker, 1978), i.e., where λ (dimensionless) is the bed porosity, and E s ≡ F za /w s is a dimensionless sediment entrainment rate or entrainment parameter (Garcia and Parker, 1991). At steady state, Eq. (3) reduces to C a = E s ; thus, the near-bed entrainment rate, E s , is necessary to predict both the vertical distribution of suspended sediment at steady state (i.e., C a in Eq. 1) and the rate of erosion and deposition for disequilibrium suspensions (Eq. 3). The application of Eqs.
(1)-(3) requires the specification of β, E s and a, and the approach in previous work has been to identify important dimensionless variables and fit data from flume experiments and natural rivers to these variables (e.g., Smith and McLean, 1977;van Rijn, 1984;Garcia and Parker, 1991). Owing to differences in datasets analyzed and the dimensionless variables explored, there is considerable deviation in the form of previous models and their predictions (Figs. 2 and 3, Tables 1 and 2). For example, Einstein (1950), Engelund and Fredsøe (1976), van Rijn (1984), and Smith and McLean (1977) used a relation for a that depends primarily on grain diameter (D; L) and secondarily on u * /w s or Shields number, τ * = τ b / (ρ s − ρ w ) gD (where τ b (ML −1 T −2 ) is bed stress, ρ s (M L −3 ) is sediment density, and ρ w (M L −3 ) is fluid density). Their rationale was based on the idea that there is a well-mixed near-bed zone of bedload transport and that suspended sediment is entrained from this zone (Fig. 1b). Thus, they developed relations for a that scale with the height of bedload saltation. In contrast, Garcia and Parker (1991) and Wright and Parker (2004b) argued for a simpler approach and used a reference height that is a small fraction of the flow depth, and they proposed a = 0.05H as a useful reference height, with no dependence on D.
In practice the factor β is often neglected (i.e., β = 1) under the assumption that the parabolic eddy viscosity assumption used to derive Eq. (1) represents the turbulent transport of fine suspended sediment well. β < 1 in Eq. (2) has been attributed to damping of turbulence due to sediment-induced stratification, which alters the eddy viscosity (Einstein, 1955;Graf and Cellino, 2002;Wright and Parker, 2004a), and to flocculation, which increases the settling velocity of finegrained sediment (e.g., Bouchez et al., 2011;Droppo and   Existing relations for the factor β (Eq. 2) as a function of (a) shear velocity normalized by setting velocity u * /w s (van Rijn, 1984;Graf and Celino, 2002;Wright and Parker, 2004b), (b) near-bed concentration normalized by slope C a /S (Wright and Parker, 2004b), and (c) w s /u * (H /D) 0.6 as proposed by Santini et al. (2019). In panel (a), the Wright and Parker relation (denoted W & P) used C a calculated from the entrainment relation of Wright and Parker (2004b). Ongley, 1994). Some formulas (van Rijn, 1984;Santini et al., 2019) and datasets (Graf and Cellino, 2002; indicate β > 1, which implies enhanced mixing of sediment relative to momentum (Table 1; Fig. 2). Graf and Cellino (2002) showed that turbulence generated by bedform roughness results in better sediment mixing and thus in a higher β. Field and flume data (Chien, 1954;Coleman, 1970;van Rijn, 1984) indicate that β is greater than unity for the coarse grain size fraction of the suspended material (Greimann and Holly, 2001). Nielsen and Teakle (2004) argued that the Fickian diffusion model in the derivation of the Rouse equation is not valid for steep concentration gradients, resulting in β > 1. Most previous relations show a trend of decreasing β with increasing u * /w s (van Rijn, 1984;Graf and Cellino, 2002;Wright and Parker, 2004b) (Fig. 2a), which could be due to turbulence damping associated with high concentration suspensions. Wright and Parker (2004b) showed that β is also a function of reference concentration divided by slope (C a /S) (Fig. 2b), which they attributed to sediment stratification. Santini et al. (2019) found that β is a function of u * /w s and the ratio between flow depth and bed grain size (H /D) (Fig. 2c).
Compared to β, even larger uncertainty exists in the dimensionless entrainment rate, E s . In previous work, E s was assessed by measuring the near-bed concentration for equilibrium suspensions, since E s = C a under those conditions (i.e., dz b dt = 0 in Eq. 3). Previous work found that C a was a function of bed stress and grain size ( Fig. 3; Table 2) (Einstein, 1950;Engelund and Fredsoe, 1976;Smith and McLean, 1977;van Rijn, 1984;Cedik and Rodi, 1984;Akiyama and Fukushima, 1986;Garcia and Parker, 1991;Wright and Parker, 2004b). While holding bed stress constant, most relations predict smaller C a with increasing grain size across the sand regime, as expected, but some surprisingly show increasing C a for coarse sands and gravel (Akiyama and Fukushima, 1986;Garcia and Parker, 1991;Figure 3. Existing relations for near-bed suspended sediment concentration under equilibrium conditions (i.e., E s = C a ) at the reference level, z = a, specified by each study as a function of (a) grain size and (b) bed stress. In panels (c) and (d) the predicted concentrations are interpolated to a common reference level, a = 0.05H , using the Rouse equation with β = 1 to better compare models. for 0.1 < w s u * < 1 u * /w s Graf and Cellino (2002) β = 3 10 + 3 4 w s u * for 0.2 < w s u * < 0.6 and no bedforms u * /w s β = 1 + 2 w s u * 2 for 0.1 < w s u * < 1 and bedforms Wright and Parker (2004b) Rose and Thorne, 2001) Wright and Parker, 2004b) (Fig. 3a, c). Due to the greater weight and settling velocity of larger particles, this behavior is unrealistic and likely occurs because these coarse particles are outside the data range used to fit the relations. Despite the importance of gravel suspension in steep mountain rivers and megafloods, data do not exist to evaluate the mod-els in the gravel regime. All relations show increasing C a with bed stress, as expected, but there are orders of magnitude differences in the predictions for C a (Fig. 3b). Part of this deviation is due to differences in the reference level (z = a) (Fig. 3b), but even when using a common reference level (a = 0.05H ; Fig. 3d), significant variability still exists. Engelund and E s = 0.65 Brownlie (1981):  Some of the variance is also due to extrapolating the models beyond the range in which they were calibrated; for example, van Rijn (1984) predicts concentrations greater than 100 % at high shear stresses. We revisited the problem of sediment entrainment and the suspension of cohesionless bed sediment by compiling a large database of sediment-size-specific data for bedsediment mixtures, testing existing relations against the database, and proposing improved relations for E s and P . Several existing formulas for E s are applicable only to uniform-sized bed sediment (Akiyama and Fukushima, 1986;Celik and Rodi, 1984;Einstein, 1950;Engelund and Fredsøe, 1976;van Rijn, 1984;Smith and Mclean, 1977) ( Table 2). However, given the strong control of grain size on near-bed concentrations, accurate formulas likely need to make grain-size-specific predictions for sediment mixtures (Garcia and Parker, 1991;McLean, 1992;Wright and Parker, 2004b). Our database takes advantage of high-resolution grain size measurements using laser diffraction Gitto et al., 2017;Haught et al., 2017;Santini et al., 2019), which allows a single concentration profile to be separated into many grain-size-specific concentration profiles. Our dataset also contains a wide range grain sizes, extending into the silt regime (median bed sizes range from 44 to 517 µm), and expands on the number of field measurements compared to previous efforts. We are unaware of studies on C a for gravel, but data do exist for saltation heights, velocities, and bedload fluxes for gravel, and calibrated relations for these variables also exist (Chatanantavet et al., 2013;Sklar and Dietrich, 2004). To attempt to better constrain suspension of very coarse sand and gravel, we used existing semiempirical saltation theory for gravel to check for consistency between sand suspension data and what might be expected for near-bed gravel concentrations.

Suspended sediment profiles
Based on previous theory, we searched for available datasets from rivers and flume experiments that had suspended sediment profiles (C(z)), depth-averaged flow velocity (U ), flow depth (H ), channel bed slope (S), bed material grain size (D), and the grain size distribution of the bed material and suspended sediment samples (Tables 3 and S1 in the Supplement). Some of the experimental studies used a narrow grain size distribution, and, like previous work, we assumed that the sediment distribution from these studies was uniform. Many of the older datasets were used in empirical regressions from previous relations (e.g., Garcia and Parker, 1991;van Rijn, 1984;Wright and Parker, 2004b). In addition, we used a river dataset from the Yellow River (Moodie, 2019), which provides a fine-grained end-member. In total, our database contains 180 concentration profiles from eight rivers and 62 profiles from six different experimental studies. We analyzed only the grain fractions coarser than 62.5 µm (i.e., sand). The mud fraction was present in the bed only in small amounts, and following previous work, mud was assumed to require a different approach, potentially due to supply limitation (Garcia, 2008), cohesion, or flocculation. Grain size distributions in older studies were determined from sieve analysis of bed material and suspended sediment samples. The more recent studies used laser diffraction techniques, which have the advantage that a larger number of grain size classes can be distinguished. We calculated the grain-size-specific suspended sediment concentration (C i ) using where f i (dimensionless) is the mass fraction of grains of the ith size and C tot is the total suspended sediment concentration for all sizes. In addition, we computed the D 50 (median grain size) and D 84 of the bed material using linear interpolation between the logarithm of D and the cumulative size distribution. Concentration profiles in the database typically contain three to eight measurements in the vertical dimension. The Rouse profile was fitted to the profile data for each grain size class in log-transformed space using linear least squares to find P i (Fig. 4a). Confidence bounds (68 %; 1σ ) for the fitted coefficients were obtained using the inverse R factor from QR decomposition of the Jacobian. Data were excluded from further steps in the analysis if the ratio between the upper and lower bound of the confidence interval was greater than 10, as these data do not follow a Rouse relation for unknown reasons (e.g., measurement error) and would appear as sparse outliers. We also eliminated points with P i < 0.01, as these profiles are nearly vertical, which poorly constrains settling velocity to be very near zero. Of the data analyzed, 201 points (15 %) were excluded based on these criteria. Some studies (Bennett et al., 1998;Muste et al., 2005) have shown that P (or β) can vary over the flow depth, but this effect cannot be incorporated into our approach; instead we found one value of P i that best fit the concentration profile for each grain size class.
We used the Rouse equation (Eq. 1) for each grain size class to extrapolate or interpolate the concentration to a reference level at 10 % of the flow depth; i.e., we set a = 0.1H and C ai = C i (z = a). Extrapolation to very near the bed can be difficult due to large concentration gradients and poorly constrained near-bed processes such as interactions with bedforms. However, a reference level that is too far away from the bed may poorly capture the exchange of sediment with the bed. Previous researchers used reference levels that were either at some fraction of the flow depth (Akiyama and Fukushima, 1986;Celik and Rodi, 1984;Garcia and Parker, 1991;Wright and Parker, 2004b) or that were related to the bed roughness height (Einstein, 1950;Engelund and Fredsøe, 1976;Garcia and Parker, 1991;Smith and Mclean, 1977) ( Table 2). We explored the collapse of the entrainment data for both types of reference levels. For each reference level,  we fit our preferred and best-fitting two-parameter entrainment relation (as described in the Results section) to the data with u * skin /w si as the first parameter and the Froude number as the second parameter. We found that a reference level at a fraction of the flow depth gave a better collapse of the entrainment data than a reference level related to the saltation layer height. Furthermore, we also tested different flow depth fractions and found that the fit improved, in a least-squared sense, as the reference level moved to a larger fraction of the flow depth (Fig. 4b). However, there is little change in r 2 once the reference level is higher than ∼ 10 % of the flow depth. Therefore, we used a reference level at 10 % of the flow depth for all results shown below. For sediment mixtures, the grain-size-specific near-bed concentration is partially controlled by the fraction of each grain size class in the surface bed material. To account for this effect, Garcia and Parker (1991) introduced an entrainment rate (E si ) for each grain size class that is linearly weighted by the fraction of that material in the bed: where C ai is the near-bed concentration of that grain size class and F bi is the mass fraction of bed material that falls in that grain size class. For uniform sediment, the entrainment rate (E si ) is equivalent to the near-bed concentration (C ai ).

Independent parameters and profile fitting
Here we review the independent parameters that we evaluated for dependencies with a dimensionless entrainment rate, E si , and for β i in the Rouse number. The primary group of parameters describes the ratio between bed stress and grain size or grain settling velocity. These parameters include the ratio between shear velocity and settling velocity (u * /w s ), whereby we evaluated the total shear velocity as assuming steady, uniform unidirectional flow. Others have proposed that entrainment depends on the skin-friction portion of the total shear stress, u * skin , rather that the total shear stress. To estimate u * skin , we used the Manning-Strickler relation, where k s = 3D 84 is the grain roughness on the bed, H sk is the depth due to skin friction, and u * skin = √ gH sk S (e.g., Wright and Parker, 2004b). To calculate the particle settling velocity, we followed Ferguson and Church (2004) for each grain size class, in which R = (ρ s − ρ w ) /ρ w is the submerged specific density of sediment, ν (L 2 T −1 ) is the kinematic viscosity of the fluid, C 1 = 18 and C 2 = 1 are constants set for natural sediment, and D i is the grain diameter within the size class of interest. Another parameter that relates to the ratio between bed stress and gravity acting on the grains is the Shields number, and we again assumed steady, uniform flow to find τ b = ρgH S. Similar to the shear velocity, it is also possible to calculate a Shields number for the skin-friction component of the total shear stress, where τ skin = ρu 2 * skin by definition. Shields numbers can be rewritten in terms of u * /w si through the use of a particle drag coefficient, which we also evaluated (i.e., τ * = The next group of parameters describes dimensionless particle sizes, including the particle Reynolds number, Likewise, this parameter can also be calculated with the skinfriction component of the shear velocity, A particle Reynolds number can be defined without shear velocity as For sediment mixtures, the relative particle size might play a role due to hiding and exposure effects (Garcia and Parker, 1991;Wright and Parker, 2004b); this effect can be captured with D i D 50 . Sediment-induced density stratification can decrease entrainment by damping near-bed turbulence, and this effect is thought to be most important in deep, lowgradient rivers (Wright and Parker, 2004b). Wright and Parker (2004b) proposed that the ratio of near-bed concentration to bed slope is a good predictor for stratification, C a S , for which they used C a at 5 % of the flow depth. Large, lowgradient rivers also have small Froude numbers and low bed slopes, so we evaluated the Froude number and slope as additional parameters. The Froude number was calculated as The entrainment rate could also be affected by turbulence or changes to the boundary layer from bed roughness or bedforms, which tend to correlate with a flow resistance friction coefficient (e.g., Engelund and Hansen, 1967), We also evaluated H /D 50 as a proxy for flow resistance due to grain roughness. In order to find relations that explain the variation in our best-fit E si and P i values from the vertical concentration profiles, we regressed the E si and P i values against the 12 variables described above (Eqs. 6-16). In some applications, like reconstructing flow conditions from sedimentary strata, it is useful to have an entrainment relation that depends on u * skin , while for forward modeling a relation based on u * is preferred. The two shear velocities are highly correlated; therefore, we explored two versions of the fit relations using either u * or u * skin , but not both at the same time. Because the Rouse parameter, P i , by definition depends inversely on u * /w si (or u * skin /w si ) (Eq. 2), we found the best-fit relations for P i rather than β i to avoid spurious correlation. We then solved for the equivalent relation for β i using the definition of P i (Eq. 2). We started the analysis by testing all models with one explanatory variable and ranked the models according to the coefficient of determination from linear least-squares regression (r 2 ) evaluated in log-log space. Next, a second parameter was used in addition to the first best-fitting parameter, and the resulting two-parameter models were ranked according to r 2 . The procedure was repeated with additional parameters until the increase in r 2 was smaller than 0.04. For the fitting of multiparameter models, we varied the exponents on each parameter in the model simultaneously to find the combination of exponents that yielded the best fit. This approach gave a higher r 2 compared to the stepwise approach used in previous work (e.g., Garcia and Parker, 1991) of first fitting the dominant variable and then fitting the secondary variables to the residuals. In addition, we tested fitting with the York method (Table S2), which gives less weight to data with large errors (York, 1968), but found only minor differences, so all results presented use the simpler linear least-squares method. All parameters were used to evaluate relations for E si and β i except for D i /D 50 , which we used only for E si since it is relevant for particle-particle interactions at the bed surface. In the results we report two versions of the best-fitting one-, two-, and three-parameter models: one version that is based on the total bed shear stress and one version that is based on the skin-friction component of the bed shear stress. Model fits using all possible combinations of the input parameters are given in Table S2.

Comparison to theory for gravel
Although gravel suspension is important in bedrock and steep mountain rivers, as well as during megafloods, we are not aware of datasets of near-bed concentration in the gravel range. Following previous work (McLean, 1992;Lamb et al., 2008a), our approach was to derive the near-bed concentration in the bedload layer and then use Rouse theory to predict that concentration at 0.1H to compare with the sand dataset (e.g., Fig. 1b). The near-bed volumetric concentration within the bedload layer can be calculated by continuity as where q b (L 2 T −1 ) is the volumetric bedload flux per unit width, H b (L) is the bedload-layer thickness, and U b (L T −1 ) is the bedload velocity. Most relations for bedload flux take the form where m and n are empirical constants, which we set to m = 5.7 and n = 1.5 (Fernandez Luque and van Beek, 1976), and τ * c is the critical Shields number at initial motion, which we set according to Lamb et al. (2008b): The bedload-layer height and velocity were determined from Chatanantavet et al. (2013). They compiled a large dataset of gravel saltation observations and found a good fit with the following relations: Equations (17)- (21) were combined with a flow resistance relation (Eq. 7, assuming no form drag) and C d = 0.7 for gravel (Lamb et al., 2017) to calculate C b in the bedload layer using a numerical iterative scheme. Manipulating the equations revealed that C b is a function of only τ * and Fr.
To calculate C a for gravel, we extrapolated the concentration profile (Fig. 1b) from the top of the saltation layer to 0.1H using the Rouse equation (Eq. 1). To obtain the Rouse number, we used our best-fit one-parameter model that uses total shear velocity (u * ) (see Table 4). We then used a wide range of input parameters (0.1 < Fr < 1 and 1 < τ * < 1000) relevant to the suspension of gravel in mountain rivers and large floods to predict a range of expected values of C a for gravel. Although Eqs. (20) and (21) have not been tested for high Shields numbers in the suspension regime, they have reasonable limiting values (U b = 0.6U , H b = H ) and provide a starting place to compare sand entrainment and gravel saltation theories.

Rouse number, P i
Equation (1) was fit to the concentration profile data wherein P i was treated as a fitting parameter. Then P i was regressed against the hypothesized controlling variables (Eqs. 11-16). Figure 5 and Table 4 show the results for P i including our best-fitting one-, two-, and three-parameter models. The variable with the most explanatory power is u * /w si (P i = (u * /w si ) −0.45 ; r 2 = 0.33). The best-fit two-parameter model that includes C f is better than the predictions of the oneparameter model (r 2 = 0.40). Going from a two-parameter model to a three-parameter model brings a smaller improvement (r 2 = 0.43). Therefore, we recommend the twoparameter model for combined accuracy and simplicity: although there are a number of models using different variables that yield similar r 2 values (Table S2), and some combinations might be preferred over others depending on the application. Equation (22) indicates that sediment is better mixed in the water column with larger u * skin /w si and with a larger bed roughness coefficient, C f . Equation (22) performs  Table 4. The plotting is done such that the exponent on the last parameter is unity; e.g., Eq. 23 in panel (b) was rewritten as P i = 0.145 1/−0.3 u * skin to allow for comparison to plots of Garcia and Parker (1991) and Wright and Parker (2004b). Yellow filled symbols are the geometric mean, and the error bars indicate the geometric standard deviation within each arbitrarily spaced bin. 3 . The two best-fitting versions of each model are shown, one with u * and one with u * skin , and the bold value is the best fitting of the two (highest r 2 ). See Table S2 for all model fits. well compared to previous relations, as is shown by a box plot of measured-to-predicted ratios (Fig. 6). For the best-fit two-parameter model, the measured-to-predicted ratio falls between 0.74 and 1.29 for 50 % of the data. Using a Rouse number with a constant β = 0.94 also provides a reasonably good fit and an improvement over several more involved relations (Fig. 6). Equation (22) can be rewritten for β i using Eq. (2), and by assuming that u * in Eq. (2) is actually u * skin , as

Number of model
Because some of the dimensional quantities appear in multiple dimensional variables and because the dimensional variables are not necessarily independent from each other, we tested for spurious correlations by rearranging the twoparameter relation for β i to isolate the dimensional dependencies on grain size and skin-friction shear velocity (Fig. 7). The data and our model show a decrease in suspended sediment mixing in the water column with larger Figure 6. Box plot of predicted versus measured Rouse number, P i , for our models and previous models. The best-fit one-, two-, and three-parameter models correspond to the equations in bold in Table 4. grain sizes, as expected (Fig. 7a). Previous relations show the same trend, but the relations of van Rijn (1984) and Graf and Cellino (2002) have stronger dependencies in the sand regime. Suspended sediments are better mixed with increasing skin-friction shear velocity (Fig. 7b). The data suggest that P i varies proportionally to u −0.4 * , whereas standard Rouse theory (Eq. 2) indicates that P i is proportional to u −1 * .

Dimensionless entrainment rate, E si
Results for the best-fit dimensionless entrainment rate are shown in Fig. 8 and Table 5, and all possible variable combinations are given in Table S1. The following one-parameter relation gives the best fit with the data for E si (r 2 = 0.61) (Fig. 8a): The Froude number was the most significant second parameter: which has a significantly better fit (r 2 = 0.79) than the bestfitting one-parameter relation (Fig. 8b). The addition of a third variable to the model gives little further improvement of the fit (Fig. 8c; r 2 = 0.80), and many of the variables used as a third parameter give a similar level of improvement (Table S1). Using Eq. (25), 80 % of the entrainment data are predicted within a factor of 3 (Figs. 8b, 9). Along with our proposed new relation (Eq. 25), we also compared the dataset against the previous relations (Fig. 9). The box plots in Fig. 9 highlight that some relations systematically underpredict (Wright and Parker, 2004b) or overpredict the entrainment rate (Garcia and Parker, 1991). In addition, previous relations have a larger spread in measured-to-predicted ratios than Eq. (25). To check for spurious correlation in the dimensionless variables, we rearranged our two-parameter entrainment relation (Eq. 25) to isolate the dependencies of entrainment rate on grain size (Fig. 10a) and skin-friction shear velocity (Fig. 10b). Our relation indicates that entrainment depends on grain size to the −3.1 power in the sand range, whereas previous relations suggest a much weaker dependence. Compared to previous work, our relation also suggests a relatively weak dependence on skin-friction shear velocity (E si ∝ u 1.77 * skin ). Similar to Garcia and Parker (1991) and Wright and Parker (2004b), we modified Eq. (25) such that the predicted entrainment rate is limited at 0.3, as total suspended sediment concentrations greater than 30 % by volume are not physically reasonable for dilute, turbulent flows. In addition, a threshold (0.015), best fit by eye, was added to the entrainment relation because the concentration data fall below the trend of the regression relation at the lower flow strengths (Fig. 8b), suggesting a threshold of significant sediment entrainment at the reference level. The resulting equation has the following form (Fig. 8b): The equivalent formula for the best-fit two-parameter model without a form drag correction (Table 5) is . Best-fit models for the dimensionless entrainment rate, E si , for grain-size-specific data with one parameter (a), two parameters (b), and three parameters (c); equations are in bold in Table 5. The plotting is done such that the exponent on the last parameter is unity to allow for comparison to plots of Garcia and Parker (1991) and Wright and Parker (2004b). Yellow filled symbols are the geometric mean, and the error bars indicate the geometric standard deviation within each bin. Table 5. One-, two-, and three-parameter models for the entrainment parameter of the form E si = AP e 1 1 P e 2 2 P e 3 3 . The two best-fitting versions of each model are shown, one with u * and one with u * skin , and the bold value is the best fitting of the two (highest r 2 ). See Table S2 9. Box plot of predicted versus measured entrainment parameters, E si , for our models and previous models. The best-fit one-, two-, and three-parameter models correspond to the equations in bold in Table 5.

Predicting sediment concentration
In Sect. 3.1 and 3.2 we found the best-fit models for the dimensionless entrainment rate (E si ) and Rouse number (P i ). However, ultimately, we want to predict the sediment concentration throughout the water column, and the best-fit models for (E si ) and (P i ) do not necessarily combine to yield the best-fit model for sediment concentration owing to nonlinearity in E si , P i , and the Rouse equation (Eq. 1). Here we used different combinations of our preferred one-, two-, and three-parameter models for entrainment (E si ) and the Rouse number (P i ) to predict the grainsize-specific concentrations at each data point in the water column for all our entries in the database and assessed model performance (Table 6). The concentration predictions improve as more parameters are added to the entrainment model, whereas a Rouse model with more than one parameter makes the predictions worse (Table 6). Our preferred model for concentration throughout the water depth uses the two-parameter model to predict the entrainment rate E si = 4.74 × 10 −4 u * skin w si 1.77 Fr 1.18 and the one- Figure 10. Entrainment parameter, E si , normalized to isolate the effect of (a) grain size and (b) skin-friction shear velocity to check for spurious correlation. Yellow filled symbols are the geometric mean, and the error bars indicate the geometric standard deviation within each bin. The best-fit models are Eqs. (24) and (25). parameter model for the Rouse number P i = (u * /w si ) −0.45 (Fig. 11). Such a model gives significantly better predictions than the most basic formulation that uses one-parameter models for entrainment and the Rouse number (r 2 of 0.87 versus 0.65; Table 6). The goodness of fit is also fairly constant over the height of the flow; the upper 1/3 of the flow has a slightly lower r 2 (0.80) than the predictions from the lowest 1/3 of the flow (r 2 = 0.89) (Fig. 11); 95 % of the data are predicted within a factor of 9.

Extension to gravel
The suspended sediment data that we used to calibrate the entrainment relation cover material in the sand range. To evaluate how our entrainment relation performs for coarser suspended sediment, we used the empirical saltation equations for gravel to infer bedload-layer concentrations, C b , and interpolated these to the reference level (0.1H ) to infer C a for gravel (Sect. 2.3). Importantly, the gravel concentrations at the reference height can be predicted from the saltation model (Sect. 2.3) using only the independent parameters of Fr and τ * , similar to our best-fitting two-parameter entrainment model. To compare to gravel, we assumed a uniform sediment size and used the version of our best-fit twoparameter model that does not include a form drag correction (Eq. 26b). This was done because gravel saltation studies were typically performed over a planar bed, whereas form drag likely played a significant role in the sand datasets owing to dunes. The marked parameter space in Fig. 12 shows the expected range of C b and C a for a wide range of model input parameters for gravel: 0.1 < Fr < 1 and 1 < τ * < 1000. Predicted concentrations in the gravel bedload layer are up to several orders of magnitude higher than the predictions from our entrainment relation at 10 % of the flow depth (Fig. 12). However, due to the rapid decrease in sediment concentration away from the bed predicted by the Rouse profile, the concentration inferred at 0.1H for gravel overlaps the empirical relation for sand, implying that Eq. (27) might also be a good predictor of gravel entrainment.

Physical rationale for model dependencies
The explanatory variables in the models for P i and β i are u * skin /w si and C f . The dependency on these parameters Table 6. Coefficient of determination (r 2 ) of grain-size-specific sediment concentration at each sample elevation in the water column for different combinations of models for the entrainment, E si , and Rouse, P i , parameters. The favored model combination is highlighted in bold font.  Figure 12. Best-fit two-parameter model for sediment entrainment without form drag, Eq. (26b), presented here for uniform gravel. The plotting is done such that the exponent on the last parameter is unity to allow for comparison to plots of Garcia and Parker (1991) and Wright and Parker (2004b). The colored regions are prediction envelopes for synthetic gravel data using the bedload-layer equations to estimate the bedload volumetric concentration, C b , and the Rouse profile to interpolate higher in the water column to the reference level at z = 0.1H (C a ). The synthetic data envelope represents a wide range of parameter space: 0.1 < Fr < 1 and 1 < τ * < 1000. could reflect processes that result in the turbulent diffusivity for suspended sediment being different from the parabolic relation used in the derivation of the Rouse profile (e.g., Bennett et al., 1998). Compared to standard Rouse theory, in which the Rouse parameter is inversely dependent on u * skin /w si (or u * /w si ), our results indicate a significant non-linearity. Our one-parameter model indicates that concentration profiles are better mixed than standard theory for small u * /w si and are more stratified for larger u * /w si , with a transition point at about u * /w si = 7, corresponding to β i = 1. Sediment-induced stratification is often cited (Winterwerp, 2006;Wright and Parker, 2004b, a) as a factor that decreases the mixing of sediment (i.e., β < 1). This effect is particularly important when the absolute concentration is high and may help explain why our best-fit model is more stratified (i.e., β < 1) than standard Rouse theory for large u * /w si . An alternative hypothesis by Nielsen and Teakle (2004) is that for steep concentration gradients, the size of the turbulent eddies is large relative to the mean height of sediment in the flow. Under these circumstances, large eddies might more effectively mix sediment that is concentrated close to the bed. This process may explain the better mixed concentration profiles we observed (i.e., β i > 1) compared to standard theory at small u * /w si when near-bed concentration gradients are large.
Our relation implies that β i correlates positively with the flow friction coefficient C f . This dependency could be from bedforms; bedforms increase C f due to form drag (Engelund and Hansen, 1967), and they may also increase the vertical mixing of sediment by deflecting transport paths up the stoss side of dunes and mixing suspended sediment in the turbulent wake in the lee of dunes. Similarly, Santini et al. (2019) found that β correlated positively with H /D (Fig. 2c), which is another measure for bed roughness in flows without bedforms. In agreement, Graf and Cellino (2002) reviewed a number of experimental studies and found that β < 1 for ex-periments without bedforms and β > 1 for experiments with bedforms. Flow resistance (C f ) can also be smaller in flows with sediment-induced stratification, which correlates with smaller β i (e.g., Wright and Parker, 2004a).
Our new empirical relation for sediment entrainment suggests that only u * skin /w si (or u * /w si ) and Fr are needed to predict entrainment rate, similar to the forward model we developed for the suspension of gravel (Sect. 2.3). The ratio u * skin /w si describes the fluid forces relative to gravitational settling; similar parameters have appeared in all previous relations that we reviewed ( Table 2). The reason for the increase in entrainment with the Froude number is less clear. A small Froude number implies a deep and low-gradient flow, and Froude numbers are typically smaller in natural rivers compared to flume experiments. Wright and Parker (2004b) introduced an entrainment relation with a bed slope dependency and argued that entrainment in large low-sloping rivers is reduced due to stratification effects. Sediment-induced stratification causing the damping of turbulence might be the cause of the Fr dependency in our relation. Regardless, we found a better fit with the data using Fr than using C a S or S (Table S2), the parameters suggested by Wright and Parker (2004b). The Froude number might also influence the size and shape of bedforms (Vanoni, 1974), which can affect boundary layer dynamics and near-bed turbulence, as well as surface waves. The Froude number is the ratio between inertial to gravitational forces or, more formally, the ratio of a unidirectional flow velocity to the celerity of a shallow water wave (i.e., Fr = U/ √ gH ). However, we do not think these are the reasons why the dependence of E si on Fr emerges in our analysis; rather, our forward model (Sect. 2.3) suggests that it emerges as a controlling parameter because of the role of U and H in determining bedload-layer concentrations, which is discussed in more detail in Sect. 4.3.
Many of the dimensionless parameters we evaluated correlate with each other in rivers. While u * skin /w si or u * /w si was consistently the dominant variable, several of the possible secondary variables had similar explanatory power as the ones given in our preferred models. For sediment entrainment, some of the more highly ranked variables are C a /S and τ * skin in one-parameter models and H /D 50 , S and D i /D 50 in two-parameter models; for the Rouse number they include R p,skin and τ * for one-parameter models and Re p and C a /S in two-parameter models (Table S2). In addition, there is some systematic deviation between datasets for different rivers (Fig. 8), and different parameters and exponents might better minimize residuals at specific locations. While our model has been fit to a large dataset, it has not been validated with independent data. More data are clearly needed on grain-size-specific concentration profiles for equilibrium suspensions under a wide range of u * /w si , Fr, and particle sizes, including coarse sand and gravel. Figure 13. Grain-size-specific concentration, C i , from the lowest 10 % of the flow (i.e., below the reference level z = 0.1H ) relative to the concentration at the reference level, C ai , as a function of depth above the bed, z, relative to the flow depth, H . Yellow filled symbols are the geometric mean; the error bars indicate the geometric standard deviation within each elevation bin, and the black line is the geometric mean of all data representing C i = 1.66C ai .

Predicting sediment concentration below the reference level
Combining the predictions for the grain-size-specific reference concentration and Rouse parameter allows for the calculation of the grain-size-specific sediment concentration throughout the water column (Fig. 11). Our relation shows that the grain-size-specific sediment concentration at any given elevation can be predicted to relatively high accuracy (r 2 = 0.87) using the preferred combination of a twoparameter entrainment relation and the one-parameter Rouse number relation (Sect. 3.3). However, it is unclear how to evaluate the sediment concentration below the reference level, which could constitute a significant portion of the sediment load. Bedload fluxes are notoriously difficult to estimate; further, it is unclear if the region below the reference level is entirely bedload or if it is bedload and suspended load (e.g., Fig. 1b). Ashley et al. (2020) showed how bedload sediment fluxes can be estimated from discharge-averaged suspended sand concentrations, but their method does not predict the concentration profile. One approach might be to assume that the sediment concentration is uniform below the reference level at C i = C ai , as might be the case in a wellmixed bedload layer (e.g., McLean, 1992) (Fig. 1b). However, this assumption is inconsistent with our analysis of the saltation equations, which shows that bedload layers should have greater concentrations and be less than 10 % of the flow depth (Fig. 12), especially for sand in deep flows. A second approach is to use the Rouse profile to extrapolate towards the bed or towards the top of the bedload layer; however, this approach is also problematic because the Rouse profile predicts infinitely large concentrations at the bed. Some of our datasets had concentration measurements below the reference level (Fig. 13). The data are scattered, but the binned data points suggest a nearly constant concentration with elevation below the reference level. For lack of a better method, we propose using the geometric mean of these data as C i = 1.66C ai to approximate the concentration-depth profile below the reference level (Fig. 13).

Extension to gravel and bedload-layer theory
The transport of sand and gravel is often modeled using different empirical formulas, which hinders the modeling of systems of mixed gravel-sand transport (Wilcock and Crowe, 2003) and gravel-sand transitions (Paola et al., 1992;Lamb and Venditti, 2016). Gravel can also be in suspension in bedrock and steep mountain rivers, as well as during large floods (Hartshorn et al., 2002;Larsen and Lamb, 2016), and it would be useful to have an entrainment relation to model sediment transport and bedrock erosion in these settings (Lamb et al., 2008a;Scheingross et al., 2014). Unfortunately, sand and gravel are often treated separately, which limits our ability to develop a unified sediment transport theory. Our entrainment relation for sand matches expectations from gravel saltation models (Fig. 12), suggesting that the entrainment relation may be used for sand, gravel, or mixtures of the two. However, we currently lack data on gravel suspension profiles, and developing methods to acquire such data should be the focus of future efforts. It is also unclear how entrainment is complicated by bimodal mixtures (Wilcock and Crowe, 2003) and complex flow hydraulics in steep rivers with large roughness that can significantly affect lift forces and near-bed turbulence (Lamb et al., 2017).
The good fit between the modeled gravel concentrations and the measured sand data for near-bed sediment concentration suggests that the bedload-layer equations (Sect. 2.3) might also be used to generate a forward model for nearbed sediment concentration that works for sand and gravel systems, similarly to previous efforts (e.g., McLean, 1992). To evaluate this possibility, we used the saltation equations (Eqs. 17-21) to calculate sediment concentrations within the bedload layer for conditions corresponding to our dataset entries for sand-bedded rivers. To extend Eq. (18) to grain size mixtures, we let D = D i and τ * c = τ * ci and used a hiding function (Parker et al., 1982) so that where τ * c50 is the critical Shields number for the mediansized bed sediment. Equation (28) accounts for the hiding of smaller grains between larger grains, which renders the smaller grains less mobile, and the exposure of larger grains into the flow, which renders them more mobile. For γ = 1, all grains in a bed mixture move at the same bed stress, while for γ = 0, the critical stress for motion is proportional to grain weight. Gravel-bedded rivers typically have γ = 0.9 (Parker, 1990). For sand, we evaluated the critical Shields number following Brownlie (1981), rather than Eq. (19), which is intended for gravel only. We then calculated the concentrations at 10 % of the flow depth using the Rouse equation (Eq. 1) with the best-fit oneparameter model for the Rouse number. Although the saltation equations in Sect. 2.3 were calibrated for gravel, similar relations have also been used for sand (e.g., Lamb et al., 2008a). Surprisingly, the near-bed sand concentrations for the entries in our database (Sect. 2.1) are relatively accurately predicted by the bedload forward model without any parameter fitting (Fig. 14b). In fact, the predictions by the bedload-layer forward model are only slightly worse than the predictions by our preferred two-parameter entrainment model (Fig. 14a) that was fit to the data (r 2 of 0.68 vs. 0.87). Importantly, the forward model yields the same controlling parameters as the empirical model, namely the Shields number, the Froude number, and the bed grain size distribution. The forward model suggests that the dependence of E si on u * skin /w si (or u * /w si ) is predominantly because larger u * skin /w si correlates with larger τ * and larger bedload-layer concentrations (Eqs. 17 and 18), as well as more efficient mixing of the bedload-layer sediment up to the reference height (Eq. 1). Less intuitive is the dependence of E si on the Froude number. In the forward model, this dependence is because larger Fr for the same τ * correlates with larger D i /H (which can be shown by manipulating Eq. 7). In turn, larger D i /H with constant τ * correlates with greater bedload fluxes (q b ∝ D 3/2 i in Eq. 18), smaller bedload-layer heights (due to smaller H in Eq. 21), and slower bedload-layer velocities (due to smaller U in Eq. 20), all of which increase the sediment concentration in the bedload layer (Eq. 17). Thus, the forward model indicates that the Fr dependency on E si emerges because bedload-layer dynamics depend on U and H , and explanations for this dependency that rely on stratification, bedforms, or surface waves are not necessary. The bedload-layer model, while slightly more complicated to implement, may provide a more robust solution when working outside the parameter space used to derive the empirical model, since it has a more physical basis. For example, the forward model can explicitly account for gravity and other physical properties of the sediment and fluid. In addition, the bedload-layer model allows for a more mechanistic link between the bedload and suspended load and avoids uncertainty in how to evaluate the sediment concentration below the reference level. The model might likely be improved by accounting for bedform-induced form drag, especially for sand-bedded rivers. That said, more accurate predictions are still achieved with the empirical entrainment relation. Figure 14. Measured versus predicted grain-size-specific near-bed concentration, C ai , at the reference level z = 0.1H for (a) our twoparameter empirical relation (Eq. 25) and (b) for the forward modeling approach using the bedload-layer equations (Sect. 4.3). Data (blue markers) are from the entire database of suspended sand described in Sect. 2.1 (Table S1) at z = 0.1H . Yellow symbols are the geometric mean, and the error bars indicate the geometric standard deviation within each bin.

Conclusions
We proposed new empirical models for the entrainment of bed material into suspension and for the shape of the concentration profile governed by the Rouse parameter. The models were obtained by regression against suspended sediment data from eight different rivers and six experimental studies. The data cover a wide range of bed material grain sizes (44-517 µm) and flow depths (0.06-32 m) and include grainsize-specific data with up to 60 size classes. Our analysis of these data suggests that the near-bed sediment concentration increases with the ratio between shear velocity or settling velocity (u * skin /w s or u * /w si ) and the Froude number -both parameters also emerge as the key controlling variables in a forward model based on bedload-layer concentrations. A parameter such as u * /w si , which represents the ratio of fluid force to particle settling, was also present in previous relations. The Froude number dependence is less clear; it could be due to stronger sediment-induced stratification in large low-Fr rivers, but our forward model indicates that it emerges because flow velocity and flow depth impact bedload-layer concentrations. Our preferred Rouse parameter model for the shape of the concentration profile suggests that the sediment concentration is better mixed in the water column with larger u * skin /w si and a larger bed friction coefficient, C f . The Rouse number is not inversely proportional to u * /w si , unlike standard Rouse theory, indicating that sediment is more stratified than expected with u * /w si >∼ 7 and better mixed than expected with u * /w si <∼ 7, possibly due to the competing effects of sediment-induced stratification when absolute concentrations are large (corresponding to large u * /w si ) and enhanced turbulent mixing when concentration gradients are steep (corresponding to small u * /w si ). The dependence of the Rouse number on the bed friction coefficient might result from increased turbulence close to the bed in rivers with large bed roughness or bedforms. We also demonstrated that nearbed concentrations can be accurately predicted with saltation equations that have been tested previously for gravel, suggesting a unified framework to model sand and gravel transport in rivers.
Data availability. All data used in the analysis are provided in Table S1 of the Supplement. Model results for all possible variable combinations are given in Table S2. Table S1 contains all the suspended sediment data that were used to find the empirical relations. Table S2 contains all entrainment and Rouse number models ranked according to goodness of fit as indicated by r 2 . The supplement related to this article is available online at: https://doi.org/10.5194/esurf-8-485-2020supplement.

Supplement.
Author contributions. MPL and GP conceived the study. JdL compiled data and led data analysis with input from MPL. MPL and JdL wrote the initial paper. AJM, DH, JGV, and JAN supplied suspended sediment data and contributed to the final paper.
Financial support. This research has been supported by the US National Science Foundation (grant no. EAR 1427262), the Caltech Discovery Fund, and the NSF Graduate Research Fellowship (grant no. EAR 1842494).
Review statement. This paper was edited by Eric Lajeunesse and reviewed by two anonymous referees.