Modeling memory in gravel-bed rivers: a flow-history-dependent relation for evolving thresholds of motion

Thresholds of motion (${\mathit{\tau }}_{\mathrm{c}}^{\ast }$) strongly control bedload transport in gravel-bed rivers. Uncertainty in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ limits the accuracy of predictions of transport and morphologic change. To improve our quantitative understanding of morphodynamic feedbacks in rivers, we propose a flow-history-dependent model where ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ evolves temporally as a function of bed shear stress. Relatively low shear stresses strengthen the bed, increasing ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ and reducing transport. Larger floods rapidly weaken the bed, decreasing ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ and increasing transport. We calibrate the model to a 23-year record of flow and bedload transport from the Erlenbach torrent, Switzerland, and find that the model predicts the field-based ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ record more accurately than assuming a constant ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$. Calibrated parameters describing strengthening are more tightly distributed than weakening parameters, which suggests that magnitudes of bed weakening may be more variable and difficult to accurately predict as a function of flood characteristics than bed strengthening during lower flows.

1 Introduction

Erosion, deposition, and morphological change in gravel-bed rivers result from bedload transport. Gravel transport rates are notoriously difficult to accurately predict in natural rivers because they are influenced by a wide variety of factors including water discharge, local channel morphology, and upstream sediment supply. These factors can vary both spatially and temporally. Nonetheless, most equations to predict bedload flux simplify this complexity using a deterministic parameterization which only considers flow intensity and a transport threshold. For example, the classic Meyer-Peter and Müller (1948) equation (MPM) can be expressed simply as $q_{\mathrm{s}}^{\ast }=\mathrm{4}{\left({\mathit{\tau }}^{\ast }-{\mathit{\tau }}_{\mathrm{c}}^{\ast }\right)}^{\mathrm{1.5}}$ for ${\mathit{\tau }}^{\ast }\ge {\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (Wong and Parker, 2006), where $q_{\mathrm{s}}^{\ast }$ is dimensionless sediment flux per unit channel width (Einstein number), τ^∗ is nondimensional shear stress (Shields stress), and threshold parameter ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ is the critical Shields stress. Shields stress is ${\mathit{\tau }}^{\ast }=\mathit{\tau }/\left({\mathit{\rho }}_{\mathrm{s}}-\mathit{\rho }\right)gD$, where τ, ρ_s,ρ, g, and D are dimensional bed shear stress (Pa), sediment density (kg m⁻³), water density (kg m⁻³), gravitational acceleration (m s⁻²), and median sediment diameter (m). Einstein number is defined as $q_{\mathrm{s}}^{\ast }=q_{\mathrm{s}}/\sqrt{\left({\mathit{\rho }}_{\mathrm{s}}/\mathit{\rho }-\mathrm{1}\right)g{D}^{\mathrm{3}}}$, where q_s is the volumetric bedload transport rate per unit width (m³ s⁻¹ m⁻¹). While the threshold parameter (${\mathit{\tau }}_{\mathrm{c}}^{\ast })$ describes the nondimensional shear stress at the onset of sediment motion, it is also often used as a physically meaningful, but empirically determined, fitting parameter in bedload transport models (e.g., Engelund and Fredsoe, 1976; Luque and Van Beek, 1976; Meyer-Peter and Müller, 1948; Wong and Parker, 2006; Shields, 1936; Wiberg and Smith, 1987; Wilcock and Crowe, 2003). For bedload equations such as the MPM, ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ is a lumped parameter that implicitly accounts for alłof the factors that influence bedload transport rates apart from τ^∗.

Even during large floods, bedload transport often occurs at shear stresses only slightly exceeding threshold conditions, making transport rate predictions particularly sensitive to threshold values (e.g., Parker, 1978; Phillips and Jerolmack, 2016; Phillips et al., 2022; Pretzlav et al., 2020). Transport thresholds also strongly influence gravel-bed channel reach morphologies (e.g., Parker, 1978; Phillips et al., 2022) and modulate the mapping of climate onto fluvial processes, informing short- and long-term sediment fluxes and the relative importance (or unimportance) of extreme events for channel evolution (e.g., Blom et al., 2017; DiBiase and Whipple, 2011; Lague et al., 2005; Shobe et al., 2018; Tucker and Bras, 2000). For these reasons, improving our ability to calculate thresholds of motion is critical not only for predicting transport rates but also for predicting mountain river morphodynamics and channel stability.

Early work assumed that thresholds were primarily controlled by grain weight relative to the fluid. By accounting for these variables through nondimensionalization, ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ was initially thought to be approximately constant for typical conditions in gravel-bed rivers (e.g., Buffington and Montgomery, 1997; Shields, 1936 ). However, Buffington and Montgomery (1997) showed that ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ varied systematically with the ratio of the median grain size to flow depth, independently of shear stress. More recent work has explored how both flow and grain interactions lead to inherent variability in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$. Thresholds vary spatially with the surrounding grain size distribution (Parker, 1990; Chen and Stone, 2008), reach slope (Lamb et al., 2008; Mueller et al., 2005), bed morphology (Dietrich et al., 1989; Gran et al., 2006; Powell et al., 2016; Monsalve and Yager, 2017; Powell and Ashworth, 1995; Roberts et al., 2020), and changes in riverbed microtopography (Brayshaw, 1985; Hodge et al., 2019; Kirchner et al., 1990; Masteller and Finnegan, 2017; Yager et al., 2018).

Thresholds for motion also evolve over time. For example, hysteresis in bedload transport rates is often observed between the rising and falling limbs of individual floods (Hsu et al., 2011; Mao, 2018; Mao et al., 2014; Pretzlav et al., 2020; Reid et al., 1985; Roth et al., 2017). Predicting hysteresis using the MPM and similar bedload models requires the threshold parameter to evolve over the course of a flood event (assuming that the prefactor and exponent remain constant). Changes in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ over multiple events have also been observed, and, in most cases, ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ values remain correlated across events, indicating a memory of past conditions (Downs and Soar, 2021; Hassan et al., 2020; Johnson, 2016; Lenzi et al., 2004; Mao, 2018; Masteller et al., 2019; Rickenmann, 2018, 2020; Saletti et al., 2015; Turowski et al., 2011).

Variable flow strength influences threshold evolution through time. Reid et al. (1985) first suggested the influence of antecedent flows based on field-based bedload transport monitoring, hypothesizing that longer inter-flood durations led to increases in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ and reduced sediment transport rates. Experiments have confirmed that the magnitude of inter-event flow affects ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ evolution (Haynes and Pender, 2005; Masteller and Finnegan, 2017; Monteith and Pender, 2005; Ockelford et al., 2019; Ockelford and Haynes, 2013; Paphitis and Collins, 2005). With little to no active sediment transport, grain-scale changes in interlocking and surface reorganization increase particle resistance to motion (Masteller and Finnegan, 2017; Ockelford and Haynes, 2013; Yager et al., 2018). Pretzlav et al. (2020) documented systematic discharge-dependent increases and decreases in motion thresholds and associated diurnal transport hysteresis during several weeks of snowmelt flooding, using instrumented “smartrocks” to measure transport. Reduction in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ following larger floods has been attributed to significant reorganization of the riverbed (Lenzi et al., 2004; Turowski et al., 2009).

Sediment supply also influences threshold evolution. Hysteresis can be caused by sediment supply variations through time (Moog and Whiting, 1998; Mao et al., 2014), changing thresholds. Increased sediment supply from channel banks and hillslopes can be important in destabilizing the bed surface or introducing mobile unconsolidated material, reducing thresholds of motion (Turowski et al., 2011; Recking et al., 2012; Rickenmann, 2020). Building on observations by Recking et al. (2012), Johnson (2016) developed a model in which ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ evolves as a function of net erosion or deposition, which are controlled by sediment supply in relation to transport capacity. After calibration to laboratory experiments, the evolving ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ model successfully predicted how transport rates responded to pulses in sediment supply.

Decades of monitoring data from the Erlenbach torrent in Switzerland similarly provide evidence for transport thresholds evolving with both sediment supply and discharge variability. Rickenmann (2020) showed that variations in sediment availability on the bed correlated with sediment transport rate fluctuations and evolving thresholds, suggesting that thresholds depend on upstream sediment supply. For the same stream, Masteller et al. (2019) showed that the magnitude of antecedent flows also influenced the evolution of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ for individual years. Consistent with experiments, Masteller et al. (2019) observed that the start of transport events showed increases in critical Shields stress with increasing inter-event flow magnitude (herein termed “strengthening”) for an intermediate range of flows spanning inter-event periods and floods with observable sediment transport. However, following even higher magnitude flows, the threshold for motion decreased (herein termed “weakening”). Masteller et al. (2019) hypothesized that the transition from bed strengthening to bed weakening was associated with a transition from local rearrangement of particles to more intense transport disrupting bed structure via particle collisions, and/or enhanced upstream sediment supply through upstream bed erosion (Yager et al., 2012), and/or enhanced hillslope–channel coupling (Golly et al., 2017). Thus, both flow strength and sediment supply likely influence thresholds of motion in the Erlenbach torrent (Rickenmann, 2020; Masteller et al., 2019; Turowski et al., 2011).

The ability to accurately predict threshold evolution through time – thereby improving bedload transport rate predictions – remains elusive due to a lack of validated models. Equations have been proposed to describe temporal bed strengthening as a function of the duration of bed exposure to a constant inter-event flow magnitude and an initial ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ based on experimental data (e.g., Ockelford et al., 2019; Paphitis and Collins, 2005). However, because these models only focus on inter-event strengthening effects, they cannot capture decreases in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$. The model of Johnson (2016) predicts ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ evolution as a function of changing sediment supply. Nonetheless, this model is an incomplete description of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ evolution because it does not account for riverbed strengthening or weakening directly caused by the flow. Notably, to our knowledge, none of these equations have been used to describe field observations of temporally varying ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$.

Our goals in the present work are (i) to propose a new model in which ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ evolves as a function of flow magnitude and encapsulates some memory of past shear stresses as reflected in the changing state of the riverbed and (ii) to evaluate whether this discharge-dependent model can capture annual strengthening and weakening trends observed in Erlenbach field data (Masteller et al., 2019). While sediment supply variations can also influence threshold evolution in general (Johnson, 2016), and in the Erlenbach in particular (Rickenmann, 2020), we explore how well threshold evolution can be predicted using only a time series of river discharge.

2 Model development

Johnson (2016) argued that ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ is a “state variable” for gravel-bed river morphodynamics because it simultaneously controls transport rates and evolves due to feedbacks with fluid shear stresses and transport rates. Our new equations take a similar form to Johnson (2016), where changes in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ depend not only on discharge-dependent shear stress but also on the current state of the transport system as characterized by ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ itself. The rate of change in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ depends on two terms, which both evolve as a function of the transport capacity, ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }$. Conceptually, the first right-hand-side term (starting with k₁) represents strengthening processes that increase ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$, while the second term (starting with k₂) represents weakening processes that reduce ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\displaystyle \frac{\partial {\mathit{\tau }}_{\mathrm{c}}^{\ast }}{\partial t}}=k_{\mathrm{1}}B{\left(\mathrm{1}+{\left({\displaystyle \frac{{\mathit{\tau }}^{\ast }}{{\mathit{\tau }}_{\mathrm{c}}^{\ast }}}\right)}^{-\mathit{\gamma }}\right)}^{-\mathrm{1}}-k_{\mathrm{2}}\left(\mathrm{1}-B\right){\left({\displaystyle \frac{{\mathit{\tau }}^{\ast }}{{\mathit{\tau }}_{\mathrm{c}}^{\ast }}}-\mathrm{1}\right)}^{\mathit{\epsilon }}\\ & H\left[{\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }-\mathrm{1}\right],\end{array}\end{array}$$

where t is time, and

$$\begin{array}{}\text{(2)}& B={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{cmax}}^{\ast }-{\mathit{\tau }}_{\mathrm{c}}^{\ast }}{{\mathit{\tau }}_{\mathrm{cmax}}^{\ast }-{\mathit{\tau }}_{\mathrm{cmin}}^{\ast }}}\text{for}{\mathit{\tau }}_{\mathrm{cmin}}^{\ast }<{\mathit{\tau }}_{\mathrm{c}}^{\ast }<{\mathit{\tau }}_{\mathrm{cmax}}^{\ast }.\end{array}$$

The scaling factor k₁ and exponent γ influence the form and magnitude of the strengthening term, while k₂ and ε do the same for the weakening term (Fig. 1); these parameters are empirically calibrated below. Both k₁ and k₂ have units of 1/time, where the units of time will depend on the time step of the discharge time series used as model input. H is the Heaviside step function ($H[{\mathit{\tau }}^{*}/{\mathit{\tau }}_{\mathrm{c}}^{*}-\mathrm{1}]=\mathrm{0}$ for ${\mathit{\tau }}^{*}/{\mathit{\tau }}_{\mathrm{c}}^{*}<\mathrm{1}$; $H[{\mathit{\tau }}^{*}/{\mathit{\tau }}_{\mathrm{c}}^{*}-\mathrm{1}]=\mathrm{1}$ for ${\mathit{\tau }}^{*}/{\mathit{\tau }}_{\mathrm{c}}^{*}>=\mathrm{1}$) such that weakening only occurs when transport occurs (${\mathit{\tau }}^{\ast }>{\mathit{\tau }}_{\mathrm{c}}^{\ast }$) (Fig. 1a). B is called the “feedback parameter,” described below. ${\mathit{\tau }}_{\mathrm{cmin}}^{\ast }$ and ${\mathit{\tau }}_{\mathrm{cmax}}^{\ast }$ are upper and lower bounds imposed on ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$, representing physical limits for how loosely packed and mobile, or tightly packed and immobile, the bed surface can become (Johnson, 2016).

The strengthening and weakening terms combine to cause increases and decreases in $\partial {\mathit{\tau }}_{\mathrm{c}}^{\ast }/\partial t$ (Eq. 1; Fig. 1). The strengthening term is generally sigmoidal for γ>1; it goes to zero as τ^∗ approaches zero, and asymptotes to a value of k₁B for ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }\gg \mathrm{1}$ (Fig. 1a). We chose a sigmoidal form to allow strengthening over a wide range of flows but also limit the amount of incremental strengthening that can result from changes in grain organization at higher transport capacities. When ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }<\mathrm{1}$, flow causes the bed to become stronger but not weaker, consistent with previous observations (Haynes and Pender, 2005; Masteller et al., 2019; Masteller and Finnegan, 2017; Monteith and Pender, 2005; Ockelford et al., 2019). Strengthening increases as ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }$ approaches 1, consistent with some (Paphitis and Collins, 2005) but not all previous work (Haynes and Pender, 2007). Strengthening increases further for ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }>\mathrm{1}$, consistent with protrusion-dependent thresholds (Masteller and Finnegan, 2017, Yager et al., 2018; Masteller et al., 2019) and with coarse grain clustering, which increases bed stability and requires transport to develop (Brayshaw, 1985; Church et al., 1998; Hassan et al., 2020; Johnson, 2017; Strom et al., 2004).

At the same time, as ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }$ exceeds 1, the weakening term becomes increasingly important (Eq. 1; Fig. 1). In the absence of other constraints on the functional form of weakening with increasing ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }$, we chose a power-law relation for simplicity. It seems likely to us that beds rapidly lose their strength as transport rate increases and fewer grains are interlocked through intergranular friction (Yager et al., 2018). Higher shear stresses capable of mobilizing more sediment grains can destabilize a larger fraction of the bed. Impacts from transported grains may also directly contribute to destabilization (Ancey and Heyman, 2014; Heyman et al., 2014; Lee and Jerolmack, 2018; Martin et al., 2014). Nonetheless, we note that Eq. (1) is agnostic towards any specific processes driving strengthening and weakening. The combination of terms results in the transition from strengthening to weakening occurring at different ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }$, depending on γ, ε, k₁, k₂, and ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (Fig. 1b).

B is considered a “feedback parameter” because it contributes to $\partial {\mathit{\tau }}_{\mathrm{c}}^{\ast }/\partial t$ being a function of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (Johnson, 2016). B has a value between 0 and 1, and changes the importance of the strengthening and weakening terms, depending on the current value of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ relative to ${\mathit{\tau }}_{\mathrm{cmin}}^{\ast }$ and ${\mathit{\tau }}_{\mathrm{cmax}}^{\ast }$. A loosely packed bed, with ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ close to ${\mathit{\tau }}_{\mathrm{cmin}}^{\ast }$ and B close to 1, can strengthen significantly in a low flow (increasing ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$), but a high flow would not cause a significant decrease in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ because the bed is already relatively weak. Conversely, a bed that was nearly as tightly packed as physically possible, with ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ close to ${\mathit{\tau }}_{\mathrm{cmax}}^{\ast }$ and B close to 0, will minimally increase ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ in response to a low flow, but a destabilizing flood would cause a significant decrease in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (Johnson, 2016).

To estimate the range of possible ${\mathit{\tau }}_{\mathrm{cmin}}^{\ast }$ and ${\mathit{\tau }}_{\mathrm{cmax}}^{\ast }$ values, we use a compilation of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ as a function of channel reach slope (Johnson, 2016; Lamb et al., 2008; Prancevic and Lamb, 2015). A best-fit power-law regression of the compiled data gives ${\mathit{\tau }}_{\mathrm{c}}^{\ast }=\mathrm{0.42}{S}^{\mathrm{0.7}}+\mathrm{0.03}$, with R²=0.41 where S is channel reach slope (an approximation of energy or water surface slope). The constant is included so that the function asymptotes to a physically reasonable value (i.e., ${\mathit{\tau }}_{\mathrm{c}}^{\ast }=\mathrm{0.03}$) as slope approaches zero. Field and flume data were weighted equally in the nonlinear best-fit regression, removing possible bias from there being ∼3.5 times more flume data points. We then visually determined the following minimum and maximum bounds to accommodate almost all compilation data points, assuming the same best-fit exponent (0.7).

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{cmin}}^{\ast }=\mathrm{0.14}{S}^{\mathrm{0.7}}+\mathrm{0.0075}\text{(4)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{cmax}}^{\ast }=\mathrm{1.4}{S}^{\mathrm{0.7}}+\mathrm{0.075}\end{array}$$

These empirical ${\mathit{\tau }}_{\mathrm{cmin}}^{\ast }$ and ${\mathit{\tau }}_{\mathrm{cmax}}^{\ast }$ relations capture the slope dependence of gravel thresholds of motion compiled in both flume and field settings. Equations (1) and (2) assume that ${\mathit{\tau }}_{\mathrm{cmin}}^{\ast }$ and ${\mathit{\tau }}_{\mathrm{cmax}}^{\ast }$ are limits that the threshold can approach but does not reach due to the feedbacks implemented by the B parameter. Equations (1) and (2) are only defined between these bounds. We acknowledge that many factors beyond slope may influence the possible range of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ at a given site, including grain size distributions (e.g., Parker, 1990) and differences in relative roughness (Buffington and Montgomery, 1997; Schneider et al., 2015). In the absence of independent constraints on ${\mathit{\tau }}_{\mathrm{cmin}}^{\ast }$ and ${\mathit{\tau }}_{\mathrm{cmax}}^{\ast }$ that could be used to describe a particular field site, the compilation should reasonably represent the range of possible values.

Figure 1 (a) Predicted change in threshold for strengthening (blue) and weakening (red) terms for a range of γ and ε values. k₁, k₂, and B are constant. (b) Predicted change in threshold for the full model for different B and ε. k₁, k₂, and γ are held constant. This example uses Erlenbach values calculated as ${\mathit{\tau }}_{\mathrm{cmin}}^{\ast }=\mathrm{0.036}$ and ${\mathit{\tau }}_{\mathrm{cmax}}^{\ast }=\mathrm{0.36}$ using Eqs. (3) and (4).

Figure 2 (a) Compiled ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ data as a function of slope from both field studies (gray diamonds) and flume experiments (open circles). The range of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ from the Erlenbach field site is plotted in red. Compiled data are limited to slopes S<0.2 and median grain sizes D₅₀≥2 mm. Most data were compiled by Prancevic and Lamb (2015), building on Buffington and Montgomery (1997), with additional data from Olinde (2015) and Lenzi et al. (2006). (b) Calibration of strengthening term exponent, γ, based on Paphitis and Collins (2005) for both the shortest conditioning time (5 min) and the longest conditioning time (120 min) spanned by the Paphitis and Collins (2005) data. Regressions to their best-fit empirical equation give gamma exponents within uncertainty of each other.

We calibrated γ using experimental data from Paphitis and Collins (2005) to reduce the number of free parameters in the model (i.e., γ, ε, k₁, k₂). Paphitis and Collins (2005) conducted experiments using fine, medium, and coarse sand, in which they systematically varied the conditioning time (E_D, the duration of flow below the threshold condition) between 5 and 120 min and the ratio of shear velocity (u_τ) to initial critical shear velocity (u_τci). The ratios of shear velocities they explored were between 70 % and 95 % of critical, which corresponds to conditioning flow shear stresses between 50 % and 90 % of the initial critical shear stress. Based on their experimental results, they presented the following equation to describe their data:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\displaystyle \frac{u_{\mathit{\tau }\mathrm{c}\left(t\right)}}{u_{\mathit{\tau }\mathrm{ci}}}}=\mathrm{1.05}\left[\mathrm{1}-\mathrm{0.01}{e}^{\left(-\mathrm{0.005}{E}_{\mathrm{D}}\right)}\right]\\ & +\left[\mathrm{0.005}+\mathrm{0.1}\left({\displaystyle \frac{u_{\mathit{\tau }}}{u_{\mathit{\tau }\mathrm{ci}}}}-\mathrm{0.7}\right)\right]\ln \left(E_{\mathrm{D}}\right)\\ & +\mathrm{0.06}\left[{\mathrm{10}}^{-\mathrm{7}\left(\mathrm{0.97}-\frac{u_{\mathit{\tau }}}{u_{\mathit{\tau }\mathrm{ci}}}\right)}\right]\end{array}\end{array}$$

for $\mathrm{0.7}\le \frac{u_{\mathit{\tau }}}{u_{\mathit{\tau }\mathrm{ci}}}\le \mathrm{0.95}$ and E_D≤120 min, where u_τc(t) is the critical shear velocity following low-flow conditioning. Paphitis and Collins (2005) reported that this function fit their experimental data with a correlation coefficient of 0.83 (i.e., R²=0.69).

To calibrate γ, we calculate $u_{\mathit{\tau }\mathrm{c}\left(t\right)}/u_{\mathit{\tau }\mathrm{ci}}$, for a range of E_D and u_τ, using values of u_τci reported for their coarse sand (D=0.774 mm) experiments, corresponding to an initial critical shear velocity of u_τci=0.0195 m s⁻¹. We then square $u_{\mathit{\tau }\mathrm{c}\left(t\right)}/u_{\mathit{\tau }\mathrm{ci}}$ to convert to a Shields stress ratio (i.e., ${\left(u_{\mathit{\tau }\mathrm{c}\left(t\right)}/u_{\mathit{\tau }\mathrm{ci}}\right)}^{\mathrm{2}}={\mathit{\tau }}_{\mathrm{c}}/{\mathit{\tau }}_{\mathrm{ci}}={\mathit{\tau }}_{\mathrm{c}}^{\ast }/{\mathit{\tau }}_{\mathrm{ci}}^{\ast }$, where the i subscript indicates the initial value and τ_c and ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ evolve through time). We then numerically calculate the partial derivative of Eq. (5) with respect to time, $\partial {\mathit{\tau }}_{\mathrm{c}}^{\ast }/\partial t$. Figure 2b shows a nonlinear regression (using MATLAB's cftool) of the strengthening term in our model (Eq. 1) to $\partial {\mathit{\tau }}_{\mathrm{c}}^{\ast }/\partial t$ calculated from Eq. (5). This regression provides a best-fit estimate of γ, including empirical regression uncertainties. For below-threshold conditions explored in the Paphitis and Collins (2005) experiments, the weakening term in Eq. (1) is zero. Given this, the calibration of γ described here is not influenced by other model parameters, particularly the weakening exponent ε. Our reported 95 % confidence interval on γ only represents the empirical regression uncertainty when fitting our function to Eq. (5). Therefore, it is likely that a somewhat wider range of γ may be able to fit the range of the experimental data from Paphitis and Collins (2005), and the true range of possible values may be somewhat larger than 2.5±0.32.

We use the data and fitting function of Paphitis and Collins (2005) to calibrate γ because it is the most complete and internally consistent dataset that we are aware of with sufficient constraints to describe the evolution of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ as a function of both transport capacity and time. Nonetheless, a possible limitation of applying these experimental data to calibrate our model is that the Paphitis and Collins (2005) experiments were conducted with unimodal sand. Specifically, boundary Reynolds numbers in their experiments are transitional between hydraulically smooth and hydraulically rough flow. For the coarsest grains they use (D₅₀=0.0774 mm), the boundary Reynolds number $R{e}_w=u_{\mathit{\tau }}{k}_{\mathrm{s}}/\mathit{\upsilon }\approx \mathrm{15}$ , where k_s is a roughness length scale assumed to be D₅₀ and υ is the kinematic viscosity of water. If k_s were instead assumed to be a multiple of D₅₀ (such as k_s=3.5D₈₄), then Re_w would be closer to the hydraulically rough flow criteria of Re_w≥100. It is also worth noting that grain size did not explicitly factor into Eq. (5) beyond its implicit control on u_τci. The insensitivity of their results to grain size suggests that the results may not depend significantly on grain size or on hydraulically rough flow being fully developed. Furthermore, converting critical shear velocities from their coarse sand experiments to critical Shields stress yields ${\mathit{\tau }}_{\mathrm{c}}^{\ast }\approx \mathrm{0.031}$, consistent with typical critical Shields stresses reported at low slopes (consistent with their experiments) for gravel river datasets (Fig. 2a).

3 Field application

The Erlenbach is a small (0.7 km² watershed), steep (10 % grade) channel in the Swiss Prealps. Bedload transport has been actively monitored for over 30 years by a variety of methods (Beer et al., 2015; Rickenmann, 2020; Rickenmann et al., 2012; Rickenmann and McArdell, 2007). Previous analyses have shown that the threshold for motion varies over 1 order of magnitude across the span of the record (Masteller et al., 2019; Turowski et al., 2011). Masteller et al. (2019) demonstrated that seasonal trends in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ were unlikely to be random and that threshold evolution depends, in part, on the magnitude of past flows.

Our goal is to evaluate how well discharge-dependent Shields stress variations alone (Eqs. 1 and 2) capture first-order seasonal trends in evolving ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ from well-constrained field data. We utilize publicly available 10 min interval discharge and bedload transport records from the Erlenbach over the last 23 years (Rickenmann et al., 2020). To calculate thresholds of motion for model comparison, we measure the discharge when the in-channel impact plate system registers grain collisions near the beginning of a flow event, following Turowski et al. (2011). Thus, the threshold data are independent of any bedload transport model. Using a rating curve from discharge to shear stress developed by Yager (2006) and a median grain size, D₅₀=8 cm (Wyss et al., 2016), the flow and transport time series were nondimensionalized to Shields stress. The median grain size is assumed not to have changed systematically across the Erlenbach record (as discussed by Masteller et al., 2019), although hillslope sediment supply may cause grain size variability both during and in between transporting events.

Critical Shields stresses at the start of transport vary by almost 1 order of magnitude in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (0.03 to 0.26). The temporal variability in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ observed at the Erlenbach is equivalent to the full range of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ observed in flume and field data for equivalent slopes (Fig. 2a). Strengthening (i.e., a systematic increase in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ for at least some portion of a given year) is dominant in 10 of the 23 years (see Masteller et al., 2019) and weakening is dominant in 3 years (1992, 2014, 2015), while the remaining 10 years show both behaviors (Fig. 3e–h).

4 Model parameterization and application

We implement the model separately to each year's flow time series, from the first transporting event in the spring through the fall (following Masteller et al., 2019). We do not calibrate ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ to the single continuous multi-year discharge and transport record because the bulk of landsliding occurs during the winter months, supplying largely unconstrained amounts of sediment to the channel bed from hillslope processes (Schuerch et al., 2006). Hillslope sediment supply variations also occur during the rest of the year and likely influence thresholds and transport rates both during and in between the transporting events we consider (e.g., Rickenmann, 2020, 2024). As possible evidence of sediment supply effects during inter-event periods, Turowski et al. (2011) found that threshold discharges were often, though not always, lower at the start of a given flow event compared to the discharge at the cessation of bedload transport at the end of the previous event. This inter-event weakening cannot be captured by the model. Because the model cannot predict every trend in the field data, we focus on the start of events only to evaluate how well discharge variations alone can improve transport predictions over seasons, consistent with the analysis of Masteller et al. (2019). Future analyses could focus on threshold evaluation and model calibration during individual flood events.

Equation (1) has four free parameters: k₁, k₂, γ, and ε. We assign γ=2.5 based on our calibration to Paphitis and Collins (2005) (Fig. 2b). This leaves three parameters that require calibration. For each year, we explored a range of parameter combinations for k₁ ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ to $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ min⁻¹, n=40), k₂ ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ to $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ min⁻¹, n=40), and ε (1 to 10, n=10). k₁ and k₂ were varied with log-spacing to explore all orders of magnitude equivalently. For each year of the dataset, we ran 16 000 forward simulations, reflecting all unique parameter combinations of k₁, k₂, and ε. For each year, we assign an initial ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ value as equal to the observed ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ at the first transport event and calculate changes in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ based on 10 min discharge data until the end of the final observed transport event.

We determined the best-fit parameter combinations that minimized mean absolute error (MAE) to ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ data for each year. We use MAE (rather than RMSE) to reduce the influence of any single large difference between the continuous model predictions and field-based data points which exhibit large amounts of scatter. The field data points only represent the discrete start of each event. “Annual” calibrations represent the best-fit parameters for each year. The “combined” calibration represents the single best-fit k₁, k₂, and ε values which minimize MAE when MAE is averaged across all sample years, with each year weighted equally. We compare both the annual and combined best-fit model to a constant ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (mean Erlenbach ${\mathit{\tau }}_{\mathrm{c}}^{\ast }=\mathrm{0.1548}$, SE = 0.0014).

5 Results

The annual calibrations show that, for all 23 years, Eq. (1) provides a better fit to the data (lower MAE) than the mean Erlenbach threshold, ${\mathit{\tau }}_{\mathrm{c}}^{\ast }=\mathrm{0.1547}$ (Fig. 3a). Annual best-fit MAE ranges from 0.0046 (2015) to 0.0293 (1990). The mean MAE = 0.0149 from the annual calibrations is less than the mean MAE = 0.0254 when applying a constant ${\mathit{\tau }}_{\mathrm{c}}^{\ast }=\mathrm{0.1547}$. Annual best-fit values for k₁ ranged from $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ to $\mathrm{2.73}\times {\mathrm{10}}^{-\mathrm{5}}$ min⁻¹ with a mean value of $\mathrm{6.55}\times {\mathrm{10}}^{-\mathrm{6}}$ min⁻¹ (Fig. 3b). In contrast, annual best-fit k₂ values spanned the entire parameter range ($\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ to 10⁻² min⁻¹) with mean $k_{\mathrm{2}}=\mathrm{2.63}\times {\mathrm{10}}^{-\mathrm{4}}$ min⁻¹ (Fig. 3c). Best-fit annual ε values also spanned the full range of parameter values explored (1 to 10) (Fig. 2c), with mean ε=6 when calculated with each year weighted equally.

Figure 3 (a) Comparison of MAE fits for constant and evolving ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ values compared to ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ data from the Erlenbach for both annually calibrated and combined calibration models. (b–d) Median and interquartile range of MAE for model runs binned by (b) k₁, (c) k₂, and (d) ε parameters. All annual minimum MAE values and associated parameters are indicated by gray stars. In the event of multiple parameter combinations yielding the minimum MAE (see Discussion), all parameter combinations are plotted. (e–h) Comparison of field data (black dots represent ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ at the start of each transporting event) to the best-fit model(s). Example years with (e) strengthening only (1988); (f) dominantly weakening, with strengthening as well (2007); (g) dominantly strengthening, with weakening as well (2010); and (h) dominantly weakening (2014). Upper panels show the flow time series for each year in blue; the model parameter ${\mathit{\tau }}_{\mathrm{cmax}}^{\ast }=\mathrm{0.36}$ is also indicated. The feedback parameter B is shown for each year (right-hand y axis). “n” gives the number of models (i.e., parameter combinations) that minimize MAE and provide an equivalent best fit to the data. Annual best-fit parameters are also specified.

The single best-fit value from the combined MAE analysis was also ε=6, with best-fit $k_{\mathrm{1}}=\mathrm{5.22}\times {\mathrm{10}}^{-\mathrm{6}}$ min⁻¹ and $k_{\mathrm{2}}=\mathrm{4.12}\times {\mathrm{10}}^{-\mathrm{6}}$ min⁻¹. These best-fit ε, k₁, and k₂ values are consistent when averaging MAE across all years and when only averaging for the 8 years with predicted weakening behavior (defined as years with step drops in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }>\mathrm{2}$ %). Model runs using the combined best-fit parameters only perform better than the constant ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ assumption in 12 out of 23 sample years (mean MAE = 0.217) (Fig. 3a). Both annual and combined model performance are generally better for years which only have ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ strengthening, as evidenced by annual mean MAE = 0.0135 (combined mean MAE = 0.0202) for these years compared to years with observed weakening (annual mean MAE = 0.0163; combined mean MAE = 0.0301). In 11 years, multiple combinations of parameters resulted in the same minimum MAE value. Most of these years only have systematic strengthening, not weakening, so models are relatively insensitive to k₂ and ε values, allowing a range of best-fit model parameterizations. The field data tend to be more variable in years with weakening; nonetheless, the calibrated model captures first-order annual trends across a range of scenarios (Figs. 3e–h, S1 in the Supplement). Across all examples, the dominance of the strengthening effects is demonstrated by B>0.5.

Comparison of mean MAE values of all model runs as a function (k₁, k₂, and ε) elucidates the relative sensitivity of model performance to each parameter (Fig. 3b–d). Model performance was most sensitive to k₁ for the parameter space we explored, with mean MAE values ranging from a minimum value of MAE = 0.0438 at $k_{\mathrm{1}}=\mathrm{8.38}\times {\mathrm{10}}^{-\mathrm{6}}$ min⁻¹ to MAE = 0.181 at $k_{\mathrm{1}}={\mathrm{10}}^{-\mathrm{2}}$ min⁻¹, reflecting the clustering of best-fit k₁ values (Fig. 3b). We note that γ=2.5 was independently calibrated; it is possible that, if our analysis also explored a range of γ values, the range of acceptable k₁ may be broader. In contrast, mean MAE values are higher and less variable when binned by k₂ (MAE = 0.105–0.122) and ε (MAE = 0.109–0.114), suggesting that annual model performance is less sensitive to variations in these parameters (Fig. 3b, c).

6 Discussion

6.1 Model performance

Our flow history model for ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (Eqs. 1, 2) performs better than a constant entrainment threshold, as indicated by lower MAE between data and model (Fig. 3a). The model captures both first-order strengthening and weakening trends seen in the field data. These include progressive increases in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ from lower discharges (Fig. 3e), to sudden decreases in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ following a large flood early in the season (Fig. 3f), to a smaller decrease in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ following a late-season flood after the riverbed may have had more time to strengthen (Fig. 3g), to intermittent but repeated weakening events across a season (Fig. 3h). When each year is calibrated separately (“annual” calibration), the model unsurprisingly performs better than when using the single set of parameters that minimizes MAE averaged across all years (“combined” calibration). Nonetheless, the “combined” best-fit parameters still outperformed a constant ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ assumption in a majority of years – particularly those with seasonal strengthening trends.

Calibrated model performance varies most with k₁, which governs the efficacy of strengthening processes (Fig. 3b, Eq. 1). γ also influences strengthening but was independently calibrated and held constant in our analysis (γ=2.5; Fig. 2b). At the Erlenbach, sediment-transporting flood events only comprise about 2 % of the discharge record (Masteller et al., 2019). Thus, weakening parameters k₂ and ε can only influence ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ evolution during this portion of the record, when ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }>\mathrm{1}$ (Eq. 1). More generally, the transport capacity (${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast })$ at which the model terms combine to transition from overall strengthening to weakening varies for different parameter combinations (Fig. 1). In years without large floods, the strengthening term of Eq. (1) dominates for most of the year, resulting in steady increases in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$, such as in 1988, when floods did not exceed ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }=\mathrm{1.64}$ (Fig. 3e). Therefore, best-fit models for 1988 were insensitive to k₂ and ε, resulting in 248 parameter combinations that minimized MAE (Fig. 3e).

Bedload is transported infrequently in gravel-bed rivers because transport thresholds are rarely exceeded. Much of bedload transport occurs during discharges that are often relatively close to bankfull and when shear stresses only slightly to moderately exceed ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (e.g., Emmett and Wolman, 2001; Parker, 1978; Phillips and Jerolmack, 2016; Pretzlav et al., 2020; Whiting et al., 1999). Therefore, weakening processes, which in our model only occur when ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }\ge \mathrm{1}$, can only reduce ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ for a limited fraction of the full discharge record. Thus, we may expect more generally that strengthening processes are dominant for most of the time in gravel riverbeds relative to weakening processes that may only occur during floods. However, despite strengthening being active for the vast majority of the discharge record, the approximately normal distribution of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ suggests that weakening processes must act more rapidly than strengthening processes to maintain values of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ intermediate between ${\mathit{\tau }}_{\mathrm{cmin}}^{\ast }$ and ${\mathit{\tau }}_{\mathrm{cmax}}^{\ast }$. Future work could further explore the consequences of this difference in the total time over which strengthening and weakening processes may occur and their resultant impact on the time-averaged state of the riverbed, as reflected by distributions of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$. At the Erlenbach, the model is largely successful in matching annual strengthening trends (Fig. 3a, e–h) with a narrow distribution of best-fit k₁ values and a single, independently calibrated value of γ=2.5. This suggests that the physical processes that are encapsulated in the strengthening term of Eq. (1) may lead to predictable changes in similar field settings (e.g., Church et al., 1998; Masteller and Finnegan, 2017; Ockelford and Haynes, 2013).

In contrast, calibrations of weakening parameters find much more variability in best-fit ε and k₂, and MAE is higher for weakening years (Fig. 2a). Annual best-fit results span the full range of both parameters (Fig.  3c, d). High-magnitude floods that cause weakening are relatively rare and short-lived. It is possible that a dataset with many more discharge-driven weakening events could more narrowly constrain these variables. However, a simpler interpretation consistent with our analysis is that weakening is inherently less deterministic than strengthening and therefore more difficult to predict accurately. Weakening events can be de-coupled from flow, for example, if hillslopes supply sediment to the channel during inter-event periods, leading to a reduction in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ at low discharge. Indeed, Masteller et al. (2019) identified a minimum discharge for inter-event strengthening, below which ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ becomes uncorrelated with flow magnitude. This loss of correlation could reflect instances where supply effects introduce event-scale variability in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$. While we focus on evaluating a discharge-driven model for ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$, it is not intended to fully address all factors influencing variability in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$. The discrepancies between observed ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ and model predictions may highlight conditions where sediment supply significantly alters bed mobility, outweighing the flow history effects that are addressed here.

Figure 3e–h illustrates how feedback parameter B controls how ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ changes in response to a given shear stress. Low ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (such as at the starts of 1988 and 2010 and following the large 2007 flood) corresponds to high B, which increases the strengthening term and reduces the weakening term (Eqs. 1, 2). As relatively smaller shear stresses lead to progressive strengthening, ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ increases, B decreases, and less strengthening occurs for a given increment of Shields stress, resulting in a gradual rollover in the rate of strengthening through time (e.g., Fig. 3e). Model response (i.e., $\partial {\mathit{\tau }}_{\mathrm{c}}^{\ast }/\partial t$) is also influenced by ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ through changes in transport capacity. For example, under strengthening conditions, even if τ^∗ remains constant, increasing ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ would cause a gradual decrease in ${\mathit{\tau }}^{\ast }/{\mathit{\tau }}_{\mathrm{c}}^{\ast }$, slowing the rate of strengthening.

Thus, the “memory” in this model is represented by the value of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$, which integrates the effects of the history of both flow conditions and channel bed conditions. Model memories tend to be asymmetric through time in the sense that floods large enough to cause significant weakening will rapidly reset the memory to lower ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ values; strengthening can only occur gradually, as it requires the cumulative effects of lower discharges over time. Conceptually, these memory effects relate to ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ being a state variable for gravel-bed channels (Johnson, 2016). For the Erlenbach, our results using the calibrated model demonstrate that knowing ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ prior to a given flood improves the prediction of transport during that flood. We view this model as a step towards a more complete understanding of mountain river morphodynamics. Our calibrated ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ equation should be useful for improved modeling of channel transport and evolution and as a component of landscape evolution modeling. When high-resolution discharge data are available for field sites, incorporation of a flow-dependent ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ may improve quantitative predictions of transport in gravel-bed rivers, although calibration to local conditions is likely necessary.

6.2 Model implementation beyond the Erlenbach

The calibration that we have performed here leverages an extensive dataset of direct measurements of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$. While similar datasets are available for a small subset of rivers (e.g., Turowski et al., 2011), most gravel-bed rivers lack time series data of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$. However, we find that the range of temporal variability observed in our calibration dataset is consistent with the existing data compilations of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ across a range of slopes (Fig. 2a), suggesting that these data compilations may provide reasonable preliminary constraints of minimum and maximum bounds on ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ in future applications of the model. Reach-averaged starting values of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ could be estimated based on bed grain size and bankfull geometry. Additional model calibration will vary depending on the intended application of the equation. Calibration of the model over approximately 30 transport events may be needed to reliably capture the expected variability in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$. This assumes that ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ are normally distributed, as observed at the Erlenbach by Masteller et al. (2019). However, the commonly used minimum sample size of 30 to characterize normal distributions also assumes independent observations, which does not apply here. An alternative approach could be to calibrate the model based on the number of subsequent events over which ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ remains correlated. Masteller et al. (2019) also found a loss of correlation between ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ values after 10–13 transport events. This number of events may be sufficient to calibrate the model to capture the trajectory of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ over time. We recognize that requiring 10–30 measurements of ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ may not always be feasible. Future studies should assess the necessary level of calibration for different applications. Further research is needed to evaluate the potential variability in the k₁, k₂, and ε parameters especially.

Perhaps the biggest mechanistic limitation of our model is that it only accounts for discharge controls on evolving thresholds, even though sediment supply has also been shown to explicitly influence transport rates in the Erlenbach data (Turowski et al., 2011, Rickenmann, 2020, 2024). In flumes, it is straightforward to impose the upstream sediment supply, measure the flux exiting the flume, and simultaneously measure changes along the flume bed (e.g., surface grain size distributions), allowing thresholds to be evaluated through time as a function of supply (e.g., Johnson, 2016). While it is possible to constrain temporal variations in upstream sediment supply in field settings (e.g., Hassan and Church, 2001; Rickenmann, 2020), these data are far more difficult to measure and less widely available than discharge time series for gravel-bed rivers. We are unaware of field monitoring sites that directly measure comparable time series of transport data in sequential channel reaches, making it difficult to directly isolate supply controls on threshold evolution in gravel-bed rivers. Some sources of sediment supply into a given channel reach, such as shallow landslides a short distance upstream not triggered by recent precipitation, may be uncorrelated with channel discharge. However, the timing and magnitude of many processes that supply sediment to channels, such as bank failures, debris flows, and shallow landslides driven by recent precipitation, are likely correlated with time series of channel discharge (Turowski et al., 2013). In addition, sediment supplied from farther upstream in a watershed is transported into a given reach by channel flow. Seasonal trends in supply (such as from increased hillslope deposition during winter months), followed by subsequent snowmelt or storm flow that progressively transports the sediment (e.g., Moog and Whiting, 1998; Mao et al., 2014), may cause threshold evolution that correlates with cumulative seasonal discharge (e.g., Pretzlav et al., 2020). In other words, discharge and upstream sediment supply are not entirely independent over the timescales of threshold evolution due to floods. Given these correlations between sediment supply delivered to a channel reach from upstream and discharge, the local discharge time series may be able to implicitly account for some temporal variations in local supply and therefore may be able to explain some supply-dependent ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ variability. The degree of correlation between supply time series and discharge time series would likely vary among watersheds based on dominant processes. Future work should attempt to disentangle how sediment supply influences parameter calibrations.

Nonlocal controls on sediment transport will fundamentally limit how well local discharge time series can explain local sediment transport rate. Previous work demonstrates that bedload transport is a nonlocal process because the flux at a location within the channel reflects not only local conditions but also spatial and temporal variations in the flow and sediment flux from upstream (e.g., Foufoula-Georgiou and Stark, 2010; Furbish et al., 2017; Martin et al., 2012). Variability from nonlocality limits the accuracy of all models for calculating bedload flux at a specific location based on local shear stress. This challenge is not unique to our threshold evolution equation. Local flux models also cannot capture spatial and temporal grain dispersion, which is as important as advection for understanding bedload transport through river networks and responses to perturbations (e.g., Bradley, 2017; Pretzlav et al., 2021; Fan et al., 2016). Our model could be applied to better determine the extent to which threshold variability (and associated transport rate variability) is a deterministic function of discharge (as Eq. 1 attempts to represent) and how much of the local transport rate signal is stochastic variability, influenced by a variety of interrelated factors including nonlocality and sediment supply. Future work could also explore how the “memory” of past conditions at a given location, imperfectly encoded in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$, depends on both local discharge variability and nonlocal supply effects.

7 Conclusions

Our study presents a flow-history-dependent model for critical Shields stress ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (Eqs. 1, 2). Calibrated using a 23-year record of river flow and the onset of sediment transport events, our model successfully captures observed trends in ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ evolution at our field site, including seasonal strengthening, rapid weakening following large floods, and gradual riverbed recovery over time. While strengthening processes appear to be relatively predictable across different years, weakening mechanisms exhibit greater variability, likely due to the stochastic nature of high-magnitude floods and external sediment supply effects. Despite these uncertainties, our results demonstrate that accounting for flow history effects by evolving ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ has the potential to significantly improve predictions of sediment mobility and bedload transport compared to a constant transport threshold.

Building on the development and calibration of the model described in this contribution, the best-fit model parameter values found in our study could be used as specific predictions to be independently tested using other field and flume data; we do not yet know how consistent model parameters may be across gravel-bed channels. Model performance could also be assessed at the scale of individual events using continuous bedload measurements, rather than just thresholds at the start of events as done here. The model does not try to isolate granular interaction-based processes that likely cause strengthening and weakening but rather lumps processes together using empirical parameters. Quantifying the systematics, inherent variability, and dominant processes involved in bed weakening warrants additional study. We suggest that a combination of discharge-based controls on ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (as explored here) and sediment supply controls on ${\mathit{\tau }}_{\mathrm{c}}^{\ast }$ (e.g., Recking, 2012; Johnson, 2016; Rickenmann, 2020; Rickenmann, 2024) may be able to explain much of the deterministic variability in threshold evolution and sediment transport rates in gravel-bed rivers.