A comparison of 1D and 2D bedload transport functions under high excess shear stress conditions in laterally-constrained gravel-bed rivers: a laboratory study

. Channel processes under high magnitude flow events are of central interest to river science and management as they may produce large volumes of sediment transport and geomorphic work. However, bedload transport processes under these conditions are poorly understood due to data collection limitations and the prevalence of physical models that restrict feedbacks surrounding morphologic adjustment. The extension of mechanistic bedload transport equations to gravel-bed rivers has emphasised the importance of variance in both entraining (shear stress) and resisting (grain size) forces, especially at low 5 excess shear stresses. Using a fixed-bank laboratory model, we tested the hypothesis that bedload transport in rivers collapses to a more simple function (i.e., with mean shear stress and median grain size) under high excess shear stress conditions. Bedload transport was well explained by the mean shear stress (1D approach) calculated using the depth-slope product. Numerically modelling shear stress to account for the variance in shear stress (2D) did not substantially improve the correlation. Critical dimensionless shear stress values were back-calculated and were higher for the 2D approach compared to the 1D. This result 10 suggests that 2D critical values account for the relatively greater influence of high shear stresses, whereas the 1D approach assumes that the mean shear stress is sufficient to mobilise the median grain size. While the 2D approach may have a stronger conceptual basis, the 1D approach performs unreasonably well under high excess shear stress conditions. Further work is required to substantiate these findings in laterally adjustable channels.


15
The adjustment of rivers to the imposed valley gradient, sediment supply, and discharge is of central interest to geomorphology and has implications for understanding and managing natural hazards and ecological habitats. In alluvial channels, the adjustment is facilitated by the movement of sediment, arising via the interaction between the flow and deformable boundary (Bridge and Jarvis, 1982;Dietrich and Smith, 1983;Church, 2010;Church and Ferguson, 2015). Despite there being no strict correlation between the magnitudes of perturbation and geomorphic effect (Lisenby et al., 2018), larger-than-average flows 20 (i.e., floods) are typically associated with channel adjustment and relatively large volumes of geomorphic work (Wolman and Miller, 1960). Extreme events may exert disproportionate control over the channel planform (Eaton and Lapointe, 2001). The study of sediment transport processes under these relatively high discharge events is central to understanding river behaviour.
Researchers have dedicated considerable effort to deriving mechanistic bedload transport functions -typically empiricallycalibrated -that relate the rate of movement to a force-balance between the flow and individual particles. Other approaches 25 exist, for example, non-threshold approaches that do not utilise a critical shear stress (Recking, 2013a). One of the most simple and widely used threshold relations is the Meyer-Peter and Müller (1948) equation that estimates bedload transport as a function of mean excess bed shear stress (τ − τ c , where τ c is critical shear stress) for a given grain diameter, typically the median (i.e., τ c50 for the D 50 grain). The extension of 1D bedload transport functions to gravel-bed rivers, typically characterised by a wide range of grain sizes, necessitated several modifications that accounted for the differential mobility of grain sizes, hiding 30 and exposure Parker, 1990;Recking, 2013b;Wilcock and Crowe, 2003). Further research emphasised that at conditions whereτ ≈ τ c50 , bedload transport is affected by the spatial variance in shear stress (Paola and Seal, 1995;Paola, 1996;Nicholas, 2000;Ferguson, 2003;Bertoldi et al., 2009;Francalanci et al., 2012;Recking et al., 2016). More recently, Monsalve et al. (2020) proposed a 2D bedload transport function that integrates across the distribution of shear stresses and can predict transport at lower flow conditions whereτ < τ c50 . In concert, these advances demonstrate a 35 consistent trend: with decreasing excess shear stress more information regarding grain size and shear stress (i.e., resisting and driving forces) is required to predict bedload transport.
Considerably less is known about rivers under high relative shear stress conditionsτ ≫ τ c50 , where most channel change occurs. This is primarily due to practical limitations. Dangers associated with floods and erosion mean that researchers may collect data before and after an event, but not during. Laboratory experiments (flumes) typically do not incorporate key degrees-40 of-freedom for morphologic adjustment that are available to alluvial channels, and thus do not model the full range of feedbacks between bedload transport and the deformable boundary. Subsequently, the notion that bedload transport in rivers collapses to a more simple function (i.e., with mean shear stress and median grain size) under high excess shear stress conditions is yet to be conclusively demonstrated. If verified, it would serve as a highly convenient assumption in understanding landscape evolution and river management. Smaller-scale laboratory experiments provide an opportunity to test this hypothesis as they 45 model larger-scale bed and ideally bank adjustments.
We test the relative performance of 1D and 2D bedload transport functions under high relative shear stress conditions in a Froude-scaled physical model. The experiments have a widely-graded sediment mixture and develop alternate bars under pseudo-recirculating conditions at a range of widths and discharges. We record total bedload volumes, bathymetry, and perform 2D hydraulic modelling to apply several transport functions akin to Meyer-Peter and Müller (1948) (i.e., based on median 50 grain size) that capture different levels of information regarding shear stress. The results highlight the effectiveness of simple threshold-based bedload transport functions under high relative shear stress conditions in laterally constrained channels, as well as differences between 1D and 2D conceptualisations of excess shear stress and bedload transport.

Methodology
Experiments were performed in the Adjustable-Boundary Experimental System (A-BES) at the University of British Columbia  The experiments utilised interlocking landscaping bricks to constrict the channel to various widths W between approximately 0.30-0.60 m. In addition to the various channel widths, four different unit discharges (q = Q/W ) were used across the experiments (i.e., discharge was scaled by width) that increased by a factor of 1.5 (Table 1). Two constant-discharge runs  3.00, 4.50 8, 4, 4, 4 20, 16, 16, 16 used the middle two discharges, and one multi-discharge run consisted of the four discharge phases in increasing order. A full list of experiments is provided in Table 2.
At the beginning of each experiment the bulk mixture was mixed by hand to minimise lateral and downstream sorting, and then the in-channel area was screeded to the height of weirs at the upstream and downstream end using a tool that rolled along the brick surface. The flow was run at a low rate with little-to-no movement of sediment until the bed was fully saturated, and 70 was then rapidly increased to the target flow.
Three different types of data were collected throughout each experiment; surface photos, stream gauge measurements, and sediment output. A rolling camera rig positioned atop the A-BES consisted of five Canon EOS Rebel T6i DSLRs with EF-S 18-55 mm lenses (set at 30 mm) positioned at varying oblique angles in the cross-stream direction to maximise coverage of the bed, and five LED lights. Photos were taken in RAW format at 0.2 m downstream intervals, providing a stereographic 75 overlap of over two-thirds. Ten water stage gauges comprised of a measuring tape on flat boards were located along the inner edge of the bricks every 1 m. To minimise edge effects, gauges were not placed within 0.60 m of either the inlet or the outlet.
The gauges were read at an almost horizontal angle which, in conjunction with the dyed blue water, minimised systematic bias towards higher readings due to surface tension effects.
The data collection procedure was designed to maximise measurement accuracy as much as reasonably possible. Given 80 that stream gauge data would later be paired with topographic data, the timing of gauge readings needed to closely coincide with surface photography. Every time photos were taken the bed was drained, as the surface water would distort the photos.
These constraints necessitated a procedure in which manual stream gauge readings (to the nearest 1 mm) were taken 30-40 seconds before the bed was rapidly drained, around the minimum time it would take to obtain the readings. The bed was then photographed and gradually re-saturated before resuming the experiment, approximately 10 minutes.

85
Each discharge phase was divided into a series of segments between which the data was collected. The procedure occurred in 5, 10, 15, 30, 60, and 120 minute segments with four repeats of each (i.e., 4 x 5 min, 4 x 10 min, and so on), which was designed to reflect the relatively rapid rate of morphologic change at the beginning of each phase. For example, in wider channels, alternate bars developed within an hour, and there was relatively little morphologic change in the following hours (Adams and Zampiron, 2020;Adams, 2021a).
Throughout the experiments, sediment falling over the downstream weir was collected in a mesh bucket, drained of excess water, weighed damp to the nearest 0.2 kg, placed on the conveyor belt at the upstream end, and gradually recirculated at the same rate it was output, as opposed to a 'slug' (i.e., all at once) injection. Based on a range of samples collected across the experiments, we determined the weight proportion of water to be approximately 5.8 percent and applied this correction factor to obtain approximate dry weights. There was no initial feed of sediment, although this no-feed period was only 5 minutes.

95
The experiments are best described as pseudo-recirculating as sediment was measured and recirculated at the end of the 5 and 10-minute segments, and for longer segments (i.e. 30, 60, 120 minutes), every 15 minutes.

Data processing
Using the images, point clouds were produced using structure-from-motion photogrammetry in Agisoft MetaShape Professional 1.6.2 at the highest resolution, yielding an average point spacing of around 0.25 mm. Twelve spatially-referenced control 100 points and additional unreferenced ones were distributed throughout the A-BES, which placed photogrammetric reconstructions within a local coordinates system and aided in the photo-alignment process. Using inverse distance weighting, the point clouds were converted to digital elevation models (DEMs) at 1 mm horizontal resolution.
Despite the use of control points, the DEMs contained a slight arch effect in a downstream direction whereby the middle of the model was bowed upwards, which was an artefact of the photogrammetric reconstruction (see doming: James and Robson, 105 2014). This effect was first quantified by applying a quadratic function (on average: p < 0.001, r 2 = 0.999, RMSE = 0.0017 mm) along the length of the bricks, which represent an approximately linear reference elevation (brick elevations vary randomly by ± 2 mm). The arch was then removed by determining correction values along the length of the DEM using the residuals, which were then applied across the width of the model. The final least-squares linear fit along the brick surface was homoscedastic with an average RMSE of 0.0018 mm (around the maximum height difference between adjacent bricks), indicating that the 110 DEM was successfully corrected.
For each DEM, ten wetted cross-sections were reconstructed using the water surface elevation data, which were then used to estimate reach-averaged hydraulics. For more detailed spatial analysis, the flow conditions of water depth, and shear stress were reconstructed using a 2D numerical flow model (Nays2DH) to the final DEM of each discharge phase. The selection of the final DEM was arbitrary as any DEM over the steady-state portion of the experiment could have been selected. Nays2DH is 115 a two-dimensional, depth-averaged, unsteady flow model that solves the Saint-Venant equations of free surface flow with finite differencing based on a general curvilinear coordinate system (further details can be found in Nelson et al., 2016). Notably, local shear stress is calculated using the bed friction coefficient and depth-averaged flow velocity components. Key input boundary conditions are the channel DEM, an initial estimate of reach-averaged Manning's n, cell resolution, and the water discharge. We selected an n value of 0.045 based on the channel conditions, a cell resolution equivalent to 5 mm, and a flow 120 duration of 200 s was sufficient to establish convergence. After an initial model run, we back-calculated a spatially variable value using the flow resistance law presented by Ferguson (2007) that accounts for the influence of relative roughness where f is the Darcy-Weisbach friction factor, a 1 and a 2 are empirically-derived coefficients, d is local flow depth, and the representative roughness length k = D 84 . The spatially variable roughness value was used as an input to run the solver again.

125
To minimise rounding errors associated with the relatively shallow depths in our experiments, the DEM size and discharge were adjusted to the prototype scale (i.e., using a length scale ratio of 25). The estimated water depths, shear stresses and velocities from Nays2DH were then back-transformed to the model scale (Table 3). We removed cells with relatively shallow flows defined arbitrarily as depths less than 2D 84 (6.4 mm) as they contributed a large peak in the frequency distribution of flow depths and likely account for a small proportion of bedload activity. Across the flow models, grid cells with flows less than this 130 threshold accounted for 20-63 percent of the channel area where d > 0, but only 1-21 percent of the total cross-sectional flow area (mean = 11 percent). This is consistent with visual observations of dispersive and stagnant flow at the channel margins.
We defined areas of the bed with flows above the 2D 84 threshold as 'wetted'. The mean-normalised (i.e., local value divided by reach-average) frequency distributions of flow depths and shear stresses were fitted with gamma and Gaussian distributions (coefficients in Table 3), where the goodness-of-fit was assessed using both Kolmogorov-Smirnov and Anderson-Darling tests.

135
The results of the flow model were quantitatively validated by comparing measured reach-averaged hydraulic depths   Table 4.

Determining a representative sediment transport rate
The channels were formed under constant discharge conditions for 4-16 hours, beginning from either a screeded bed or a morphology developed at a lower discharge. Each experimental phase comprised an initial adjustment period during which 145 morphology, hydraulics, and sediment transport were non-stationary. This adjustment period, which varied from minutes to an hour, was followed by a steady-state period where these characteristics fluctuated around a mean value (see Adams, 2020;Adams and Zampiron, 2020). Under recirculating conditions, the stationarity of bedload transport represents a condition in which there is no net aggradation or degradation over time. In Figure 3   transport, followed by fluctuations around a mean value. These fluctuations were explained by processes such as bar reshaping and sediment waves (e.g., Dhont and Ancey, 2018), which are outside the scope of this study.
We determined a representative sediment transport rate for each experimental phase by averaging output over the final threehour period (Table 3), thus removing the initial adjustment period. There was little difference between averaging over the final hour versus the final three hours, with almost all average sediment output values being ± 12.5 percent. There were three 155 instances where these two averaging windows yielded values differing by 15-25 percent due to high-magnitude fluctuations around an otherwise stationary bedload transport rate. In addition, we calculated the standard deviation of sediment output over the final three-hour period.

1D and 2D excess shear stress
We examined the correlation between the observed representative sediment transport rate and two formulations of excess shear stress based on the Meyer-Peter and Müller (1948) equation where q b is width-averaged bedload transport, k accounts for flow resistance and the relative density of sediment, and the exponent 1.6 is based on Wong and Parker (2006). The value of k is highly variable across empirical datasets, whereas the exponent is relatively consistent (Gomez and Church, 1989). The critical shear stress value for the D 50 (τ c50 ) is estimated by where τ * c is the dimensionless critical shear stress, g is gravity, ρ is the density of water, ρ s is the density of sediment We aimed to investigate the concepts underlying 1D and 2D bedload transport equations, rather than refine them. Subsequently, we ignored the parameter k that typically varies across channels and simplified Equation 2 to express the correlation between observed sediment transport and mean excess shear stress (raised to the exponent): This equation was modified to integrate across the distribution of local shear stresses where τ (x) is local bed shear stress and A is the total bed area. Equations 3 and 4 are 1D and 2D approaches to correlating observed transport capacity with excess shear stress. We applied both equations using shear stress values calculated in two 175 ways: 1) the depth-slope product (τ = ρgdS) based on numerically modelled flow depths, and 2) numerically modelled shear stresses, thus yielding four different approaches (Table 4). Each of these approaches are summarised in Appendix A. We intentionally did not account for sinuosity or side-wall effects in the depth-slope product approach. In the case of the 1D depthslope approach, depth was calculated using the mean depth and mean channel gradient, whereas in the 2D depth-slope we varied depth but the gradient remained constant. For each approach, we back-calculated the optimal value of τ * c by systematically 180 varying it and finding the strongest correlation (least-squares linear fit) between q b and excess shear stress (i.e., [τ − τ c50 ] 1.6 or Σ[τ x − τ c50 ] 1.6 /A), indexed by root-mean-square-error (RMSE), which is shown in Figure 4. We report optimised values of τ * c and least-squares goodness-of-fit statistics in Table 4, and also include values obtained using the exponent 1.5 in each equation.  (Table 4).

185
Under the imposed channel widths (0.30-0.60 m) and unit discharges (2.22-7.50 L/m/s) all channels developed an alternate bar morphology with pools, bars, and riffles (see Figure 5 for an example). Especially at low unit discharges, wetted areas (d > 2D 84 ) on average occupied only a portion of the total available width, between 52 and 95 percent. When unit discharge was calculated using the wetted width, it was closely correlated with mean shear stress based on least-squares linear regression (Figure 6a), indicating a coupled adjustment between active width and shear stress. The depth-slope method of calculating mean shear stress estimated higher values compared to the numerical model (7-23 percent higher), and also higher values of critical dimensionless shear stress in the corresponding transport functions (τ * c = 0.066 and 0.050, respectively, Table 4). Both methods yielded similar estimates of excess shear stress (τ / τ c50 = 1.36-2.11 and 1.56-2.53, respectively). The strong correlation between the two estimates of shear stress supports the assumption that at the reach-scaleτ ≈ ρgdS. Estimated values of τ * c using the 2D approaches were consistently higher than the values obtained using the 1D approaches, but were slightly less sensitive to how shear stress was calculated (τ * c ≈ 0.095 for both methods). Based on the 2D approach, the proportion of the wetted bed area experiencing excess shear stress was linearly related to unit discharge and ranged between 37-84 percent (Figure 6b). In several experiments 2D estimates of τ c50 were higher thanτ .
Local shear stresses at or below the mean were estimated to exceed τ c50 only at unit discharges exceeding approximately contributed a second mode of flow depths corresponding to dispersive flow or stagnant water at the channel margins. In the case of the shear stress distributions, the shape parameter α was linearly related to unit discharge based on least-squares regression (RMSE = 0.69, r 2 = 0.39, p < 0.01), and the scale parameter β was negatively correlated (RMSE = 0.58, r 2 = 0.32, p < 0.01).

210
The parameters of the gamma distribution indicate that with increasing unit discharge the distribution of shear stress became more concentrated and less positively skewed.
Despite following similar frequency distributions, modelled local flow depth and shear stress were not strongly coupled spatially (Figure 8). These two parameters were roughly correlated but with considerable scatter, whereby for a given grid cell mean-normalised shear stress was commonly more than a factor-of-two greater or less than normalised flow depth (i.e., high We present the correlation between bedload transport and the four different representations of excess shear stress in Figure 9. These represent combinations of two different methods of calculating bed shear stress, depth-slope product and numerically 220 modelled, against 1D and 2D representations of excess shear stress (Table 4). All four methods yielded similar correlations between excess shear stress and observed bedload transport, indicated by RMSE values between 0.38 and 0.51, where these end-values correspond to the 2D modelled shear stress (B2) and 1D depth-slope product approach (A1), respectively. Changing the exponent from 1.6 to 1.5 in Equations 3 and 4 had almost no effect on the estimated values of τ * c or the prediction errors. Altering the representative grain size from D 50 to D 84 had no effect on the correlaation between q b and excess shear stress (i.e., identical RMSE), and merely reduced the back-calculated estimates of τ * c . especially under conditions where the larger-than-average grain size is only partially mobile (MacKenzie and Eaton, 2017;MacKenzie et al., 2018;Booker and Eaton, 2020;Adams, 2021a). We measured total bedload volumes and adjustments to bed topography during flood stages, which is not possible in the field or in many recirculating experiments. The applied flows were longer and more constant than floods typically observed in nature (4-16 hours experimental time or 20-80 hours in the field 235 prototype), which allowed the experiments to reach an idealised steady-state whereby morphology, hydraulics, and bedload fluctuated around a mean condition (Figure 3). These characteristics make the experimental dataset appropriate for investigat- Note the absence of shallow depths (d < 2D84). The average gamma distribution fit for the normalised shear stress distribution is included (α = 3.30, β = 0.30), as well as the average Gaussian fitted distribution (σ = 0.47). b) Cumulative frequency distribution of non-normalised modelled shear stresses in experimental phases with highest (Exp1c(4)) and lowest (Exp1c(1)) mean shear stress (circled points in Figure 6).
Estimates of τc50 using 1D and 2D approaches (B1 and B2, respectively) are indicated by dashed lines, and the horizontal line is the median shear stress, which closely corresponds to the mean.
ing the effectiveness of bedload transport equations in laterally-constrained gravel-bed rivers under high relative shear stress conditions.
We evaluated four different bedload transport functions based on the correlation between excess shear stress and observed 240 volumes of bedload transport, averaged over the final three hours of each experimental phase. We first focus our discussion on three of these approaches in increasing order of sophistication (A1, B1, then B2), and then explain their relative effectiveness.
Finally, we discuss the conceptual differences between 1D and 2D bedload transport functions.

Comparison between prediction errors
Most bedload transport functions index the applied excess shear stress using the mean depth-slope product as this data is 245 relatively easy to collect in field contexts (Gomez and Church, 1989;Barry et al., 2004;Recking, 2013b). This approach   Table 4. The dashed black line is the least-squares best fit, and solid black lines indicate ± 1 RMSE, and whiskers indicate ± 1 standard deviation over final 3 hours of sediment output measurements.
relies on the assumption that local variations in channel gradient and flow depth cancel out, such that mean flow depth is proportional to mean shear stress (Nicholas, 2000;Ferguson, 2003). We did indeed observe this condition whereby meannormalised flow depth and shear stress followed similar frequency distributions (Figure 7a), despite being decoupled spatially ( Figure 8). The approach A1 (1D depth-slope product) in our analysis was the most simplistic, and in addition, did not account 250 for sinuosity (note the slight sinuosity in Figure 5 that reduced the mean channel gradient), flow resistance, or energy losses to the channel banks. The strength of the correlation between excess shear stress and bedload transport (RMSE = 0.51) provides an approximate reference point for other approaches.
In recent decades, technological advancements in remote sensing and hydraulic modelling have allowed researchers to directly model bed shear stress, thus providing a potentially more accurate estimate. This advancement is utilised in the B1 255 approach (1D modelled shear stress), which accounted for the effect of both sinuosity, flow resistance, and energy losses to the channel banks. Accounting for these additional factors may explain the 13 percent reduction in RMSE (0.44) compared to approach A1. Further advancements have led to the proliferation of 2D hydraulic models and some 2D bedload transport equations, which aim to account for the proportion of the bed participating in transport and the spatial variation in shear stress (Monsalve et al., 2020). The B2 approach (2D modelled shear stress) that integrates across the frequency distribution of shear 260 stresses did not significantly improve upon approach A1, with a similar RMSE (0.38) to approach B1.
Numerical modelling of shear stress and accounting for its frequency distribution led to similarly strong correlations between bedload transport and excess shear stress, compared to the mean depth-slope product method. The ability of the mean shear stress to effectively capture variation in bedload transport is consistent with empirical evidence. In a re-analysis of data from Oak Creek, OR, Monsalve et al. (2020) compared the  equation to a modified 2D version and 265 found that accounting for the distribution of shear stresses reduced prediction error by only 13 percent. Their study modelled a range of flows to the same bathymetry, and we obtained a similar result when the bed was allowed to fully adjust to the imposed flow. Using numerical and analytical models, several studies have predicted that variance in shear stress may enhance bedload transport but that this effect rapidly diminishes whenτ ≫ τ c (Ferguson, 2003;Francalanci et al., 2012;Recking, 2013a). The most probable reason for this sensitivity is the nonlinearity of the bedload transport law, which means that around 270τ ≈ τ c small increases in τ produce relatively large increases in bedload transport. The similar effectiveness of 1D and 2D functions herein provides empirical evidence that bedload transport is less sensitive to the shape of the shear stress distribution under high relative shear stress conditions.

Comparison between 1D and 2D approaches
The four approaches demonstrated key differences based on how shear stress was calculated (depth-slope product vs numeri-275 cally modelled) and more importantly the formulation (1D vs 2D). Both estimates of mean shear stress were linearly related to unit discharge but those based on the depth-slope product were 7-23 percent higher (Figure 6), which is consistent with findings by Monsalve et al. (2020). These differences in estimated shear stress led to approximately commensurate differences in the estimated 1D values of τ * c (32 percent higher). Both 1D estimates of τ * c were relatively high for gravel-bed rivers but were within the range of reported estimates from both field and laboratory channels (Buffington and Montgomery, 1997).
Despite having similar prediction errors, the 1D and 2D functions provided considerably different estimates of critical dimensionless shear stress. Using the 2D approach, estimates of τ * c were 48 and 72 percent higher than the 1D depth-slope and modelled shear stress methods, respectively. In several channels, the estimated critical shear stress was greater than the mean shear stress but bedload transport was observed and well predicted by the model (Figure 6), which in the case of thresholdbased a 1D equation, would correspond to zero estimated transport. This is a distinct advantage of 2D equations at low flows, 285 as they can account for flows where excess shear stress occupies only a fraction of the bed (Monsalve et al., 2020).
The differences between estimates of τ * c arise from differences in how the equations conceptualise excess shear stress. In a 1D equation, when bedload transport data is available, τ c may be back-calculated from the mean shear stress, as is done herein.
The value of τ * c is adjusted until excess shear stress explains the observed bedload transport, assuming thatτ is responsible for all entrainment. In contrast, the 2D equation does not assume that the mean shear stress participates in bedload entrainment.

290
Based on the 2D approach, we estimated that the mean shear stress did not exceed the estimated critical value for the D 50 until a certain discharge threshold (5 L/m/s), and even under the highest flows these areas (i.e., τ c50 < τ <τ ) characterised a maximum of 37 percent of the wetted area. We did not validate these values as we did not anticipate the need for observation, although the estimates appear reasonable compared to our visual observations of the experiments. This result suggests that the mean shear stress is far less significant for bedload transport compared to the larger-than-average stresses, which is intuitive 295 especially given that these are the first stresses to entrain bed material as the flow is increased.
By conceptualising transport as a function of mean shear stress, 1D equations may inflate the importance of relatively moderate shear stresses and deflate values of τ * c . This insight is based on back-calculated values rather than measurements of incipient motion, although it is important to note that studies measuring incipient motion have also been based on the mean shear stress and therefore this 1D paradigm is subsumed within the results (Gilbert, 1914;Kramer, 1935;Neill and 300 Yalin, 1969;Wilcock, 1988). We also relied on spatio-temporally integrated rather than instantaneous local shear stresses that promote entrainment (e.g., Nelson et al., 1995). Nevertheless, the higher estimates of critical dimensionless shear stress using the 2D approach, evaluated by considering the relative importance of shear stresses across the frequency distribution, may have a stronger conceptual basis. More broadly, the results highlight that as long as τ c is back-calculated, its value will be highly dependent on how shear stress is estimated, and whether its distribution is treated one-or two-dimensionally.

305
The results may have implications for non-threshold approaches to predicting bedload transport in natural gravel-bed rivers Parker, 1990;Wilcock and Crowe, 2003;Recking, 2013a). These approaches recognise that usage of a single critical shear stress is ineffective at low flows and is always an approximation, especially in the case of partial transport conditions where not all grain sizes (or even grains of a given size) are equally mobile (Wilcock and McArdell, 1993). The effectiveness of threshold-based approaches under high excess shear stresses suggests that in channels with fully-developed 310 morphology and a wide range of grain sizes, non-threshold-based approaches may not render an improvement. Also, the results challenge recent critiques of bedload transport predictions based on mean shear stress, and particularly the depth-slope assumption (Yager et al., 2018). Although there is a poor mechanistic link between shear stress and bedload transport (e.g., Nelson et al., 1995), this approximation may be unreasonably effective when applied at a sufficiently large spatio-temporal scale or excess shear stress.
Further work is required to investigate differences in 1D and 2D estimates of τ * c under lower excess shear stress conditions. If broadly applicable, the effectiveness of highly reductionist bedload transport functions based only on median grain size and mean shear stress would present a convenient assumption for researchers and practitioners interested in channel-forming flows. More research is required to substantiate this approach under supply-limited conditions and realistic hydrographs that enable both upward and downward adjustments with inherited channel conditions. Given that our experiments do not allow 320 for significant lateral adjustment and meandering, the results are most applicable to channels confined by bedrock, or with cohesive or highly vegetated banks. Fully alluvial channels comprise additional feedbacks that are worthy of investigation, and the extent to which these affect reach-averaged bedload transport remains poorly understood.

Conclusions
We investigated the performance of 1D and 2D bedload transport functions under high relative shear stress conditions in a 325 Froude-scaled physical model. The analysis highlights the effectiveness of highly reductionist bedload transport functions based only on median grain size and mean shear stress calculated using the depth-slope product. Numerically modelling shear stress to account for flow resistance and energy losses from the channel planform and banks did not substantially reduce prediction error, nor did accounting for the relative importance of shear stresses across the frequency distribution. The results suggest that bedload transport may collapse to a more simple function (i.e., with average shear stress and grain size) under high 330 excess shear stress conditions. Given the channels herein have limited lateral mobility, our conclusions are most applicable to channels where lateral adjustment is suppressed. Further work is required to examine the effect of planform adjustments (widening, meandering), where small-scale laboratory experiments serve as an effective research tool. The 1D and 2D approaches provided substantially different estimates of critical dimensionless shear stress, reflecting differences in how these approaches conceptualise excess shear stress. Estimates of τ * c from 2D functions may have a stronger conceptual basis, as they 335 are derived by considering the relative importance of shear stresses across the frequency distribution, and do not assume that the mean shear stress is sufficient to mobilise the median grain size.
Code and data availability. Raw hydraulic and sediment transport data from Tables 3-4  Approach A1: mean shear stress assuming depth-slope product where q b is unit bedload transport,τ is estimated with ρgdS (where g is gravity, d is mean flow depth from a 2D flow model, S is the channel gradient, and τ c50 is estimated by τ * c g(ρ s − ρ)D, where τ * c is the dimensionless critical shear stress chosen using best-fit (Table 4), ρ is the density of water, ρ s is the density of sediment.

345
Approach A2: distribution of shear stresses assuming depth-slope product where τ (x) is the local bed shear stress calculated using local depth from a 2D flow model (ρgd (x) S) and A is the total active bed area (defined as areas where d > 2D 84 ).

350
Approach B1: mean shear stress based on numerical model whereτ is estimated using a 2D flow model.

355
Approach B2: distribution of shear stresses based on numerical model where τ (x) values are based on 2D flow model.
Author contributions. DA was responsible for conceptualisation, data collection, formal analysis, visualisation, and writing. BE was responsible for supervision, review, and editing.