Articles | Volume 8, issue 3
Earth Surf. Dynam., 8, 825–839, 2020
Earth Surf. Dynam., 8, 825–839, 2020

Research article 29 Sep 2020

Research article | 29 Sep 2020

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

A bed load transport equation based on the spatial distribution of shear stress – Oak Creek revisited
Angel Monsalve1, Catalina Segura2, Nicole Hucke1, and Scott Katz2,3 Angel Monsalve et al.
  • 1Departamento de Ingeniería de Obras Civiles, Universidad de la Frontera, Francisco Salazar 1145, Temuco, Chile
  • 2Forest Engineering, Resources & Management, Oregon State University, 201 Peavy Hall, Corvallis, OR 97331, USA
  • 3Natural Systems Design, Bellingham, WA 98225, USA

Correspondence: Angel Monsalve (


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 m3 s−1 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 (Qb) based on a single grain size (the median – D50), 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 Qb 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 Qb 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 1Location of the study reach in Oak Creek, Oregon (442319.092′′ N, 1231951.312′′ W). Contours every 0.1 m are indicated.

Figure 2Surface 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 km2 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 (Sb) of 0.014 m m−1 and bankfull dimensions of 6 m in width and 0.46 m in depth. Recent field observations indicate that bankfull discharge (Qbf) is 3.4 m3 s−1 (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 D50 is 45 mm, while the subsurface D50_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 m3 s−1, equivalent to 0.12 Qbf to Qbf) 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 (Cd). 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, Cd 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 R2>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 Cd, the vertically averaged streamwise (u) and cross-stream (v) velocities, and water density (ρ), assumed as 1000 kg m−3.

(1) τ x y = ρ C d u x y 2 + v x y 2 ,

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):

(2) f ( τ ) = α α ( τ / τ ) ( α - 1 ) e - α ( τ / τ ) τ Γ ( α ) ,

where Γ is the standard gamma function, α is the shape parameter, and β=τ/α 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 (χv2), defined as chi square (χ2) divided by the number of degrees of freedom, according to

(3) χ 2 = f k - f x k 2 σ k 2 ,

where and fk and f(xk) are observed and predicted mean-normalized shear stress frequencies in a given bin interval, k. The uncertainty associated with the observed frequencies, σk2, 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 χv21 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 m3 s−1 or higher to calibrate the new equations (for calibration purposes, our lower discharge was 0.99 m3 s−1). 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 (τ*):

(4) τ * = τ ρ s - ρ g D ,

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),

(5) G = W * W r * = 5.6 × 10 - 3 1 - 0.853 ϕ 4.5 ,

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

(6) ϕ = τ * / τ r * .

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

(7) W * = ( s - 1 ) g q ( τ / ρ ) 1.5 ,

where s is the specific gravity of sediment (s=ρs/ρ).

Figure 3Frequency 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 (Di; subscript i denotes the size range) and ϕi>0.95. In the most general form, the equation for the dimensionless transport rate is Wi*=0.0025Gi, where the constant is the reference transport rate of Parker and Klingeman (1982) (Wr*=0.0025), and Gi is the new (modified) transport relation. The proposed relation is a two-part equation applicable to sediment mixtures:

(8) W i * = 0.0025 × 10 - 3 exp 26.6 ϕ i - 1 - 19.53 ϕ i - 1 2 , for 0.95 < ϕ i < 1.65 W i * = 0.57 1 - 0.853 ϕ i 4.5 , for ϕ i 1.65 .

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

(9) τ r i * τ r 50 * = D i D 50 - 0.982 ,

where τr50*=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 Di and for the entire distribution of shear stress values (i.e., every τi*), was rewritten as

(10) ϕ i = τ i * / τ r i * .

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−2). For all discharges, τj was defined such that it ranges from zero to the maximum predicted shear stress value. For a given Di and τj, the volumetric transport rate per unit width (qij) is

(11) q i j = τ j ρ 1.5 F i W i j ( s - 1 ) g f τ j ,

where Fi is the volume fraction of the ith grain-size class in the subsurface GSD, Wij is calculated using Eq. (8) for each τj, and fτ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

(12) q b = b i j q i j ,

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−2 for flows between 0.12 and 1.0 Qbf (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 1Summary of model shear stress distributions including reach-averaged shear stress (τ〉), mean modeled water depth (h), and total shear stress (τt=ρghSb); measured bed load transport rate (qb_meas) (Milhous, 1973), gamma fit parameters (α and β), and goodness of fit: reduced chi-square (χv2) 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−1.

Download Print Version | Download XLSX

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


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 Wi* is consistently lower for all ϕi values (Fig. 5). The consistently lower Wi* 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 0.95<ϕi<1.65 and a power function with an exponent equal to 4.5 for ϕi≥1.65. We also maintained τr50*=0.0876 and the exponent of the hiding function (0.982) as fixed values (i.e., were not adjusted while calibrating our equations).

Figure 5Comparison 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 m3 s−1 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 Qbf). 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 m3 s−1 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 6Comparison between measured and predicted bed load transport rate for different methods and datasets. Five events (Q≥0.99 m3 s−1) 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 m3 s−1) 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 2Modeled (qb_pred) and observed (qb_meas) bed load transport rates and modeled (D50_pred), Nash–Sutcliffe efficiency index, and observed (D50_meas) median grain size for the events used in the calibration of Eq. (8).

Download Print Version | Download XLSX

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 m3 s−1 the efficiency index was greater than or equal to 0.85 (Table 2). The difference between predicted and observed bed load median grain sizes (D50measD50pred) was lower than 10 mm in all these cases. For the Q=0.99 m3 s−1 event, the error in the median grain size was larger (12.3 mm), with predicted grain size values consistently coarser (Fig. 7).

Figure 7Measured 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 m3 s−1. 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 (qb_predqb_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 (qb_predqb_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 Qbf (Q=0.99 m3 s−1), 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 Qbf (Q=1.46 m3 s−1) and 17.5 % for 0.56 Qbf (Q=1.91 m3 s−1). At bankfull flow conditions (Q=3.40 m3 s−1), 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 8Predicted local bed load transport rate per unit width (qb) for four flows levels between 0.29 and 1.00 Qbf. 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 Qbf) 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 m3 s−1. This equation performed well over a wider range of flows, including those between 0.4 and 0.64 m3 s−1.

(16) W i * = 0.38 1 - 1.5 ϕ i 1.5 , for ϕ i 1.5

The performance of Eqs. (8) and (16) in terms of predicted bed load transport rates and GSD was relatively similar for Q≥0.99 m3 s−1 (Table 3; for simplicity only D50 is shown). However, Eq. (16) also predicted qb and GSD well for all discharges lower than 0.99 m3 s−1, with errors below an order of magnitude. When Eq. (8) is applied to the 0.4 m3 s−1 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 m3 s−1 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 m3 s−1. The ability of Eq. (16) to accurately capture low-flow events is explored in detail in Sect. 4.3.

Table 3Bed load transport rates (qb) and median grain size estimates (D50) using Eqs. (8) and (16).

Download Print Version | Download XLSX

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 m3 s−1 (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 m3 s−1. 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 m3 s−1. 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 m3 s−1. 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 m3 s−1 (Fig. 10a). For Q<0.8 m3 s−1, the equations behave differently. Equation (8) had consistently larger qb 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 m3 s−1 (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 qb, while Eq. (16) overpredicted qb. 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 qb 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 m3 s−1 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 qb 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 qb. 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.).

Data availability

FaSTMECH predicted distributions of depth, velocity, and shear stress for seven flow levels are available at the ScholarsArchive@OSU ( (Segura and Katz, 2020). The bed load data are published in Milhous (1973).

Author contributions

CS proposed the original idea; AM and CS led the organization of the study; CS acquired the funding. Both SK and CS contributed to the field data collection. SK conducted the flow modeling under the supervision of CS. NH contribute to the analysis under the supervision of AM. AM, CS, NK, and SK contributed to the interpretation of results and to the writing.

Competing interests

The authors declare that they have no conflict of interest.


We thank Pete Klingeman who led the Oak Creek Sediment Transport facility, which collected the database used as the basis for this study. Without Pete, none of this would have been possible. This research was funded by the National Science Foundation (award no. 1619700) and the USDA National Institute of Food and Agriculture (McIntire Stennis project OREZ-FERM-876). We also thank Allison Pfeiffer and James Pizzuto for their constructive comments that helped improve an earlier version of this paper. The authors would like to thank the McDonald Dunn Forest Director, Stephen Fitzgerald, for logistical support while setting our study site. Richard McDonald of the USGS provided very useful suggestions for the hydrodynamic model FaSTMECH.

Financial support

This research has been supported by the National Science Foundation (grant no. 1619700) and the USDA National Institute of Food and Agriculture (grant no. OREZ‐FERM‐876).

Review statement

This paper was edited by Lina Polvi Sjöberg and reviewed by James Pizzuto and Allison Pfeiffer.


Barry, J. J., Buffington, J. M., and King, J. G.: A general power equation for predicting bed load transport rates in gravel bed rivers, Water Resour. Res., 40, 1–22,, 2004. 

Barton, G. J., McDonald, R. R., Nelson, J. M., and Dinehart, R. R.: Simulation of flow and sediment mobility using a multidimensional flow model for the white sturgeon critical-habitat reach, Kootenai River near Bonners Ferry, Idaho, US Geological Survey Scientific Investigations Report 2005-5230, US Geological Survey, Reston, Virginia, p. 54, 2005. 

Beschta, R. L.: Increased bag size improves Helley-Smith bed load sampler for use in streams with high sand and organic matter transport, in: Erosion and sediment transport measurement, Proceedings of the Florence symposium IAHS, June 1981, Florence, 17–25, 1981. 

Bevington, P. R. and Robinson, D. K.: Data reduction and error analysis for the physical sciences, 3rd Edition, McGraw-Hill, New York, NY, USA, 2003. 

Bradley, D. N. and Tucker, G. E.: Measuring gravel transport and dispersion in a mountain river using passive radio tracers, Earth Surf. Proc. Land., 37, 1034–1045,, 2012. 

Bunte, K., Abt, S. R., Potyondy, J. P., and Swingle, K. W.: A comparison of coarse bedload transport measured with bedload traps and Helley-Smith samplers, Geodin. Acta, 21, 53–66,, 2008. 

Clayton, J. A. and Pitlick, J.: Spatial and temporal variations in bed load transport intensity in a gravel bed river bend, Water Resour. Res., 43, 1–13,, 2007. 

Conner, J. T. and Tonina, D.: Effect of cross-section interpolated bathymetry on 2D hydrodynamic model results in a large river, Earth Surf. Proc. Land., 39, 463–475,, 2014. 

Einstein, H. A.: The bed-load function for sediment transportation in open channel flows, US Soil Conservation Service Technical Report 1026, US Soil Conservation Service, Washington, D.C.,, 1950. 

Emmett, W. W.: A field calibration of the sediment-trapping characteristics of the Helley-Smith bedload sampler, USGS Professional paper 1139, USGS, Washington, D.C.,, 1980. 

Ferguson, R. I.: The missing dimension: Effects of lateral variation on 1-D calculations of fluvial bedload transport, Geomorphology, 56, 1–14,, 2003. 

Fernandez Luque, R. and Van Beek, R.: Erosion And transport Of bed-load sediment, J. Hydraul. Res., 14, 127–144,, 1976. 

Freedman, D. and Diaconis, P.: On the histogram as a density estimator: L2 theory, Z. Wahrscheinlichkeitstheor. Verw. Gebiet., 57, 453–476,, 1981. 

Gomez, B. and Church, M.: An assessment of bed load sediment transport formulae for gravel bed rivers, Water Resour. Res., 25, 1161–1186,, 1989. 

Hassan, M. A. and Church, M.: Sensitivity of bed load transport in Harris Creek: Seasonal and spatial variation over a cobble-gravel bar, Water Resour. Res., 37, 813–825,, 2001. 

Katz, S. B., Segura, C., and Warren, D. R.: The influence of channel bed disturbance on benthic Chlorophyll a: A high resolution perspective, Geomorphology, 305, 141–153,, 2018. 

Kinzel, P. J., Nelson, J. M., and Heckman, A. K.: Response of sandhill crane (Grus canadensis) riverine roosting habitat to changes in stage and sandbar morphology, River Res. Appl., 25, 135–152,, 2009. 

Lisle, T. E.: Particle size variation between bed load and bed material in natural gravel beds channel, Water Resour. Res., 31, 1107–1118,, 1995. 

Lisle, T. E., Nelson, J. M., Pitlick, J., Madej, M. A., and Barkett, B. L.: Variability of bed mobility in natural, gravel-bed channels and adjustments to sediment load at local and reach scales, Water Resour. Res., 36, 3743–3755,, 2000. 

Maturana, O., Tonina, D., McKean, J. A., Buffington, J. M., Luce, C. H., and Caamaño, D.: Modeling the effects of pulsed versus chronic sand inputs on salmonid spawning habitat in a low-gradient gravel-bed river, Earth Surf. Proc. Land., 39, 877–889,, 2014. 

May, C. L. and Pryor, B. S.: Initial motion and bedload transport distance determined by particle tracking in a large regulated river, River Res. Appl., 30, 508–520,, 2014. 

May, C. L., Pryor, B., Lisle, T. E., and Lang, M.: Coupling hydrodynamic modeling and empirical measures of bed mobility to predict the risk of scour and fill of salmon redds in a large regulated river, Water Resour. Res., 45, W05402,, 2009. 

McDonald, R., Nelson, J., Paragamian, V., and Barton, G.: Modeling the effect of flow and sediment transport on white sturgeon spawning habitat in the Kootenai River, Idaho, J. Hydraul. Eng., 136, 1077–1092,, 2010. 

Mcdonald, R. R., Nelson, J. M., and Bennett, J. P.: Multi-Dimensional Surface-Water Modeling System User's Guide, US Geological Survey Techniques and Methods 6-B2, US Geological Survey, Reston, VA, 2005. 

Meyer-Peter, E. and Müller, R.: Formulas for Bed-Load Transport, in: Proceedings of the 2nd Meeting of the International Association of Hydraulic Research, 7–9 June 1948, Stockholm, Sweden, 39–64, 1948. 

Milhous, R. T.: Sediment transport in a gravel-bottomed stream, PhD Thesis, Department of Civil Engineering, Oregon State University, Oregon, 1973. 

Miller, A. J. and Cluer, B. L.: Modeling Considerations for Simulation of Flow in Bedrock Channels, in: Bedrock Channels, Rivers Over Rock: Fluvial Processes in Bedrock Channels, edited by: Tinkler, K. J. and Wohl, E. E., American Geophysical Union, Washington, D.C., USA, 1998. 

Monsalve, A. D., Yager, E. M., Turowski, J. M., and Rickenmann, D.: A probabilistic formulation of bed load transport to include spatial variability of flow and surface grain size distributions, Water Resour. Res., 52, 3579–3598,, 2016. 

Mueller, E. R. and Pitlick, J.: Sediment supply and channel morphology in mountain river systems: 2. Single thread to braided transitions, J. Geophys. Res.-Earth, 119, 1516–1541,, 2014. 

Mueller, E. R., Pitlick, J., and Nelson, J. M.: Variation in the reference Shields stress for bed load transport in gravel-bed streams and rivers, Water Resour. Res., 41, 1–10,, 2005. 

Nash, J. E. and Sutcliffe, J. V.: River flow forecasting through conceptual models part I – A discussion of principles, J. Hydrol., 10, 282–290,, 1970. 

Nelson, J. M. and McDonald, R. R.: Mechanics and modeling of flow and bed evolution in lateral separation eddies, US Geological Survey, Glen Canyon Environmental Studies Report, Flagstaff, Arizona, USA, 1995. 

Nelson, J. M. and Smith, J. D.: Flow in meandering channels with natural topography, in: River Meandering, Water Resources Monographs, vol. 12, edited by: Ikeda, S. and Parker, G., American Geophysical Union, Washington, D.C., USA, 69–102, 1989. 

Nelson, J. M., Bennett, J. P., and Wiele, S. M.: Flow and Sediment-Transport Modeling, in: Tools in Fluvial Geomorphology, edited by: Kondolf, G. M. and Piegay, H., John Wiley & Sons, Ltd, Chichester, UK, 539–576, 2003. 

Nelson, P. A., Dietrich, W. E., and Venditti, J. G.: Bed topography and the development of forced bed surface patches, J. Geophys. Res.-Earth, 115, F04024,, 2010. 

Nicholas, A. P.: Modelling bedload yield braided gravel bed rivers, Geomorphology, 36, 89–106,, 2000. 

Nitsche, M., Rickenmann, D., Turowski, J. M., Badoux, A., and Kirchner, J. W.: Evaluation of bedload transport predictions using flow resistance equations to account for macro-roughness in steep mountain streams, Water Resour. Res., 47, W08513,, 2011. 

Nitsche, M., Rickenmann, D., Kirchner, J. W., Turowski, J. M., and Badoux, A.: Macroroughness and variations in reach-averaged flow resistance in steep mountain streams, Water Resour. Res., 48, W12518,, 2012. 

O'Connor, J. E., Mangano, J. F., Anderson, S. W., Wallick, J. R., Jones, K. L., and Keith, M. K.: Geologic and physiographic controls on bed-material yield, transport, and channel morphology for alluvial and bedrock rivers, western Oregon, Bull. Geol. Soc. Am., 126, 377–397,, 2014. 

O'leary, S. J. and Beschta, R. L.: Bedload transport in an Oregon coast range stream, Water Resour. Bull., 17, 886–894,, 1981. 

Olinde, L. and Johnson, J. P. L.: Using RFID and accelerometer-embedded tracers to measure probabilities of bed load transport, step lengths, and rest times in a mountain stream, Water Resour. Res., 51, 7572–7589,, 2015. 

Paola, C.: Incoherent structure: Turbulence as a metaphor for stream braiding, in: Coherent Flow Structures in Open Channels, edited by: Ashworth, P., Bennett, S. J., Best, J. L., and McLelland, S., John Wiley, Chichester, UK, 706–723, 1996. 

Parker, G.: Self-formed straight rivers with equilibrium banks and mobile bed. Part 2. The gravel river, J. Fluid Mech., 89, 127–146,, 1978. 

Parker, G.: Surface-based bedload transport relation for gravel rivers, J. Hydraul. Res., 28, 417–436,, 1990. 

Parker, G. and Klingeman, P. C.: On why gravel bed streams are paved, Water Resour. Res., 18, 1409,, 1982. 

Parker, G., Klingeman, P. C., and McLean, D. G.: Bedload and size distribution in paved gravel-bed streams, J. Hydraul. Div., 108, 544–571, 1982. 

Paustian, S. J. and Beschta, R. L.: The suspended sediment regime of an Oregon coast range stream, J. Am. Water Resour. Assoc., 15, 144–154,, 1979. 

Pitlick, J.: Variability of bed load measurement, Water Resour. Res., 24, 173–177,, 1988. 

Pitlick, J., Mueller, E. R., and Segura, C.: Differences in Sediment Supply to Braided and Single-Thread River Channels: What Do the Data Tell Us?, in: Gravel-Bed Rivers: Processes, Tools, Environments, edited by: Church, M., Biron, P. M., and Roy, A. G., John Wiley & Sons, Ltd., Chichester, UK, 502–511, 2012. 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: Numerical Recipes: The Art of Scientific Computing, 3rd Edn., Cambridge University Press, Cambridge, 2007. 

Recking, A.: A comparison between flume and field bed load transport data and consequences for surface-based bed load transport prediction, Water Resour. Res., 46, 1–16,, 2010. 

Recking, A.: An analysis of nonlinearity effects on bed load transport prediction, J. Geophys. Res.-Earth, 118, 1264–1281,, 2013a. 

Recking, A.: Simple method for calculating reach-averaged bed-load transport, J. Hydraul. Eng., 139, 70–75,, 2013b. 

Rickenmann, D. and McArdell, B. W.: Continuous measurement of sediment transport in the Erlenbach stream using piezoelectric bedload impact sensors, Earth Surf. Proc. Land., 32, 1362–1378,, 2007. 

Schmidt, K.-H. and Ergenzinger, P.: Bedload entrainment, travel lengths, step lengths, rest periods – studied with passive (iron, magnetic) and active (radio) tracer techniques, Earth Surf. Proc. Land., 17, 147–165,, 1992. 

Segura, C. and Katz, S.: FaSTMECH Modelling results for Seven Flow Levels at Oak Creek (Version 1), Oregon State University, Dataset,, 2020. 

Segura, C. and Pitlick, J.: Coupling fluvial-hydraulic models to predict gravel transport in spatially variable flows, J. Geophys. Res.-Earth, 120, 834–855,, 2015. 

Turowski, J. and Rickenmann, D.: Measuring the statistics of bed-load transport using indirect sensors, J. Hydraul. Eng., 137, 116–121,, 2011. 

Vericat, D., Church, M., and Batalla, R. J.: Bed load bias?: Comparison of measurements obtained using two (76 and 152 mm) Helley-Smith samplers in a gravel bed river, Water Resour. Res., 42, W01402,, 2006. 

Wilcock, P. R. and Crowe, J. C.: Surface-based transport model for mixed-size sediment, J. Hydraul. Eng., 129, 120–128,, 2003. 

Wilcock, P. R. and Kenworthy, S. T.: A two-fraction model for the transport of sand/gravel mixtures, Water Resour. Res., 38, 1–12,, 2002.  

Wyss, C. R., Rickenmann, D., Fritschi, B., Turowski, J. M., and Weitbrecht, V.: Laboratory flume experiments with the Swiss plate geophone bed load monitoring system: 1. Impulse counts and particle size identification, Water Resour. Res., 52, 7744–7759,, 2016a. 

Wyss, C. R., Rickenmann, D., Fritschi, B., Turowski, J. M., Weitbrecht, V., Travaglini, E., Bardou, E., and Boes, R. M.: Laboratory flume experiments with the Swiss plate geophone bed load monitoring system: 2. Application to field sites with direct bed load samples, Water Resour. Res., 52, 7760–7778,, 2016b. 

Wyss, C. R., Rickenmann, D., Fritschi, B., Turowski, J. M., Weitbrecht, V., and Boes, R. M.: Measuring bed load transport rates by grain-size fraction using the swiss plate geophone signal at the erlenbach, J. Hydraul. Eng., 142, 1–11,, 2016c. 

Yager, E. M., Kirchner, J. W., and Dietrich, W. E.: Calculating bed load transport in steep boulder bed channels, Water Resour. Res., 43, W07418,, 2007. 

Yager, E. M., Dietrich, W. E., Kirchner, J. W., and McArdell, B. W.: Patch dynamics and stability in steep, rough streams, J. Geophys. Res.-Earth, 117, 1–16,, 2012a. 

Yager, E. M., Dietrich, W. E., Kirchner, J. W., and McArdell, B. W.: Prediction of sediment transport in step-pool channels, Water Resour. Res., 48, W01541,, 2012b. 

Yager, E. M., Venditti, J. G., Smith, H. J., and Schmeeckle, M. W.: The trouble with shear stress, Geomorphology, 323, 41–50,, 2018. 

Short summary
Part of the inaccuracies when estimating bed load transport in gravel-bed rivers is because we are not considering the wide distributions of shear stress in these systems. We modified a subsurface-based bed load transport equation to include these distributions. By doing so, our approach accurately predicts bed load transport rates when the pavement layer is still present, while the original one predicts zero transport. For high flows, our method had similar performance to the original equation.