A bed load transport equation based on the spatial distribution of shear stress – Oak Creek revisited

Bed load transport formulations for gravel-bed rivers are often based on reach-averaged shear stress values. However, the complexity of the flow field in these systems results in wide distributions of shear stress, whose effects on bed load transport are not well captured by the frequently used equations, leading to inaccurate estimates of sediment transport. Here, we modified a subsurface-based bed load transport equation to include the complete distributions of shear stress generated by a given flow within a reach. The equation was calibrated and verified using bed load data measured at Oak Creek, OR. The spatially variable flow field characterization was obtained using a two-dimensional flow model calibrated over a wide range of flows between 0.1 and 1.0 of bankfull discharge. The shape of the distributions of shear stress was remarkably similar across different discharge levels, which allowed it to be parameterized in terms of discharge using a gamma function. When discharge is high enough to mobilize the pavement layer (1.0 m³ s⁻¹ in Oak Creek), the proposed transport equation had a similar performance to the original formulation based on reach-averaged shear stress values. In addition, the proposed equation predicts bed load transport rates for lower flows when the pavement layer is still present because it accounts for bed load transport occurring in a small fraction of the channel bed that experiences high values of shear stress. This is an improvement over the original equation, which fails to estimate this bed load flux by relying solely on reach-average shear stress values.

1 Introduction

Predicting bed load is both expensive and practically challenging, as data from a wide range of flows are required to develop robust relationships between discharge and load. In addition, characterizing bed load at high flow levels – that transport the majority of the sediment – is often dangerous (Bunte et al., 2008). Samples collected using handheld devices can be widely variable due to factors related to variations in their orifice size and the sampling time (Beschta, 1981; Emmett, 1980; Pitlick, 1988; Vericat et al., 2006). While advances in safe, accurate sediment sampling technology such as bed load traps (Bunte et al., 2008), radio tracers (Bradley and Tucker, 2012; May and Pryor, 2014; Olinde and Johnson, 2015; Schmidt and Ergenzinger, 1992), and acoustic impact methods (Rickenmann and McArdell, 2007; Turowski and Rickenmann, 2011; Wyss et al., 2016a–c; Yager et al., 2012b) provide possible alternatives to handheld samplers, field efforts remain expensive and out of reach for many practical applications.

Bed load modeling can be a convenient alternative to measuring bed load in the field. The development of empirical bed load relationships has progressed significantly over the past 3 decades such that many formulations allow for the estimation of bed load based on hydraulic and grain size information. In general sediment transport equations are based on reach-averaged 1D shear stress estimates and the surface (e.g., Barry et al., 2004; Parker, 1990; Recking, 2013b; Wilcock and Crowe, 2003) or subsurface (Parker et al., 1982; Parker and Klingeman, 1982) grain size information.

Many of these sediment transport equations (e.g., Parker et al., 1982; Parker and Klingeman, 1982) were developed based on data from Oak Creek, a steep, coarse, gravel-bed stream in the Oregon Coast Range (Milhous, 1973). The Oak Creek dataset was collected using a vortex sampler between 1969 and 1990; data from 1971 were published in the thesis work of Milhous (1973). The dataset is unique because the vortex sampling method enabled capturing of the entire bed load flux of sand–cobble size particles for a wide range of flows over long time periods, reducing the error associated with handheld samplers (Parker et al., 1982). Although it has been reported that the efficiency of the vortex sampler decreased for smaller grain sizes (Milhous, 1973; O'leary and Beschta, 1981), the Oak Creek dataset remains one of the most comprehensive to date. The Oak Creek-based transport equations were developed by collapsing the relations between reference conditions for the motion of different grain sizes into single functions (i.e., a similarity collapse) (Einstein, 1950; Parker, 1990; Parker et al., 1982; Parker and Klingeman, 1982). Both Parker and Klingeman (1982) and Parker et al. (1982) limited their analysis to flows during which the surface channel layer was mobilized (“pavement” was broken) to develop their transport functions. Parker et al. (1982) computes total bed load (Q_b) based on a single grain size (the median – D₅₀), whereas Parker and Klingeman (1982) expands that relationship to the entire grain size distribution (GSD). This is accomplished by introducing a hiding function that accounts for differences in the hiding and exposure of particles to the flow in mixed-sized beds. Additionally, Parker and Klingeman (1982) incorporated a low-flow transport relation to estimate the GSD of Q_b at a full range of flows. Later, Parker (1990) modified Parker and Klingeman (1982) equations to be based in the surface grain size distribution.

Although the transport relations of Parker et al. (1982) and Parker and Klingeman (1982) have been successfully applied to many rivers, the work of Recking (2013b) highlighted the variability that can be incorporated into Q_b estimates due to uncertainty in input shear stress (τ) values. The high spatial variability in τ throughout a river reach has been well documented (Clayton and Pitlick, 2007; Katz et al., 2018; Lisle et al., 2000; May et al., 2009; McDonald et al., 2010; Monsalve et al., 2016; Recking, 2013a; Segura and Pitlick, 2015; Yager et al., 2018). However, most transport functions, including Parker and Klingeman (1982) and Parker et al. (1982), utilized reach-averaged estimates of τ in their calculations and are highly sensitive to uncertainties in these values due to the nonlinear exponents on each function (Recking, 2013a). Significant differences in bed load estimates computed using τ from one-dimensional (1D) and two-dimensional (2D) approximations have been found because of the spatial variability in τ (Ferguson, 2003; Gomez and Church, 1989; Recking, 2013a). Thus, the simplification of τ to a 1D variable may not capture spatial changes in bed load associated with localized values of high τ (Segura and Pitlick, 2015).

The objective of this study is to develop a bed load transport equation based on the subsurface GSD that uses the complete shear stress distribution for different discharge levels within a specific reach. Although a surface-based equation (e.g., Parker, 1990) would have been stronger at estimating bed load transport at low flows, we adopted a subsurface-based equation (Parker and Klingeman, 1982) because it has fewer parameters to adjust. A subsurface-based equation also allows the consideration of sand sizes, which are commonly found in the bed load (Clayton and Pitlick, 2007; Hassan and Church, 2001; Lisle, 1995; Mueller et al., 2005; Recking, 2010; Segura and Pitlick, 2015). This new approach is developed using field measurements of bed load transport rates and GSD, river topography, and 2D flow modeling. The performance of the new equation is then tested using the historic Oak Creek dataset (Milhous, 1973). Specific objectives of our study are as follows:

• i.

  to analyze the characteristics of shear stress distributions over a wide range of discharge levels

• ii.

  to generate synthetic shear stress distributions based solely on discharge

• iii.

  to modify a reach-averaged subsurface-based equation (Parker and Klingeman, 1982) developed for Oak Creek to use complete shear stress distributions

• iv.

  to test the performance of the proposed equation for a wide range of discharge levels.

Figure 1 Location of the study reach in Oak Creek, Oregon (44^∘23^′19.092^′′ N, 123^∘19^′51.312^′′ W). Contours every 0.1 m are indicated.

Figure 2 Surface and subsurface grain size distribution (GSD). The average surface GSD is based on 23 cross sections (XS) and average subsurface GSD based on two samples of the substrate collected from exposed bars.

2 Methods

2.1 Study area

This study was conducted in Oak Creek, a cobble–gravel stream located in the Oregon Coast Range (Milhous, 1973, Figs. 1 and 2). The catchment drains 7 km² of forest land underlain by basaltic lithology (Milhous, 1973; O'Connor et al., 2014). The climate is Mediterranean with wet winters and cool/mild summers. Elevations within the Oak Creek watershed range from 143 to 664 m (Paustian and Beschta, 1979). The basin is located in the McDonald-Dunn Forest, which is owned and managed by the College of Forestry at Oregon State University and dominated by Douglas fir (Pseudotsuga menziesii) and Oregon white oak (Quercus sp.). In the riparian area, vegetation is dominated by alder (Alnus sp.), black cottonwood (Populus trichocarpa), and big leaf maple (Acer macrophyllum), with lower densities of Douglas fir and Oregon white oak. The 150 m study reach has a pool–riffle sequence in the upstream end and a relatively straight section in the downstream section (Katz et al., 2018) and is located directly upstream from a historic sediment-transport sampling facility where bed load samples were collected between 1969 and 1973 (Milhous, 1973). The site has a rectangular cement weir in which a stage–discharge (Q) relationship has been developed (Katz et al., 2018). The stream has a slope (S_b) of 0.014 m m⁻¹ and bankfull dimensions of 6 m in width and 0.46 m in depth. Recent field observations indicate that bankfull discharge (Q_bf) is 3.4 m³ s⁻¹ (Katz et al., 2018), which is similar to the bankfull discharge reported almost 40 years ago by Milhous (1973). The stream bed is armored with coarser surface overlying a finer subsurface. The surface D₅₀ is 45 mm, while the subsurface D_{50_s} is 21 mm (Katz et al., 2018) (Fig. 2).

2.2 Two-dimensional modeling

Spatial distributions of the flow field, in particular local shear stresses, were estimated for seven discharges (0.4 to 3.4 m³ s⁻¹, equivalent to 0.12 Q_bf to Q_bf) using the Flow and Sediment Transport with Morphological Evolution of Channels (FaSTMECH) 2D flow solver (McDonald et al., 2010). Specific details of the modeling effort can be found in Katz et al. (2018). The model has also been described and used in several studies (e.g., Clayton and Pitlick, 2007; Conner and Tonina, 2014; Kinzel et al., 2009; Lisle et al., 2000; Maturana et al., 2014; Mcdonald et al., 2005; Monsalve et al., 2016; Mueller and Pitlick, 2014; Nelson and McDonald, 1995; Nelson and Smith, 1989; Nelson et al., 2010; Segura and Pitlick, 2015) – therefore, only the most relevant characteristics of it are described here. The model uses a finite-difference solution to the vertically integrated conservation of mass and momentum equations (Nelson et al., 2003) with calculations performed in an orthogonal curvilinear grid that follows the surveyed planform topography of the channel (Nelson and Smith, 1989). Roughness is included using a unitless drag coefficient (C_d). A zero-equation model for the lateral eddy viscosity (LEV) that assumes homogeneous and isotropic turbulence is used for turbulence closure (Barton et al., 2005; Miller and Cluer, 1998; Nelson et al., 2003). For our models, C_d ranged from 0.017 to 0.04, and LEV ranged 0.0010 to 0.0032 (Katz et al., 2018). The calibration indicated strong model fits in terms of water surface elevation, with root mean square errors (RMSEs) between 0.025 and 0.048 m and R²>0.99 (Katz et al., 2018).

The local shear stress (τ_xy) was calculated at every grid node in the model domain as a function of C_d, the vertically averaged streamwise (u) and cross-stream (v) velocities, and water density (ρ), assumed as 1000 kg m⁻³.

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{xy}=\mathit{\rho }{C}_{\mathrm{d}}\left(u_{xy}^{\mathrm{2}}+v_{xy}^{\mathrm{2}}\right),\end{array}$$

where the subscripts x and y correspond to the stream-wise and cross-stream directions.

2.3 Shear stress distribution analysis

Characteristics of the distributions of predicted τ_xy were analyzed as a function of discharge. We produced histograms of the mean-normalized shear stress distribution (τ∕〈τ〉) (subscripts x and y were dropped for simplicity) to compare patterns between flows. For each flow level we fitted the frequency distributions of τ∕〈τ〉 to a two-parameter gamma function (Nicholas, 2000; Paola, 1996; Pitlick et al., 2012; Recking, 2013a; Segura and Pitlick, 2015):

$$\begin{array}{}\text{(2)}& f\left(\mathit{\tau }\right)={\displaystyle \frac{{\mathit{\alpha }}^{\mathit{\alpha }}(\mathit{\tau }/\langle \mathit{\tau }\rangle {)}^{(\mathit{\alpha }-\mathrm{1})}{e}^{-\mathit{\alpha }(\mathit{\tau }/\langle \mathit{\tau }\rangle )}}{\langle \mathit{\tau }\rangle \mathrm{\Gamma }\left(\mathit{\alpha }\right)}},\end{array}$$

where Γ is the standard gamma function, α is the shape parameter, and $\mathit{\beta }=\langle \mathit{\tau }\rangle /\mathit{\alpha }$ is the scale parameter. The parameters of the gamma function that best fitted the distributions were found using the maximum likelihood estimation (MLE) method (Bevington and Robinson, 2003). We assessed the goodness of fit of the gamma function in each flow event by computing the room mean square error (RMSE) and the reduced chi square (${\mathit{\chi }}_v^{\mathrm{2}}$), defined as chi square (χ²) divided by the number of degrees of freedom, according to

$$\begin{array}{}\text{(3)}& {\mathit{\chi }}^{\mathrm{2}}=\sum {\displaystyle \frac{{\left[f_k-f\left(x_k\right)\right]}^{\mathrm{2}}}{{\mathit{\sigma }}_k^{\mathrm{2}}}},\end{array}$$

where and f_k and f(x_k) are observed and predicted mean-normalized shear stress frequencies in a given bin interval, k. The uncertainty associated with the observed frequencies, ${\mathit{\sigma }}_k^{\mathrm{2}}$, was estimated as the square of the number of observations in each bin (Bevington and Robinson, 2003; Press et al., 2007). Initially, in all cases, we specified the bin width using the Freedman–Diaconis rule (Freedman and Diaconis, 1981). To improve statistics when the number of τ values in a given bin was less than five, we joined two consecutive bins until all bins had five or more τ values. Typically, for the used goodness of fit indicators, an excellent fit is ${\mathit{\chi }}_v^{\mathrm{2}}\le \mathrm{1}$ and a RMSE of zero (Bevington and Robinson, 2003; Press et al., 2007).

2.4 Sediment transport equations

The original subsurface-based sediment transport equation of Parker and Klingeman (1982) was modified to explicitly consider the spatial distribution of shear stress. This equation was chosen because it was developed from measurements collected in the same reach as this study, it gives accurate estimates of bed load transport, and it is relatively simple to extend for our purposes (see below). The modified version of the Parker and Klingeman (1982) equation was formulated such that it accounts for the bed load transported by each increment of shear stress, which means that it considered the range of local contributions of τ across the channel bed. By doing so, all τ values, even those less-frequent, high-magnitude shear stresses, are explicitly included in the calculations. To obtain the new equation, the parameters of the Einstein bed load function (G) proposed by Parker (1978) were relaxed and fitted as new parameters. The parameter values were optimized based on the fit of volumetric transport rate per unit width of channel (q) and the bed load GSD. Like the original equation, we only consider discharges of approximately 1 m³ s⁻¹ or higher to calibrate the new equations (for calibration purposes, our lower discharge was 0.99 m³ s⁻¹). The fitting procedure of the parameters minimized the absolute error between predicted and measured q and maximized the Nash–Sutcliffe efficiency index (Nash and Sutcliffe, 1970) using the calculated and observed bed load GSD. Equal importance (equal weight) was given to the fit of q and to the fit of the bed load GSD.

The new equation was based on the locally dimensionless shear stress (τ^*):

$$\begin{array}{}\text{(4)}& {\mathit{\tau }}^{*}={\displaystyle \frac{\mathit{\tau }}{\left({\mathit{\rho }}_{\mathrm{s}}-\mathit{\rho }\right)gD}},\end{array}$$

where ρ_s is the sediment density, g is the acceleration due to gravity, and D is the grain size. Notice that for a given flow discharge τ^* has a distribution of values depending on the local τ (previously defined as τ_xy) and variations in the fraction of the GSD. The original transport relation of Parker and Klingeman (1982) (Eq. 5) is valid for uniform grain sizes and φ>1, with φ being the transport stage (Eq. 6),

$$\begin{array}{}\text{(5)}& G={\displaystyle \frac{W^{*}}{W_r^{*}}}=\mathrm{5.6}\times {\mathrm{10}}^{-\mathrm{3}}{\left(\mathrm{1}-{\displaystyle \frac{\mathrm{0.853}}{\mathit{\varphi }}}\right)}^{\mathrm{4.5}},\end{array}$$

where the subscript r denotes a reference value associated with a small but measurable transport rate. Transport stage (ϕ) is defined as

$$\begin{array}{}\text{(6)}& \mathit{\varphi }={\mathit{\tau }}^{*}/{\mathit{\tau }}_{\mathrm{r}}^{*}.\end{array}$$

The dimensionless transport rate, W^* (Eq. 5), is defined as

$$\begin{array}{}\text{(7)}& W^{*}={\displaystyle \frac{(s-\mathrm{1})gq}{(\mathit{\tau }/\mathit{\rho }{)}^{\mathrm{1.5}}}},\end{array}$$

where s is the specific gravity of sediment ($s={\mathit{\rho }}_{\mathrm{s}}/\mathit{\rho }$).

Figure 3 Frequency distributions of mean-normalized shear stress (τ∕〈τ〉) for the seven discharge levels. Fitted gamma distribution curves are shown as dashed lines. Discharges values are indicated in the upper-right corner of each panel.

We extended Eq. (5) to include all grain size fractions in the subsurface GSD (D_i; subscript i denotes the size range) and ϕ_i>0.95. In the most general form, the equation for the dimensionless transport rate is $W_i^{*}=\mathrm{0.0025}{G}_i^{\prime }$, where the constant is the reference transport rate of Parker and Klingeman (1982) ($W_{\mathrm{r}}^{*}=\mathrm{0.0025}$), and $G_i^{\prime }$ is the new (modified) transport relation. The proposed relation is a two-part equation applicable to sediment mixtures:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{W}_i^{*}=& \mathrm{0.0025}\times {\mathrm{10}}^{-\mathrm{3}}\exp \left(\mathrm{26.6}\left({\mathit{\varphi }}_i-\mathrm{1}\right)-\mathrm{19.53}{\left({\mathit{\varphi }}_i-\mathrm{1}\right)}^{\mathrm{2}}\right),\\ & \mathrm{for}\mathrm{0.95}<{\mathit{\varphi }}_i<\mathrm{1.65}\\ W_i^{*}=& \mathrm{0.57}{\left(\mathrm{1}-{\displaystyle \frac{\mathrm{0.853}}{{\mathit{\varphi }}_i}}\right)}^{\mathrm{4.5}},\\ & \mathrm{for}{\mathit{\varphi }}_i\ge \mathrm{1.65}.\end{array}\end{array}$$

To account for the mobility of individual grain sizes, we used the Parker and Klingeman (1982) hiding function:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{{\mathit{\tau }}_{\mathrm{r}i}^{*}}{{\mathit{\tau }}_{\mathrm{r}\mathrm{50}}^{*}}}={\left({\displaystyle \frac{D_i}{D_{\mathrm{50}}}}\right)}^{-\mathrm{0.982}},\end{array}$$

where ${\mathit{\tau }}_{\mathrm{r}\mathrm{50}}^{*}=\mathrm{0.0876}$ is the reference Shields stress for the median grain size of the subsurface Parker and Klingeman (1982) obtained for the same reach. The transport stage (Eq. 6), valid for any grain size D_i and for the entire distribution of shear stress values (i.e., every ${\mathit{\tau }}_i^{*}$), was rewritten as

$$\begin{array}{}\text{(10)}& {\mathit{\varphi }}_i={\mathit{\tau }}_i^{*}/{\mathit{\tau }}_{\mathrm{r}i}^{*}.\end{array}$$

To obtain the volumetric transport rate, the predicted shear stresses were grouped in a series of intervals (τ_j; subscript j denotes an interval of τ values) with a regular shear stress increment (Δτ_j=0.25 N m⁻²). For all discharges, τ_j was defined such that it ranges from zero to the maximum predicted shear stress value. For a given D_i and τ_j, the volumetric transport rate per unit width (q_ij) is

$$\begin{array}{}\text{(11)}& q_{ij}={\displaystyle \frac{{\left(\frac{{\mathit{\tau }}_j}{\mathit{\rho }}\right)}^{\mathrm{1.5}}{F}_i{W}_{ij}}{(s-\mathrm{1})g}}{f}_{{\mathit{\tau }}_j},\end{array}$$

where F_i is the volume fraction of the ith grain-size class in the subsurface GSD, W_ij is calculated using Eq. (8) for each τ_j, and $f_{{\mathit{\tau }}_j}$ is the fraction of the bed area where a certain τ_j acts. The width-integrated volumetric transport rate for a given flow event is

$$\begin{array}{}\text{(12)}& q_{\mathrm{b}}=b\sum _i\sum _j{q}_{ij},\end{array}$$

with b being the width of the gravel bed. In all bed load estimations, sand grains likely to move in suspension were excluded; thus the subsurface GSD was truncated at 2 mm.

Figure 4 (a) Relationship between the reach-averaged shear stress (〈τ〉), 1D total shear stress (τ_t), and discharge (Q). (b) Relationship between the parameters of the gamma function (α and β) and discharge. Total shear stress from depth slope product has a similar trend as 〈τ〉, but it is consistently higher. All our results are based on the 2D models, and τ_t is shown here only for comparison purposes.

3 Results

3.1 Spatial distributions of shear stress

The numerical models allowed the characterization of the spatial distribution of τ for each discharge level (Fig. 3). In terms of reach-averaged values, the predicted 〈τ〉 varied between 18.3 and 51.1 N m⁻² for flows between 0.12 and 1.0 Q_bf (Table 1). Furthermore, the mean shear stress, 〈τ〉, scaled with discharge such that an exponential function explained 97 % of its variance (Fig. 4a). The predicted 〈τ〉 were 66 % to 79 % smaller than the mean shear stress values calculated based on the depth–slope product (Table 1).

Table 1 Summary of model shear stress distributions including reach-averaged shear stress (〈τ〉), mean modeled water depth (h), and total shear stress (τ_t=ρghS_b); measured bed load transport rate (q_{b_meas}) (Milhous, 1973), gamma fit parameters (α and β), and goodness of fit: reduced chi-square (${\mathit{\chi }}_v^{\mathrm{2}}$) and root mean square error (RMSE) between the observed distribution of shear stress and the gamma fit predicted distribution.

^* We estimated the total shear stress (τ_t) assuming uniform flow (depth–slope product) and a constant energy slope of 0.014 m m⁻¹.

The shape of the distributions of τ∕〈τ〉 was remarkably similar across all modeled discharges (Fig. 3). In all cases the highest frequencies of local τ were around the mean value, and approximately 92 % of the predicted τ∕〈τ〉 were below 2. We fitted the normalized shear stress distributions to gamma functions with α parameters that varied between 7.49 and 3.60 and β parameters that varied between 0.13 and 0.27 (Table 1). These parameters, α and β, varied linearly with discharge (Fig. 4b). In both cases discharge explained more than 92 % of the variability in α and β.

The equations that relate the gamma fit parameters and the reach-averaged shear stress to the discharge are

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }=-\mathrm{1.155}Q+\mathrm{7.329}\text{(14)}& {\displaystyle }& {\displaystyle }\mathit{\beta }=\mathrm{0.044}Q+\mathrm{0.120}\text{(15)}& {\displaystyle }& {\displaystyle }\langle \mathit{\tau }\rangle =\mathrm{17.288}{e}^{\mathrm{0.319}Q}.\end{array}$$

Combining equations (Eqs. 13 and 15) with Eq. (2), an expression to estimate the distribution of τ for any given Q can be obtained. Additionally, these synthetic distributions can be used to evaluate the accuracy of our bed load transport equation for discharge levels different than those used for its calibration (see Sect. 3.3).

3.2 Characteristics of our sediment transport relation

The proposed sediment transport equation has the same shape as the Parker and Klingeman (1982) relation, but it is scaled such that $W_i^{*}$ is consistently lower for all ϕ_i values (Fig. 5). The consistently lower $W_i^{*}$ indicates that bed load transport occurs are relatively small localized areas of the bed where stress stresses are higher than the average value. While calibrating this formulation, we kept some key features of the original equation in Parker and Klingeman (1982); thus, we reduced the number of degrees of freedom. Specifically, the shape of both equations is the same, and each is valid within the same ϕ_i intervals. We used an exponential function with a second-degree polynomial function as argument for $\mathrm{0.95}<{\mathit{\varphi }}_i<\mathrm{1.65}$ and a power function with an exponent equal to 4.5 for ϕ_i≥1.65. We also maintained ${\mathit{\tau }}_{\mathrm{r}\mathrm{50}}^{*}=\mathrm{0.0876}$ and the exponent of the hiding function (0.982) as fixed values (i.e., were not adjusted while calibrating our equations).

Figure 5 Comparison between different subsurface-based sediment transport equations and the one proposed in this study. The relation of Segura and Pitlick (2015), which is also a modified version of Parker and Klingeman (1982), is shown as reference.

3.3 Sediment transport calculations

All flow events used for calibrating had an error of less than an order of magnitude between the measured and predicted bed load transport rate (Fig. 6, Table 2). In terms of sediment transport estimates, this order of error is generally considered as a relatively strong estimation (Yager et al., 2007, 2012b). Similar to the Parker and Klingeman (1982) equation, our bed load estimates for flows lower than 0.4 m³ s⁻¹ were weaker. This is not surprising given that these low flows were not used for calibration and that there are very low rates of transport at such low discharges (∼10 % of Q_bf). The equation of Parker and Klingeman (1982) was not designed to include distributions of τ. However, to have a point of comparison, we contrasted the measured and predicted bed load rates for the original Parker and Klingeman (1982) equation applied over the complete shear stress distribution predicted by the 2D numerical model instead of our formulation. While our estimated bed load rates for Q≥0.99 m³ s⁻¹ were within 1 order of magnitude of the observed value, those predicted using the Parker et al. (1982) equation were consistently overestimated by over an order of magnitude in all cases (Fig. 6).

Figure 6 Comparison between measured and predicted bed load transport rate for different methods and datasets. Five events (Q≥0.99 m³ s⁻¹) were used when calibrating Eq. (8) (black circles). Triangles represent the estimated bed load using Eq. (8) for two low-flow events (Q<0.64 m³ s⁻¹) that were not used for calibration. The equation of Parker and Klingeman (1982) applied locally to the complete shear stress distributions is shown as reference (squares). Additionally, a synthetic spatial shear distribution based on Eq. (2) and the parameters given in Eqs. (13)–(15) was used with our equation to calculate the bed load rate (grey circles). Measured field data were collected by Milhous (1973).

Table 2 Modeled (q_{b_pred}) and observed (q_{b_meas}) bed load transport rates and modeled (D_{50_pred}), Nash–Sutcliffe efficiency index, and observed (D_{50_meas}) median grain size for the events used in the calibration of Eq. (8).

In terms of prediction of the bed load GSD, the Nash–Sutcliffe efficiency index was in all cases greater than 0.65 (Table 2). For Q≥1.33 m³ s⁻¹ the efficiency index was greater than or equal to 0.85 (Table 2). The difference between predicted and observed bed load median grain sizes ($D_{{\mathrm{50}}_{\mathrm{meas}}}$–$D_{{\mathrm{50}}_{\mathrm{pred}}}$) was lower than 10 mm in all these cases. For the Q=0.99 m³ s⁻¹ event, the error in the median grain size was larger (12.3 mm), with predicted grain size values consistently coarser (Fig. 7).

Figure 7 Measured and predicted bed load grain size distributions for all events used in the calibration of Eq. (8). Flow discharges are shown in the lower-left corner of each panel.

Equation (8) is only applicable when the spatial distribution of τ is known. However, this is not the case in most studies and practical applications. In our case, given the strong correlations between discharge and reach-averaged shear stress and also between discharge and the gamma function parameters, combining Eqs. (13) and (15) with Eq. (2) allowed us to generate synthetic distributions of τ for a given flow of interest. We tested the accuracy of our equation when these synthetic distributions were used as input using a subset of the Milhous (1973) database (grey circles in Fig. 6). The scenarios considered correspond to the same 22 flow events used by Parker and Klingeman (1982) in their analysis and had flow discharges that ranged between 1.02 and 3.4 m³ s⁻¹. Using the synthetic distributions of τ, our equation predicted bed load rates within an order magnitude of error for all 22 events. Considering the logarithm of the ratio between the measured and predicted bed load transport rate, log (q_{b_pred}∕q_{b_meas}), which is a measure of the accuracy of an estimation (0 indicates perfect agreement and ±1 an error of an order of magnitude; Yager et al., 2007), the estimated bed load rates had a median log (q_{b_pred}∕q_{b_meas}) of −0.07, minimum of −0.84, 25th percentile of −0.42, 75th percentile of 0.17, and a maximum of 0.55.

Based on the local shear stress, we identified the areas of the bed where most of the bed load likely occurs for a given flow level (Fig. 8). Similar to the study of Segura and Pitlick (2015), in terms of bed load transport rate per unit width, the size of these areas increases with discharge. At 0.29 Q_bf (Q=0.99 m³ s⁻¹), most bed load transport occurs in a relatively small, localized area of the bed, occupying approximately 5.4 % of the total bed area (100 % is the wetted area under bankfull flow conditions). The percentage increases to 7.6 % at 0.43 Q_bf (Q=1.46 m³ s⁻¹) and 17.5 % for 0.56 Q_bf (Q=1.91 m³ s⁻¹). At bankfull flow conditions (Q=3.40 m³ s⁻¹), the proportion of the channel bed that is predicted to be mobile is 52.5 % and mainly concentrated along the thalweg. It is important to mention that this method only provides an approximation of the region where most bed load occurs. In this approach we are not considering bed evolution and other time-dependent processes that may alter the location of areas where bed load transport occurs.

Figure 8 Predicted local bed load transport rate per unit width (q_b) for four flows levels between 0.29 and 1.00 Q_bf. White areas in each case indicate zero predicted transport. Black lines correspond to the wetted area under bankfull flow conditions. Maps of the predicted local shear stress for the same study area and flow scenarios are available in Katz et al. (2018).

4 Discussion

4.1 Using spatial distribution of shear stress to estimate reach-average bed load rates

Existing transport equations represent the available shear stress for a given flow with a single shear stress value, usually the reach-averaged total shear stress estimated as the depth–slope product (Barry et al., 2004; Fernandez Luque and Van Beek, 1976; Meyer-Peter and Müller, 1948; Parker, 1990; Parker et al., 1982; Parker and Klingeman, 1982; Recking, 2013b; Wilcock and Crowe, 2003; Wilcock and Kenworthy, 2002). In other words, most of the available equations assume that this single shear stress value represents the entire shear stress distribution for any given flow level. While this assumption may be appropriate in certain cases, for example in straight reaches with few roughness elements, it is unlikely to represent the hydraulic conditions in complex reaches with variable planforms and roughness characteristics. In our approach, we explicitly account for the local variability in bed surface elevation, channel curvature, and roughness characteristics by including spatially variable estimates of shear stress over a range of hydraulic conditions within the reach, making it more applicable to a wide range of stream types. The main difference between the transport function proposed in Eq. (8) and those typically used when estimating bed load transport rates (e.g., Parker, 1990; Parker et al., 1982; Recking, 2013b; Wilcock and Crowe, 2003) is that Eq. (8) uses the full distribution of shear stress rather than the reach-averaged shear stress value for a given flow. In practical applications, both approaches require the same input data, specifically, a given discharge, measure of bed roughness, GSD, and bed surface elevation. While it may be enough for equations using the reach-averaged τ to define energy gradients using a longitudinal bed profile, our method requires detailed measurements of bed topography to adequately construct a numerical 2D flow model to estimate spatial shear stress distributions. Although acquiring detailed bed surface topography may be restrictive, this method offers an alternative to modern approaches that rely on detailed field measurements to estimate the τ applied to the mobile sediment fractions of a given bed. Current flow resistance and shear stress partitioning techniques used in mountain river applications require a characterization of the macro-roughness (Nitsche et al., 2011, 2012) that involves careful field measurements of the diameter, protrusion, concentration, and spacing of boulders (e.g., Monsalve et al., 2016; Yager et al., 2012a), length, slope, and height of steps (Nitsche et al., 2011), and every other source of roughness beside skin friction. Therefore, in general terms, comparable field effort is required for both modelling of shear stress and estimating of shear stress partition.

We modified a subsurface-based equation to include the spatial distribution of shear stress (i.e., Parker and Klingeman, 1982). Alternatively, we could have chosen to modify the surface-based equation of Parker (1990), also developed using data from Oak Creek, because from a mechanistic point of view, it is the bed surface that is in contact with the water, whereas the subsurface is not always directly accessed by the flow. However, the fits for the larger number of parameters in the Parker (1990) approach would have been weaker considering the small number of flow events with sufficient information of both the bed load GSD and spatial distribution of shear stress (Table 2) (see also Sect. 4.2). Nonetheless, future improvements to our approach could consider the use of a surface-based equation (e.g., Parker, 1990, or Wilcock and Crowe, 2003).

In our equation we used a reach-averaged GSD. Recent studies have shown that including the local τ^*, based on local shear stress and grain size characteristics, can improve sediment transport predictions in complex mountain rivers (e.g., Monsalve et al., 2016). However, we used a reach-averaged GSD in this study because of the following: (i) measuring local grain size distributions (or sediment patches) in a given river is practically complicated for developing a method broadly applicable – this is especially true when trying to delineate submerged sediment patches; (ii) the GSD over a reach may vary spatially, but the reach-averaged GSD of a given reach is less sensitive to changes in discharge than the shear stress – Segura and Pitlick (2015) compared the variability in the shear stress distribution and the grain size distribution and found that the shear stress distributions varies more than the GSD; and (iii) spatial-scale modeling restrictions. Two-dimensional models are not able to incorporate the effects of fine-scale variability in the surface grain size. Usually the grid cell sizes in these models are on the order of 20–50 cm. Therefore, even if a detailed grain size distribution were available, fully coupling them within a 2D approach is not yet possible.

4.2 Alternative formulations for sediment transport prediction using spatial distribution of shear stress

When calibrating Eq. (8), we used a total of five flow levels covering a wide range of discharges, from the lower limit (approximately 0.29 Q_bf) used by Parker and Klingeman (1982) up to bankfull conditions. While conducting the calibration, we found an alternative formulation defined by a single equation (instead of a two-part equation like Eq. 8), also calibrated for flows above 0.99 m³ s⁻¹. This equation performed well over a wider range of flows, including those between 0.4 and 0.64 m³ s⁻¹.

$$\begin{array}{}\text{(16)}& W_i^{*}=\mathrm{0.38}{\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1.5}}{{\mathit{\varphi }}_i}}\right)}^{\mathrm{1.5}},\mathrm{for}{\mathit{\varphi }}_i\ge \mathrm{1.5}\end{array}$$

The performance of Eqs. (8) and (16) in terms of predicted bed load transport rates and GSD was relatively similar for Q≥0.99 m³ s⁻¹ (Table 3; for simplicity only D₅₀ is shown). However, Eq. (16) also predicted q_b and GSD well for all discharges lower than 0.99 m³ s⁻¹, with errors below an order of magnitude. When Eq. (8) is applied to the 0.4 m³ s⁻¹ flow, it overestimates the measured bed load rate by 27 times (Fig. 6). It is important to remark that in the calibration process of Eqs. (8) and (16), the discharge levels of 0.4 and 0.64 m³ s⁻¹ were not used. We presented Eq. (8) in the “Results” section because it resembles the Parker and Klingeman (1982) equation to which we were comparing (Fig. 6). However, from a practical perspective, either formulation could have been used for Q>0.99 m³ s⁻¹. The ability of Eq. (16) to accurately capture low-flow events is explored in detail in Sect. 4.3.

Table 3 Bed load transport rates (q_b) and median grain size estimates (D₅₀) using Eqs. (8) and (16).

4.3 Comparison between Eqs. (8) and (16) and Parker and Klingeman (1982)

Not surprisingly, the subsurface-based sediment transport equation of Parker and Klingeman (1982) gives accurate estimates of bed load for flow events capable of breaking the pavement in a certain reach, given that the equation was exclusively developed for those conditions. Since we are presenting a new approach for estimating bed load transport rates, we compared the performance of Eqs. (8) and (16) to the Parker and Klingeman (1982) equation. First, we studied the accuracy of these three methods for 27 events with flow discharges larger than 1 m³ s⁻¹ (Fig. 9a). All approaches had practically an equal performance when predicting these sediment transport events and had estimates within an order of magnitude of error (Fig. 9b). The equations of Parker and Klingeman (1982) and Eq. (8) predicted a total of 16 events (59 %) within factor of 2 (between 0.5 and 2 times the measured bed load rate), whereas and Eq. (16) predicted 14 events within this range (52 %). Compared to Parker and Klingeman (1982) Eqs. (8) and (16) underpredicted bed load for most of the events but had a slight improvement in terms of the RMSE of the predicted bed load transport rate (Fig. 9b).

One limitation of the Parker and Klingeman (1982) equation is that it is valid only for ϕ>0.95. In practical terms, a value of ϕ=0.95 in Oak Creek is close to the already mentioned discharge of 1 m³ s⁻¹. This relatively high value introduces a practical limitation in the applicability of this method because low discharges are more frequent than high-flow events. According to 266 observed sediment transport events in Oak Creek, including the data of Milhous (1973) and measurements collected in 1978–1990, the majority of the monitored events (∼86 %) were at discharges below 1 m³ s⁻¹. In all these cases (230 events) a bed load transport rate was measured. Using this dataset we tested the performance of Eqs. (8) and (16) for predicting low- and high-flow events that vary between 0.01 and 3.4 m³ s⁻¹. Given that our equations use the distribution of shear stress, they, theoretically, should predict sediment transport even at relatively low flows, and, by doing so, they would overcome the limitation of the Parker and Klingeman (1982) formulation.

Figure 9 (a) Comparison between measured and predicted bed load transport rate using the Parker and Klingeman (1982) equation and Eqs. (8) and (16). In this case, Parker and Klingeman (1982) was applied as proposed in the original publication (i.e., using reach-averaged flow properties). Equations (8) and (16) use spatial distributions of τ obtained with a gamma function and α and β parameters varying with Q. (b) The log of the ratio of predicted to measured sediment bed load rate for the three approaches. A value of zero indicates that the measured volume was predicted exactly. The top and bottom of each box are the 25th and 75th percentiles, and the middle line inside the box is the median value. Lines extending out of the box correspond to the maximum and minimum predicted bed load ratios. The rate at the top of each box corresponds to the RMSE of the predicted bed load rate.

Figure 10 (a) Measured and predicted bed load transport rates as a function of discharge. Equations (8) and (16) can be represented as a continuous line because the spatial distributions of τ and the α and β parameters vary with Q. Predictions of Parker and Klingeman (1982) were calculated using the reach-averaged shear stress based on Milhous (1973) measurements. Therefore, shear stress does not monotonically increase with larger discharges. (b) Measured versus predicted bed load transport rate using Eq. (16). (c) The log of the ratio of predicted to measured sediment bed load rate for Eqs. (8) and (16). Grey arrows extending out of the box correspond to the number of events under- or overpredicted bed load by more than an order magnitude error.

Equations (8) and (16) predicted relatively similar bed load rates for discharges above 0.8 m³ s⁻¹ (Fig. 10a). For Q<0.8 m³ s⁻¹, the equations behave differently. Equation (8) had consistently larger q_b compared to Eq. (16). The difference between Eqs. (8) and (16) increased as flow discharge decreased, and the maximum difference was about 15 times for Q=0.01 m³ s⁻¹ (Fig. 10a). We found that 91 % of the observed sediment transport events in Oak Creek were predicted within an order magnitude (Fig. 10b) with Eq. (16). In general, Eq. (8) underpredicted q_b, while Eq. (16) overpredicted q_b. Specifically, Eq. (8) overpredicted bed load for 72 % of the discharge events, while Eq. (16) underpredicted bed load for 67 % of the discharge events (Fig. 10c). These discrepancies were also reflected in the distribution of the ratio of predicted to measured bed load rate values (Fig. 10c). Considering the logarithm of this ratio, for Eq. (8) the 25th, 50th, and 75th percentiles were 0.52, 0.92, and 0.04, while for Eq. (16) the 25, 50th, and 75th percentiles were −0.24, 0.14, and −0.58 (Fig. 10c). Contrary to Parker and Klingeman (1982), Eqs. (8) and (16) were able to predict q_b at low flows because the lower limit of ϕ=0.95 in Eq. (8) or ϕ=1.5 in Eq. (16) did not correspond to a given discharge (1 m³ s⁻¹ in the case of Parker and Klingeman, 1982). Instead, when using our Eqs. (8) and (16), ϕ varies locally with Q such that it captures high values of τ that occur even at very low flows in small portions of the bed.

4.4 Practical and management implications

The ability of the proposed transport Eqs. (8) and (16) to accurately predict bed load transport rate at a wide range of flows allows our approach to be applied across many different practical scenarios. For small streams like Oak Creek (less than ∼10 m in bankfull width) with relatively simple channel geometry and low relative roughness, Eqs. (13)–(15) can be combined with Eq. (16) to estimate q_b across a range of flow levels and without a 2D hydraulic model. Equations (13)–(15) can first be used to estimate the τ distribution for a given discharge level, and then Eq. (16) can be used with that distribution to estimate q_b. Because streams of this type are fairly ubiquitous in modern urban and suburban society, this method can be applied to a range of management situations such as addressing elevated sediment loads caused by urbanization or glacial retreat. For larger streams and rivers, our approach can be utilized in conjunction with the development of a 2D hydraulic model to accurately estimate sediment transport using either Eq. (8) or Eq. (16). In all situations, our approach is an improvement on previous methods in predicting bed load transport for lower flow levels. This is especially important because it allows for practitioners to better predict the responses of management actions on sediment transport dynamics for these more frequent flow levels. It should be noted that, although our method could be capable of predicting fluxes with better accuracy than previous approaches, all our results are based on measurements in a single river reach of Oak Creek. Therefore, we would recommend using this method with caution until it has been further tested in other systems.

5 Conclusions

Compared to traditional subsurface sediment transport equations that use reach-averaged properties, the proposed equations were able to accurately predict the observed bed load rates at a wider range of flow levels. The shape of the spatial distribution of shear stress was relatively similar for different discharges and allowed us to characterize it in terms of a gamma function. Therefore, we were able to extend our results to scenarios were no field measurements were made. Nonetheless, increasing the accuracy in bed load estimates requires additional efforts compared to the most approaches (i.e., reach-averaged equations). Specifically, the method proposed relies on detailed numerical flow modelling and field measurements, which can restrict the applicability in typical practical studies. However, this may not be a limitation for its use. Considering that realistic estimates of flow resistance in gravel-bed rivers require a characterization of the all sources of roughness, including macro-roughness elements, both approaches need similar field effort, which is, from a practical point of view, the most time-consuming process. In our method, accurate estimates of bed load transport rates at low-flow discharge were possible because we explicitly considered high values of τ, even though they occur in small portions of the bed. Future lines of work should include the extension of surface-based bed load equations and exploring how the shape of the spatial distribution of shear stress varies in other rivers with different geomorphological conditions (e.g., step-pool morphologies, steeper slopes, bed surface patchiness, etc.).