Gravel threshold of motion: A state function of sediment 1 transport disequilibrium? 2

8 In most sediment transport models, a threshold variable dictates the shear stress at which non-9 negligible bedload transport begins. Previous work has demonstrated that nondimensional 10 transport thresholds ( * c τ ) vary with many factors related not only to grain size and shape, but 11 also with characteristics of the local bed surface and sediment transport rate ( q s ). I propose a 12 new model in which q s -dependent * c τ , notated as * ) ( s q c τ , evolves as a power-law function of 13 net erosion or deposition. In the model, net entrainment is assumed to progressively remove 14 more mobile particles while leaving behind more stable grains, gradually increasing * ) ( s q c τ and 15 reducing transport rates. Net deposition tends to fill in topographic lows, progressively 16 leading to less stable distributions of surface grains, decreasing * ) ( s q c τ and increasing transport 17 rates. Model parameters are calibrated based on laboratory flume experiments that explore 18 transport disequilibrium. The * ) ( s q c τ equation is then incorporated into a simple 19 morphodynamic model. The evolution of * ) ( s q c τ is a negative feedback on morphologic 20 change, while also allowing reaches to equilibrate to sediment supply at different slopes. 21 Finally, * ) ( s q c τ is interpreted to be an important but nonunique state variable for 22 morphodynamics, in a manner consistent with state variables such as temperature in 23 thermodynamics. 24 25


Motivation
Despite over a century of quantitative study (Gilbert, 1914), it often remains challenging to predict gravel transport rates to much better than an order of magnitude because of the complexity of grain interactions with the flow and the surrounding grains (e.g., Schneider et al., 2015;Nitsche et al., 2011;Rickenmann, 2001;Wilcock and Crowe, 2003;Chen and Stone, 2008).Predictive models for complex systems often derive utility from their simplicity, as is the case with the widely-used Meyer-Peter and Müller (1948) transport equation, as modified by Wong and Parker (2006): where q * s is a nondimensional sediment transport rate per unit width, τ * is a nondimensional shear stress imparted by the fluid on the channel bed (a Shields stress), and τ * c is the nondimensional threshold stress at which grains begin to move (a critical Shields stress).Variables are nondimension-alized as follows: where q s is volume sediment transport rate per unit width (m 2 s −1 ), D is grain diameter (m), ρ s is sediment density (m 3 kg −1 ), ρ is water density (m 3 kg −1 ), g is gravitational acceleration (m s −2 ), and τ is shear stress (Pa).In principle, these nondimensionalizations should account for differences in grain size, fluid and sediment density and gravity, allowing meaningful comparisons of transport and stress across different conditions.For a given grain diameter (and constant ρ s , ρ and g assumed for terrestrial landscapes), the simplicity of Eq. ( 1) is that it predicts transport rate using just two variables, τ * (a function of flow strength) and τ * Figure 1.Threshold of motion data from both field and experimental studies.A power law regression to these data gives R 2 = 0.34, indicating that a majority of the variability is not explained by slope alone.Dotted lines indicate common range of τ * c = 0.03 to 0.06 often assumed for modeling transport, although measured data fall well out of this range.Data have been additionally filtered to only include D 50 > 2 mm (i.e., gravel) and slopes between 0.002 and 0.2.Data were compiled and provided by Prancevic and Lamb (2015), based in part on Buffington and Montgomery (1997), with additional data from Olinde (2015) and Lenzi et al. (2006).
it an empirical fitting parameter for a given transport model (e.g., Wong and Parker, 2006;Buffington and Montgomery, 1997).For example, using the original data set of Meyer-Peter and Müller (1948), τ * and q * s give best-fit τ * c = 0.0495 for Eq.(1) (Wong and Parker, 2006).Other bedload transport models have been developed that do not use an absolute threshold stress below which transport is zero, but rather a "reference" stress that corresponds to a very low but nonzero transport rate (e.g., Parker, 1990;Wilcock and Crowe, 2003).For most applications the practical difference between threshold and reference stresses are negligible (Buffington and Montgomery, 1997).In the present work, threshold and reference stresses are used interchangeably.
Thresholds of motion for gravel often span an order of magnitude or more (Fig. 1).Variability in τ * c greatly influences bedload flux predictions in mountain rivers because transport typically occurs close to thresholds conditions, even during large floods (Phillips et al., 2013;Parker et al., 1982;Parker and Klingeman, 1982).Previous work has demonstrated that a great many factors collectively cause τ * Sand content is a related GSD bed state control: increasing sand content of alluvial bed surfaces has been shown to decrease gravel thresholds of motion (e.g., Curran and Wilcock, 2005;Iseya and Ikeda, 1987;Jackson and Beschta, 1984).Wilcock and Crowe (2003) explicitly incorporated this sand dependence into their transport model: where τ * rm is a reference stress (rather than an absolute threshold) for the geometric mean diameter of the bed surface GSD, F s is the spatial fraction of sand on the bed surface, and constants c 1 ,c 2 and c 3 were empirically calibrated from flume data to be 0.021, 0.015 and 20 respectively.These values result in τ * rm varying between 0.021 and 0.036, which is in the range of typical τ * c (Fig. 1).Subsequent work has shown that the effects described by Eq. ( 4) are not unique to sand sizes only.Thresholds of motion for intermediate surface diameters (e.g., D 50 ) can similarly be reduced by grains substantially smaller than the bed surface but larger than 2 mm (Venditti et al., 2010;Sklar et al., 2009;Johnson et al., 2015).Mechanistically, the addition of sand or finer gravels smooths the bed surface by preferentially filling local topographic lows, which can affect pocket geometries (making it easier for larger grains to rotate out of a stable position), and also reduce local hydraulic roughness, increasing near-bed velocity and increasing drag on protruding grains.
Many studies have explored the bed state control of stabilizing structures formed by coarse grain clusters (e.g., Church et al., 1998;Strom and Papanicolaou, 2009).Other bed state controls include the degree of overlap, interlocking and imbrication among grains, and bed compaction or dilation (e.g., Parker, 1990;Wilcock and Crowe, 2003;Sanguinito and Johnson, 2012;Buscombe and Conley, 2012;Mao, 2012;Kirchner et al., 1990;Marquis and Roy, 2012;Powell and Ashworth, 1995;Richards and Clifford, 1991;Ockelford and Haynes, 2013).By combining experimental data and a numerical model, Measures and Tait (2008) show that increasing grain-scale bed roughness tends to shelter downstream grains, reducing entrainment.These factors show that grains can move from less stable to more stable configurations, even if grain size does not change.Coarse grain clusters can also enhance bed stability by increasing surface roughness, tending to deepen potential grain pockets.
Flow characteristics influencing τ * c include particle Reynolds number, flow depth relative to grain size, the intensity of turbulence, the history of prior flow both above and below transport thresholds, and the partitioning of stress into form drag and skin friction (e.g., Shvidchenko and Pender, 2000;Ockelford and Haynes, 2013;Schneider et al., 2015;Valyrakis et al., 2010;Celik et al., 2010).Most flowdependent controls are not independent of the bed surface controls.For example, flow depths, turbulence and form drag depend on slope and bed roughness, while the stress history influences τ * c by changing grain interlocking and surface roughness.Mao (2012) showed that thresholds of mo-tion and bed roughness both evolved during hydrograph rising and falling limbs, leading to bedload hysteresis.
Recent work also suggests that sediment transport can affect τ * c , with higher rates of upstream supply corresponding to more mobile sediment and lower τ * c (Recking, 2012;Bunte et al., 2013).If transport rate influences τ * c then there must be feedback because, by definition, τ * c influences transport rate (Eq.1).Mechanistically, mobile grains impacting stationary grains have been shown to dislodge and entrain grains into the flow (Ancey et al., 2008).Bunte et al. (2013) interpreted that lower τ * c corresponded to looser beds caused by higher rates of sediment supply from upstream, and noted that the stability of bed particles can be qualitatively assessed in the field while doing pebble counts.Yager et al. (2012b) demonstrated that in-channel sediment availability varied inversely with the degree of boulder protrusion, indicating preferential filling of topographic lows by mobile sediment.
Recking (2012) compared bed load monitoring records from steep natural channels (> 5 % slope) to differences in sediment supply interpreted from aerial photographs of surrounding hillslopes.Channels with higher supply rates had higher transport rates for a given shear stress, consistent with a dependence of transport thresholds on supply.While stating that deriving a threshold model "taking into account the sediment input as a parameter would be difficult", Recking (2012) proposed quantitative bounds on reference stress for the end-member cases of very high sediment supply (τ * mss ) and very low sediment supply (τ * m ) in steep mountain channels: (5) It should be noted that these reference stress equations describe reference stresses for the D 84 grain size (rather than D 50 ) (Recking, 2012).The ratio D 84 /D 50 is included in Eqs. ( 5) and ( 6) to represent surface armoring, which tends to vary with sediment supply (Dietrich et al., 1989), thus relating bed state controls to supply-dependent bounds.Overall, my review of previous work suggests that numerous interrelated variables influence τ * c , but also that many controls on τ * c may share similar sensitivities to changing bed roughness and sediment supply.

Morphodynamics and hypotheses
Feedback between channel morphology and bedload transport defines mountain river morphodynamics.The Exner equation of sediment mass conservation quantifies how transport changes correspond to topographic changes (Paola and Voller, 2005): where z is bed elevation (vertical position), x is horizontal position, t is time, and λ p is bed porosity.In this morphodynamic equation (presented for simplicity without an uplift or subsidence term), topographic equilibrium (∂z/∂t = 0) is attained when the sediment flux into a reach equals the sediment flux out (∂q s /∂x = 0).Channel morphology has long been recognized to influence sediment transport.Of particular relevance to the present work, Stark and Stark (2001) proposed a landscape evolution model with a variable called channelization that is defined as representing "the ease with which sediment can flux through a channel reach".(Stark and Stark, 2001).Interestingly, the above definition of channelization could also be applied to τ * c .Because of its control on transport rates, changes in τ * c should influence channel morphodynamics, both over human timescales (e.g., in response to natural and anthropogenic perturbations such as landslides, floods, post-wildfire erosion, land use, changing climate) and longer timescales (landscape evolution).
The overall goal of the present work is to understand and model possible feedbacks among thresholds of motion, changes in transport rate, and the morphological evolution of channels.First, I hypothesize that variability in gravel τ * c is physically meaningful, and that the implicit effects of multiple processes on τ * c can collectively be accounted for in terms of sediment flux dependence.Second, because changes in alluvial channel morphology are strongly coupled with sediment flux (Eq.7), I hypothesize that the evolution of τ * c can implicitly model effects of evolving channel morphology.
The paper is organized as follows.First, I propose a conceptual model for how τ * c should evolve through time as a function of sediment flux (Sect.2.1), and then translate this model into equations (Sect.2.2).Next, I describe flume experiments on disequilibrium gravel transport (Sect.2.3), and use these experiments to empirically calibrate τ * c model parameters (Sect.3.1, 3.2).After that, effects of τ * c evolution on river channel morphodynamics are explored using a simple model for river channel longitudinal profile development (Sect.3.3).Finally, I argue that τ * c is one of many morphodynamic "state variables" that describe how river channels evolve in response to external forcing and internal feedbacks, analogous to state variables in thermodynamics (Sect.4.2).

Conceptual framework for τ * c evolution
The τ * c model proposed below is designed to be applicable at the reach scale, over timescales ranging from changing discharge during floods to the morphodynamic evolution of channels and surrounding landscapes.By definition, models are useful representations of reality because many complexities are omitted.Although recent work demonstrates a richness of threshold and transport behavior caused by turbulent velocity fluctuations and the statistical mechanics of particle populations over short timescales (e.g., Schmeeckle and Nelson, 2003;Diplas et al., 2008;Furbish et al., 2012), these dynamics are not explicitly considered or parameterized in my deterministic model formulation.
Section 1.1 shows that a great many variables and processes influence τ * c .While separate models for every isolated control on τ * c would be informative, it would also be difficult to combine myriad process-specific models and still meaningfully predict the temporal evolution of τ * c for the morphodynamic evolution of channels.Rather than being a process "splitter", I approach the problem as a process "lumper".I hypothesize that many factors affecting grain mobility share common underlying dependencies on net entrainment and net deposition.
Consistent with the form of most bedload transport equations (e.g., Eq. 1), τ * c is defined as a particular Shields stress at which only the most mobile grains of that size become entrained.However, for a population of grains of a given size on the bed surface, there should actually be a distribution of τ * c -notated here as a set of values τ * c -because each individual grain has a particular pocket geometry and near-bed flow velocity at its unique location, and hence a somewhat different individual threshold.Gravel flux increases with discharge primarily because thresholds are gradually exceeded for increasing proportions of surface grains of a given size.For a given transport equation (e.g., Eq. 1), a particular low τ * c value from the distribution τ * c should best predict sediment flux.Conceptually, an underlying assumption is that net entrainment or net deposition changes the underlying τ * c distribution, and therefore changes the value of τ * c that best predicts transport rates.
In the case of a channel reach undergoing net erosion (q sout > q sin ), the most mobile individual grains -i.e., the lowest τ * c values in the τ * c distribution -would preferentially be entrained first, while the grains remaining on the bed would tend to have higher thresholds.Therefore, I hypothesize that progressive erosion tends to entrain grains from increasingly more stable positions on the bed, gradually increasing τ * c .Conversely, during net deposition (q sout < q sin ), I assume that grains tend to preferentially deposit in more stable bed positions such as local topographic lows.Continued deposition tends to occur in progressively less stable positions, gradually decreasing τ * c .These hypothesized τ * changes represent averages for the population of grains; individual grains would exhibit great variability.For example, during net deposition individual grains would also both deposit and be entrained from more and less stable positions, but grains would have a greater probability of remaining deposited in the more stable positions.Mechanistically, τ * c evolution would also be driven by changes in bed topography and roughness, grain clustering and stabilizing structures, compaction of the bed and interlocking of grains, etc.None of these physical variables are explicitly included in the model equations; instead their combined effects are assumed to vary with net erosion or deposition.Importantly, the amount by which τ * c changes should also depend on the current state of the bed surface.Starting from a relatively rough and interlocked bed surface, net deposition would initially cause relatively large decreases in bed roughness as local lows preferentially filled with loose grains, and relatively large corresponding decreases in τ * c .However, for a given surface GSD there must be physical limits for bed roughness and grain packing.If bed surface grains are already relatively loose and mobile, additional deposition would cause less decrease in τ * c , or no decrease if the bed is already as unstable as possible for a given surface GSD.Thus, the change in τ * c should also be a function of τ * c .The processes that cause changes in τ * c also place physical limits on how high, and low, τ * c can evolve.These τ * c dependencies describe negative transport feedbacks: net erosion progressively reduces rates of erosion by making grains harder to entrain, while net deposition progressively makes grains more mobile.It has long been recognized that channel reaches evolve towards steady-state configurations in which the sediment flux into a reach balances the flux exiting, leading to zero net erosion or deposition (Mackin, 1948).At this statistical steady state, τ * c should also be at equilibrium.If τ * c were still systematically evolving (e.g., from continued bed state changes), then transport rate through the reach would also change, perturbing the channel away from its statistical equilibrium.

model equations
While the above discussion makes the case that τ * c can evolve through time due to a variety of interrelated factors, the new model proposed here is specifically in terms of sediment flux.I use the notation τ * c(q s ) to distinguish this specific model from more general thresholds of motion (i.e., τ * c ) in other models and analyses.Because longitudinal coordinate x increases downstream, net erosion in a reach is indicated by ∂q s /∂x > 0 and net deposition by ∂q s /∂x < 0 (Eq.7).The following relations are proposed: where κ ent and κ dep are dimensionless exponents corresponding to entrainment and deposition, respectively, and k is a scaling factor.These three parameters will be empirically fit to experiments.τ * cmin and τ * cmax represent bounds on how low or high τ * c(q s ) can plausibly evolve (assumed to be 0.02 and 0.35 respectively).Equation ( 8) predicts that τ * c(q s ) incrementally decreases with net deposition, and incrementally increases during net erosion."Feedback factor" B has a value between 0 and 1 and makes Eq. ( 8) a differential equation.It scales the incremental change in τ * c(q s ) so that deposition on an already "loose" bed (τ * c(q s ) close to τ * cmin ) minimally decreases τ * c(q s ) , but erosion causes a larger τ * c(q s ) increase.Conversely, if τ * c(q s ) is already high (close to τ * cmax ), then erosion causes a much smaller τ * c(q s ) change than does deposition.Finally, I note that representing ∂τ * c(q s ) /∂t as a function of ∂q s /∂x (Eq.8) is broadly analogous in form to Exner (Eq.7).
A limitation of Eq. ( 8) is that, for dimensional consistency, the units of k vary with κ ent and κ dep .An improved equation replaces spatial changes in flux with spatial changes in the thickness of deposited or eroded sediment: θ s is the thickness of sediment deposited or eroded at a given location.∂θ s /∂x is a dimensionless ratio representing spatial changes in erosion and deposition.In this case, k has dimensions 1/t and scales how quickly τ * c(q s ) evolves.θ s can be calculated by integrating Eq. ( 7) over time interval t 1 to t 2 : for a generic function f ).Using discrete flume data, θ s is calculated over a measurement interval t as 1 − λ p −1 (q sout − q sin ) t/ x, where x is the length of the flume and the sediment flux terms are averaged over t.
A r is a dimensionless armoring parameter, calculated below in several ways in order to explore whether predictions can be improved by explicitly including bed surface grain size or bed roughness characteristics.Setting A r = D 50 /D 84 (the reciprocal of the Recking (2012) Figure 2. Flume experiment data (Johnson et al., 2015).(a) Sediment transport rate in (Q sfeed ) and out of the flume.The upstream sediment supply rate was zero other than during the Q sfeed period.Experiment 1 was run for a longer duration than the others but shows similar trends.Note that the outlet Q s adjusts much faster to match the increase in supply than it does to decrease during periods of no input.(b) Median bed surface grain diameters decreased during the feed of finer gravel, and then increase beyond their previous stable bed.(c) Flume-averaged bed slopes changed relatively little even as transport rates and D 50 changed greatly in response to initial bed stabilizing and supply perturbations.
are larger where D 50 is relatively closer to D 84 .Setting A r = 2D 50 / (D 84 − D 16 ) suggests that τ * c(q s ) changes should be larger when intermediate diameters are large relative to a measure of the normalized width of the bed surface GSD.I also try A r = D 50 /σ , where σ is bed surface roughness.A r = D 50 /σ suggests that, relative to topographic lows and highs, large grains cause bigger τ * c(q s ) changes than small grains.Finally, A r is simply set to 1 in many calculations below.

Experimental design
The flume experiments used to calibrate k, κ ent and κ dep were designed to explore feedback during disequilibrium transport in gravel-bed rivers.Figure 2 shows how transport rates, surface D 50 and bed slope evolved in response to fine gravel pulses.Johnson et al. (2015) provide details of the experimental conditions and how they scale to natural conditions most consistent with step-pool development, and so the summary here is brief.Four experiments were conducted in a small flume 4 m long and 10 cm wide.Experiments 1 and 4 were done at 8 % initial slope, and 2 and 3 at 12 % initial slope; slopes subsequently evolved fairly little during morphological adjustment (Fig. 2c).Water discharge was held constant throughout to better isolate the influence of sediment supply changes on transport.Sediment transported out of the flume was caught in a downstream basket, sieved and weighed.Overall sediment diameters ranged from 0.45 to 40 mm; these sizes were sorted and painted different colors based on five size classes with D 50 = 2.4, 4.5, 8.0, 15.4, and 27.2 mm (D 16 = 2.0, 3.4, 6.7, 12.4, and 24.0 mm; D 84 = 2. 8, 5.7, 10.3, 19.7, and 31.3 mm, respectively).Surface GSDs were measured using image analysis of colored bed surface grains during the experiments.Bed topography was measured using a triangulating laser, and bed roughness (σ ) was calculated from longitudinal topographic swaths as the standard deviation of detrended bed elevations.Water surface elevations were measured using an ultrasonic distance sensor, and water depths were calculated by subtracting bed elevations.Total shear stress (τ ) was calculated assuming steady uniform flow when spatially averaged over the flume: where h is water depth corrected for sidewall effects following the method of Wong and Parker (2006), and S is water surface slope.
The experiments started with mixed-size sediment screeded flat.Initially, all surface sizes were observed to be mobile (and therefore above thresholds of motion).At the beginning no sediment was fed into the upstream end (q sfeed = 0), and the bed responded by coarsening, roughening and gradually stabilizing as transport rates dropped by ≈ 3 orders of magnitude (Fig. 2a).After this initial stabilization, a stepfunction pulse of the finest gravel size (D 50 = 2.4 mm) was fed into the flume at q sfeed = 1000 g min −1 , representing an idealization of a landslide, debris flow, post-wildfire erosion, or anthropogenic gravel augmentation that would suddenly supply sediment finer than the existing bed surface.The feed rate was chosen to be similar to the high initial transport rates (Fig. 2a), while not so high as to inhibit morphodynamic feedbacks by fully burying the existing bed surface.Initially some deposition occurred, but the channel adjusted rapidly, by both entraining coarser bed surface grains and transporting most of the finer supplied gravel, so that the outlet transport rate (q sout ) approximately matched q sin .After that the sediment supply pulse q sin was again dropped to zero, and the bed gradually restabilized.Johnson et al. (2015) explained in detail how bed roughness evolved, and how the addition of finer gravels ultimately caused surface coarsening (Fig. 2b).Unbalanced transport rates into and out of the flume demonstrate disequilibrium conditions (Fig. 2a), although transport evolved towards equilibrium.

The Wilcock and Crowe (2003) transport model
To quantify thresholds of motion from these experimental data (Fig. 2) requires a transport model.The Wilcock and Crowe (2003) "Surface-based Transport Model for Mixed-Size Sediment", abbreviated as W&CM, is used for two main reasons.First, the model can, at least in principle, account for the effects of changing surface GSD on τ * c .Second, the model should also be able to account for possible effects of sand and fine gravel abundance on thresholds of motion (Eq.4).By using the W&CM to isolate and remove GSD effects, experimentally constrained thresholds of motion can then be used to evaluate the proposed τ * c(q s ) functions (Eqs.8-11).A secondary goal is to evaluate how well the W&CM predicts disequilibrium transport at steeper slopes and lower water depths than Wilcock and Crowe (2003) used in their own steady-state experiments.
A key variable in the W&CM is τ * rs50 , the nondimensional reference stress for the median surface grain size (D s50 ).τ * rs50 corresponds to a very low transport rate of W * i = 0.002.W * i is a nondimensional bedload transport rate for grain size class i, where q bi is the volumetric transport rate per unit channel width of grains of size i, F i is the fraction of size i on the bed surface, and u τ is shear velocity (u τ = √ τ/ρ).Wilcock and Crowe (2003) presented an empirical relationship between transport and shear stress: where τ ri is a dimensional reference stress for size class i, with dimensionless equivalent τ * ri (Eq.3).A "hiding function" determines how nondimensional reference stresses vary with grain size: The hiding function exponent b is calculated as b = 0.67 1 + e (1.5−D i /D s50 ) . (16) Note that Eq. ( 16) is slightly modified from the exact Wilcock and Crowe (2003) version by replacing D sm (the geometric mean surface diameter) with D s50 , to more simply use just one measure of the central tendency of the surface GSD.16) to calculate hiding function exponent b, while "Power-law fit" calculates a best-fit b along with τ * rs50 .Error bars give 95 % confidence intervals on τ * rs50 based on the regressions; although uncertainty can be broad the trends are clear and consistent.Shaded area indicates times of fine gravel addition (sediment feed) in each experiment.Data in these plots are available in Johnson (2016).

Results
In this section, the W&CM is used to calculate best-fit thresholds of motion.Next, the τ * c(q s ) model is shown to predict the experimentally constrained τ * rs50 trends after calibrating several parameters.Finally, the influence of τ * c(q s ) on morphodynamics is explored using a simple model for gravel-bed river profile evolution.

Best-fit τ * rs50 and hiding functions
The experimental data are used to determine W * i (Eq.13) and τ * ri for each of the five grain size classes (Eq.14).Bestfit τ * rs50 is then calculated in two ways.In the first approach, b is calculated using Eq. ( 16), and τ * rs50 and 95 % confidence intervals are estimated using nonlinear multiple regression in Matlab (Fig. 3, "W&CM fit").Importantly, the temporal evolution of best-fit τ * rs50 is not explained by grain size changes, because the W&CM already accounts for the effects of surface GSD (Fig. 3).
While b varies with relative grain size in the W&CM (Eq.16), other proposed hiding functions have found (or assumed) that a single b value applies to different grain sizes, at least for a given set of flow and surface conditions (Parker, 1990;Buscombe and Conley, 2012).My second approach for estimating τ * rs50 explores whether the results are sensitive to  The W&CM hiding function (Eq.16) does a good job matching the data, although it was not fit to these points.The first six measurements of each experiment (roughly the first 10 min) were excluded because of large scatter associated with the greatest bed instability.The axes reflect the left and right hand sides of Eq. ( 15), but uses dimensional stresses to be consistent with plots shown in Wilcock and Crowe (2003).
the particular form of Eq. ( 16).Rather than Eq. ( 16), nonlinear multiple regression was used to estimate both b and τ * rs50 in Eq. ( 15), with separate regressions for each time step (Fig. 3).The temporal evolution of experimental τ * rs50 is generally comparable for the two different approaches (Fig. 3).
Interestingly, Fig. 4 shows that the hiding function exponents determined using the nonlinear multiple regressions for b and τ * rs50 are consistent with Eq. ( 16) of Wilcock and Crowe (2003).In spite of substantial scatter there is a slope break which corresponds to a change in b for surface grains smaller and larger than the median, suggesting that the W&CM reasonably can describe hiding and exposure relations among grains in steeper channels and for shallower flow depths than used in the Wilcock and Crowe (2003) experiments.
3.2 Calibration of τ * c(q s ) model parameters Figure 5 compares experimentally constrained thresholds of motion to several predictions of these trends.First, I test whether surface sand fraction (Eq.4) can explain the evolution of τ * rs50 (Curran and Wilcock, 2005;Wilcock and Crowe, 2003).As described in Sect.1.1, the effect of finer grains on thresholds of motion is not limited to sand sizes alone (Venditti et al., 2010;Sklar et al., 2009;Johnson et al., 2015).In the Johnson et al. (2015) experiments, the finest grain size class has D 50 = 2.4 mm, D 16 = 2.0 mm, D 84 = 2.8 mm.Setting F s equal to the surface fraction of this size class, a nonlinear multiple regression of Eq. ( 4) to all four experiments together yielded a poor although statistically significant fit to the data (R 2 = 0.13; p = 3 × 10 −5 ; c 1 = 0.097 ± 0.04, c 2 = 0.103 ± 0.11, and c 3 = 5.6 ± 11), confirming that surface grain size changes alone cannot explain observed τ * rs50 patterns (Fig. 5, "Sand fraction").Note that I have assumed for simplicity that τ * rs50 = τ * rm , i.e., substituting the surface D 50 for the geometric mean surface diameter in Eq. ( 4).
Various τ * c(q s ) models provide better fits to experimentally constrained τ * rs50 (Fig. 5; Eqs. 8, 10).Models with A r = 1 are shown using a single set of model parameters for all four experiments ("collective best fit", k = 0.17, κ dep = 0.20, κ ent = 0.40), and also the best fit for each experiment separately.The best-fit overall model has R 2 = 0.69, suggesting statistically that effects of supply and transport disequilibrium can explain over 2/3 of the variability in τ * rs50 (Table 1).Note that τ * c(q s ) and τ * rs50 are assumed to be interchangeable.Because Eqs. ( 8) and (10) are differential equations, bestfit parameters could not be calculated using nonlinear multiple regressions.Instead, I use a brute-force approach of incrementally stepping through a wide range of k, κ dep and κ ent , and finding the combination of parameters that give the smallest root-mean-square deviation (RMSD).These calculations started at τ * rs50 = 0.036 at t = 0, which is consis- tent with the experiments, and is also the τ * rs50 proposed by Wilcock and Crowe (2003) in the absence of sand dependence.
Model fits using A r = D 50 /σ are not substantially different from A r = 1, and R 2 = 0.69 is the same (Fig. 5).Table 1 includes additional regressions for A r = D 50 /D 84 and A r = 2D 50 / (D 84 − D 16 ).These fits overlap almost perfectly with those shown on Fig. 5.As explained in Sect.2.2, A r should account for surface GSD and bed topography influences on thresholds.The fact that regressions are not improved by including these variables may suggest that transport disequilibrium is a more important control on threshold evolution over a broad range of surface GSD and bed roughness.Similarly, Eq. ( 8) (dimensional ∂q s /∂x) does not fit the data better than Eq.(10) (Table 1).Because these variants do not substantially improve τ * c(q s ) model fits, I use the simplest dimensionally consistent model (Eq. 10 with A r = 1) in the analysis below.

on morphodynamics
Next, an idealized morphodynamic model demonstrates how the proposed τ * c(q s ) relations influence the evolution of channel profiles, focusing on reach slopes and timescales of adjustment.Because the modeling goal is to isolate and understand effects of evolving τ * c(q s ) , the underlying model is arguably the simplest reasonable representation of morphodynamic feedback.Inspired by Parker (2005), the model describes a channel reach in which slope evolves through aggradation and degradation.The downstream boundary elevation is fixed (constant base level).Sediment transport and bed elevation are modeled using Eq.(1) (substituting τ * c(q s ) for τ * c as needed) and Eq. ( 7) with a single grain diameter (D).Unit water discharge q w is held constant for simplicity.Upstream sediment supply rate (q sfeed ) is imposed, and is varied to drive channels to new steady states.Relationships among flow depth, depth-averaged velocity and discharge are imposed by assuming that hydraulic roughness remains constant, parameterized though a Darcy-Weisbach hydraulic friction coefficient: For a given discharge this allows both U and h to be determined: Two model variations are compared: in the "Exneronly" morphodynamic model, τ * c is a constant.In the "Exner +τ * c(q s ) " variant, τ * c(q s ) evolves through time following Eq.( 10).At equilibrium, channel slope can be predicted for both model variants (and substituting τ * c(q s ) for τ * c where appropriate) by combining Eqs.(1), (2), ( 12) and ( 19): For a given discharge, Eq. ( 20) indicates that both sediment supply and the threshold of motion influence steady-state morphology (slope).
Away from equilibrium, rates of change of bed elevation along a river profile should depend not only on the sediment flux at a given channel cross section, but also on the average velocity at which grains move downstream.This control has occasionally been ignored in previous models of profile evolution.In my model, it is crudely incorporated by assuming that average bedload velocity is a consistent fraction of water velocity, broadly consistent with previous findings that bedload velocities are proportional to shear velocity (e.g., Martin et al., 2012).The modeling timestep is set to be equal to the time it takes sediment to move from one model node (bed location) to the next, and is adjusted during simulations.While this approach makes the temporal evolution of channel changes internally consistent within the model, timescales for model response will still be much shorter than actual adjustment times in field settings because flood intermittency is not included (the model as implemented is always at a constant flood discharge).In addition, the upstream sediment supply is imposed in the model, while in natural settings hillslope-floodplain-channel coupling could greatly affect q sfeed over time if significant aggradation or downcutting took place. .Profile evolution, comparing the morphodynamic responses of models with and without threshold evolution.The initial condition is an equilibrium channel with τ * c(q s ) = 0.06, upstream sediment supply q s = 1 × 10 −3 m 2 s −1 , and an initial equilibrium slope of 0.0147.Sediment supply is increased 5× at t = 0. Lines are each 5 model days apart, and indicate the evolution to a new transport equilibrium.
Table 2 provides parameters used for morphodynamic modeling.Although the highly simplified model is not intended for quantitative field comparisons, variables D(D 50 = 50 mm), f (0.1), and q w (1 m 2 s −1 ) were chosen to be broadly consistent with a moderate (≈ 2-3 year peak discharge recurrence interval) bedload-transporting flood in Reynolds Creek, Idaho (Olinde and Johnson, 2015).Reynolds creek is a snowmelt-dominated channel with reach slopes that vary widely from ∼ 0.005 to 0.07.In an instrumented reach with a slope of 0.02, Olinde (2015) used RFID-tagged tracers and channel-spanning RFID antennas to measure τ * rs50 ≈ 0.06.A constant τ * c = 0.06 is used for the Exner-only models, while τ * c(q s ) = 0.06 is used as the initial condition for Exner +τ * c(q s ) models.Field constraints on upstream sediment feed rates were not available, and so q sfeed values were chosen to provide reasonable model slopes.Exponents κ dep and κ ent used the experimental calibrations, while k were chosen so that changes in τ * c occurred over the same range of timescales as topographic adjustments, to better illustrate the interplay of variables in morphodynamic evolution (Table 2).

Morphodynamic model results
Figures 6 and 7 compare how longitudinal profiles respond to an increase in sediment supply, for both the Exner-only (constant τ * c ) and Exner +τ * c(q s ) models.The initial condition is a channel at equilibrium (q sout = q sfeed ).At t = 0, sediment supply is increased by a factor 5. The Exner +τ * c(q s ) model aggrades to a new equilibrium slope that is lower than the Exner-only model.This occurs because deposition (∂q s /∂x < 0) causes τ * c(q s ) to decrease over time, progressively increasing transport efficiency (i.e., higher transport rates at a lower slope) compared to constant τ * c = 0.06 (Fig. 7).Feedback causes the reverse effect for a decrease in q sfeed : τ * c(q s ) progressively increases as the bed erodes  6, t = 0 corresponds to an equilibrium condition where the initial slope and initial threshold are consistent with the initial upstream sediment supply.Slope and τ * c(q s ) were averaged over nodes 3-10, leaving out the first and last two nodes because of minor model boundary effects.and slope decreases, leading to channel re-equilibration both sooner and at a higher slope (Fig. 7).
I measure an equilibrium timescale (t eq ) as the amount of time it takes from a supply perturbation (t = 0 in these models) to the slope adjusting to be within 0.0001 of its equilibrium slope (Eq.20).In Fig. 7, t eq are substantially longer for the Exner-only models than for the Exner +τ * c(q s ) models.For Exner +τ * c(q s ) , an increase in q sfeed leads to aggradation, in turn increasing local q * s by both increasing slope and also decreasing τ * c(q s ) (Eqs. 1, 10).Both factors adjusting enable equilibrium to be reached sooner.
Over a q sfeed range of 2 orders of magnitude, equilibrium slopes change less for the Exner +τ * c(q s ) model than for Exner-only (Fig. 8a).The ratio of these equilibrium slopes illustrates the magnitude of the change, where "S eq ratio" is S eq for Exner +τ * c(q s ) divided by Exner-only S eq (Fig. 8b).An order-of-magnitude decrease in q sfeed caused Exner +τ * c(q s ) S eq to be roughly 24-36 % larger than Exner-only S eq .An order-of-magnitude increase in q sfeed caused Exner +τ * c(q s ) S eq to be roughly 20 % smaller than the constant-τ * c model S eq .Calculations are also shown for several values of scaling factor k. A larger k means that τ * c(q s ) increases or decreases more rapidly for a given amount of aggradation or degradation (Eq.10), which in general enables a new equilibrium to be reached with a smaller change in slope.
Equilibrium timescales are quite sensitive to k as well as to sediment supply rate (Fig. 8c).Similar to the S eq ratio, the "t eq ratio" is t eq for Exner +τ * c(q s ) , divided by t eq for the Exner-only model (Fig. 8d).There is an asymmetry in equilibrium times for aggradation vs. degradation; in general the difference between Exner-only and Exner +τ * c(q s ) is  20) and the morphodynamic model calculations demonstrate that the models did asymptotically attain equilibrium slopes.(b) S eq ratio is the ratio of equilibrium slopes of the Exner +τ * c(q s ) model divided by S eq for the Exner-only model, to show the relative affect that τ * c evolution has on equilibrium slopes.(c) Equilibrium timescales for model adjustment.(d) t eq ratio is the ratio of t eq for the Exner +τ * c(q s ) model divided by t eq for the Exner-only model.Values are lower than 1, indicating that the τ * c evolution has a large influence on equilibrium timescales.(e) Evolution of τ * c(q s ) .
somewhat smaller during bed aggradation, and the difference decreases with increasing q sfeed .Interestingly, the highest k (2.8 ×10 −5 ) results in a threshold-like response where the t eq ratio rapidly increases from roughly 0.01 to 0.8 (Fig. 8d).This change occurred because τ * c(q s ) "bottomed out", i.e., reached its minimum possible value (τ * c(q s ) ≈ τ * cmin = 0.02) before the equilibrium slope had been reached (Fig. 8e).At that point, τ * c(q s ) could no longer act as a buffer to reduce slope changes, and it took much longer to reach an equilibrium slope.
Finally, Fig. 9 shows that the spatial as well as temporal evolution of τ * c(q s ) can influence river profiles.The models are the same as in Fig. 6.At t = 0, the feed rate into the upstream-most node (q sfeed ) increases by a factor of 5. Therefore, the upstream end feels the supply perturbation both sooner and more strongly than downstream nodes.Aggradation from the supply perturbation increases upstream slopes first.In the Exner-only model, downstream slopes gradually catch up.Because τ * c stays constant, every location along the channel eventually asymptotes to the single slope required to persist even at equilibrium.
transport the new q sfeed at the given discharge (Fig. 9a).However, for evolving thresholds, enhanced upstream aggradation caused upstream τ * c(q s ) to decrease both more rapidly and to lower values than downstream nodes.Spatial differences in τ * c(q s ) persisted at equilibrium, resulting in spatial variations in equilibrium slope (Exner +τ * c(q s ) ; Fig. 9b, c).

Discussion
In this section, the dependence of τ * c(q s ) on sediment supply is compared to previous work.τ * c(q s ) evolution is identified as a negative feedback on morphologic change that can impart a memory of previous channel "states" to the system.Finally, τ * c(q s ) is interpreted as a channel state variable, analogous to temperature in thermodynamics.
As described in Sect.2.1, previous work on sediment supply-dependent thresholds of motion includes Recking (2012), who proposed high sediment supply (τ * mss ; Eq. 5) and low sediment supply (τ * m ; Eq. 6) bounds on thresholds of motion.Figure 10 shows how these relations compare to the experimentally constrained τ * rs50 .It should be noted again that these bounds were calibrated to the D 84 grain size rather than D 50 (Recking, 2012).While the actual values are therefore not expected to be equivalent, τ * mss and τ * m do tend to bound τ * rs50 .The low-supply bound τ * m is roughly 2-4 times larger than the experimental constraints.The high-supply bound τ * mss is similar in magnitude to τ * rs50 and predicts the decrease during the feed period.The (linear) correlation between τ * mss and τ * rs50 is weak (R 2 = 0.13) although statis- tically significant (p = 3 × 10 −5 ).Nonetheless, given that threshold of motion uncertainties are typically large, Eq. ( 5) arguably provides a reasonable independent prediction of our experimental disequilibrium transport data, based on experimental slope, D 84 and D 50 .
The τ * c(q s ) model is consistent with previous interpretations that high sediment supply corresponds to low thresholds of motion, and vice-versa (Recking, 2012;Bunte et al., 2013).In the τ * c(q s ) model (Eq.10), an increase in upstream sediment supply that causes net aggradation will lower τ * c(q s ) , unless τ * c(q s ) has already reached its lower physical limit (τ * cmin ).Conversely, a decrease in supply that causes net erosion will increase τ * c(q s ) , unless τ * c(q s ) is already high (≈ τ * cmax ).However, while the τ * c(q s ) model can thus explain an inverse relation between supply and thresholds of motion, it is worth noting that Eqs. ( 8) and (10) describe a subtly different feedback: τ * c(q s ) does not directly increase or decrease with supply, but rather with the history of sediment supply changes relative to transport capacity over time.If q sin equals q sout , τ * c(q s ) could remain constant regardless of whether q sin is high or low.

Negative feedback and asymmetric approaches to equilibrium
The evolution of τ * c(q s ) acts as a negative feedback because it reduces the morphodynamic response to perturbations.Reach slopes and τ * c(q s ) both change in the direction that brings transport back towards equilibrium, allowing smaller slope changes to accommodate supply changes (Figs. 6,7,8a,b,9).However, as with other buffered systems, there is a limit to how large of a perturbation can be accommodated by τ * c(q s ) (as illustrated by k = 2.8 × 10 −5 in Fig. 8c, d, e).The amount of possible τ * c(q s ) change depends on how close τ * c(q s ) is to τ * cmin or τ * cmax (Eq.9).When changes in τ * c(q s ) are negligible but transport and morphology are not equilibrated, then the time to equilibrium (t eq ) increases because only channel morphology can adjust (Fig. 8c, d, e).
The experiments suggest that τ * c(q s ) changes faster in response to aggradation than degradation (Figs. 2, 5).This asymmetry is expressed in the best-fit exponents: κ dep is smaller than κ ent for all scenarios tested (Table 1).Note that because ∂θ s /∂x is much smaller than 1 (i.e, spatial changes in bed elevation are small compared to the horizontal distance the change is measured over), the smaller exponent (κ dep ) corresponds to a larger change in τ * c(q s ) for a given ∂θ s /∂x (Eq.10).For a given increment of sediment thickness (θ s ), aggradation is more efficient at decreasing τ * c(q s ) than degradation is at increasing τ * c(q s ) .Future work is required to explore how specific physical processes vary during net deposition or erosion and lead to asymmetry in τ * c(q s ) change.Still, a tentative hypothesis linking bed roughness and τ * c(q s ) change asymmetry is that during deposition, clasts tend to deposit preferentially in topographic lows, because these tend to be the most sheltered locations, and simply because of the direction of gravity.Preferentially filling in lows tends to directly decrease bed roughness, in turn reducing topographic sheltering and hydraulic friction and increasing near-bed flow velocities.All of these factors decrease τ * c(q s ) .However, erosion does not simply have an opposite but symmetric effect on bed topography as deposition.Clasts are not preferentially eroded from topographic lows, as these locations tend to remain the most sheltered.Instead, the process of increasing bed roughness during erosion is more complex and results from the more gradual development of stabilizing structures around keystones, as grains are rearranged locally to positions where they protrude into the flow but remain stable due to interlocking with surrounding grains.Thus, roughness reduction and enhancement should not be equally sensitive to net erosion or deposition.Mao (2012) showed that bed roughness evolved at different rates during symmetric rising and falling limbs of hydrographs, influencing gravel transport hysteresis.Bed roughness due to sand ripple and dune evolution has also been shown to increase and decrease at different rates during hydrograph rising and falling limbs, leading to hysteresis in a sand transport system that is not threshold dominated (Martin and Jerolmack, 2013).
In Fig. 8c, the Exner +τ * c(q s ) model indicates that equilibrium timescales are longer for aggradation (q sfeed / initial q sfeed > 1) than for degradation.At first glance this seems to contradict the argument that aggradation is more efficient at decreasing τ * c(q s ) .The explanation is that the equilibrium timescale does not only depend on the exponents, but also on how much total aggradation or degradation occurs to attain equilibrium.Slope changed more during aggradation than degradation for these particular Exner +τ * c(q s ) models, even though τ * c(q s ) also tended to change more during aggradation than degradation (Fig. 8a, e).
In the experiments, average slopes changed very little in response to changes in sediment supply and transport disequilibrium, while grain size and bed roughness changed much more (Fig. 2; bed roughness is presented in detail in Johnson et al., 2015).Because the W&CM already accounts for grain size changes in determining experimental τ * rs50 (Fig. 3), bed roughness and various unquantified mechanisms (such as grain interlocking) are interpreted to have physically caused the τ * rs50 evolution.What does this suggest for k, which scales how much τ * c(q s ) changes for a given amount of aggradation or degradation?The best-fit k was 2.83 × 10 −3 s −1 , which reflects the rapid adjustment of experimental τ * c(q s ) compared to slope changes (Fig. 5, Table 1, Eq. 10).In contrast, in the morphodynamic modeling I used k values 2 to 3 orders of magnitude smaller, so that the response to a perturbation in supply would involve nonnegligible changes in slope (the only morphologic variable in the simple morphodynamic model) as well as changes in τ * c(q s ) .Higher values of k in the morphodynamic model cause τ * c(q s ) to adjust more rapidly and slope to adjust less (Fig. 8).
An implication of τ * c(q s ) evolving with reach morphodynamics is that local channel form can retain "memories" of previous conditions, which can influence local responses to subsequent forcing.In Fig. 9b and c, an increase in supply led to the temporal and spatial evolution of τ * c(q s ) , which in turn caused spatial variations in equilibrium slope.Upstream reaches acted as filters of the supply perturbation to downstream reaches.In nature, spatially and temporally averaged morphodynamic equilibrium will reflect "channelforming" discharges and a representative sediment supply from upstream, but floods, local supply perturbations and history add to spatial variability in both τ * c(q s ) and morphology.I also acknowledge that the model parameters were calibrated to flume experiments at steep 8 and 12 % slopes with a GSD that includes scaled boulders (Table 1; Johnson et al., 2015); future work is required to determine how the surface GSD influences the strength of τ * c(q s ) evolution, and how well the model predicts τ * c(q s ) changes in lower slope gravel-bed rivers.

State variable framework for modeling morphodynamics
Next, I argue that τ * c(q s ) should be redefined as a state variable (or state function) for gravel-bed channels, and outline a possible state variable approach for modeling the morphodynamic evolution of channels.The term "bed state" has long been informally used to describe collective aspects of local channel morphology, such as surface GSD and armoring and clustering, that change with relative ease and influence transport rates (e.g., Church, 2006;Gomez and Church, 1989).Although explicitly defining τ * c(q s ) evolution and related feedbacks in terms of state and path variables appears to be novel (to my knowledge), channel morphodynamics have long been implicitly described using similar ideas.For example, Phillips (2007) presented a qualitative conceptual model of landscape evolution in terms of improbable system states, arguing that although deterministic process "laws" act on topography, the actual outcome (i.e., any particular landscape) depends on initial conditions and in particular is sensitive to history.Many other works have similarly generalized complex channel process and response feedbacks to understand morphodynamics (e.g., Fonstad, 2003;Phillips, 2011Phillips, , 2009Phillips, , 1991;;Chin and Phillips, 2007;Stark and Stark, 2001;Yanites and Tucker, 2010).
State variables are integral to many disciplines, including control systems engineering and thermodynamics.Thermodynamic state variables include temperature, pressure, enthalpy and entropy.By definition state variables are pathindependent (Oxtoby et al., 2015).For example, temperature (T ) describes the amount of thermal energy per unit of a material.A change in temperature depends only on the initial and final states (i.e., T = T 2 − T 1 ), but does not depend on the path, i.e., the history of temperatures between times t 2 and t 1 .In contrast, heat -the flow (transfer) of thermal energy -is a path variable (or process variable), not a state variable.Heat flow between bodies is both controlled by and changes the temperature of those bodies, but the amount of total heat transferred does depend on the path.Three other points about state variables are relevant to morphodynamics.First, state variables are rarely independent of one another.For example, Gibbs free energy is a state variable calculated from temperature, enthalpy and entropy (Hemond and Fechner, 2014).Second, although state variables are technically only defined at equilibrium, in practice they are useful for understanding gradually evolving systems (e.g., Kleidon, 2010).Third, the evolution of systems involving multiple state variables are usually described with coupled differential equations.
Channel morphodynamics can be described by a similar framework of state and path variables.Analogous to heat, the cumulative discharges of both water and sediment are path variables that drive bed state evolution.Channel morphology can be described by numerous bed state variables, including but not limited to surface GSD, slope, width, depth, bed roughness, surface grain clustering, interlocking, overlap, imbrication, and finally τ * c(q s ) .Analogous to temperature, I explicitly define τ * c(q s ) as a state variable.The amount of change in τ * c(q s ) from time t 1 to t 2 does not depend on the progression of values in between.However, the amount of sediment transported between t 1 and t 2 does depend on the history of τ * c(q s ) , and also influences the history of τ * c(q s ) (Eqs. 8, 10).
Entropy is the state variable perhaps used most often to characterize channel systems (e.g., Chin and Phillips, 2007;Leopold and Langbein, 1962;Rodriguez-Iturbe and Rinaldo, 1997).Entropy can provide a closure for underconstrained sets of equations, by assuming that geomorphic systems inherently maximize their entropy at equilibrium (Kleidon, 2010;Chiu, 1987).A limitation of some maximum-entropy landscape models is that physically based surface processes are not always explicitly modeled, making them less useful for predicting landscape responses to environmental perturbations, even if they can create reasonable equilibrium morphologies (Paik and Kumar, 2010).In contrast to entropy, state variable τ * c(q s ) has a clear process-based meaning.I suggest that landscape evolution models could incorporate subgrid-scale channel feedbacks by treating τ * c(q s ) as a state variable (i.e., Eq. 10).Conceptually, the τ * c(q s ) model "lumps" processes related to multiple bed state variables (Sects.2.1, 3.1).Similarly, because many channel state variables influence transport and therefore are not independent of τ * c(q s ) , I hypothesize that aspects of morphology can be implicitly subsumed into evolving τ * c(q s ) for modeling purposes, because τ * c(q s ) captures essential feedbacks over spatial and temporal scales of interest.This is similar to the channelization approach of Stark and Stark (2001).

Form drag vs. parsimony
Calculations of best-fit τ * rs50 and transport rates used total shear stress (Eq.12), rather than partitioning stress into form drag and a lower effective stress for calculating transport rates (skin friction).Although not a state variable, form drag is physically justifiable because larger clasts that protrude higher into the flow (e.g., stable boulders) tend to account for a disproportionate amount of the total stress through drag, turbulence generation and pressure gradients.Form drag corrections have been incorporated into many transport models to enable reasonable transport rates to be calculated using τ * c values typical of systems without form drag (e.g., Rickenmann and Recking, 2011;David et al., 2011;Yager et al., 2012a).Conversely, another common approach (and that taken here) is simply to use higher τ * c (e.g., Bunte et al., 2013;Lenzi et al., 2006), treating τ * c as a physically meaningful fitting parameter to predict transport.Using field data, Schneider et al. (2015) recently compared gravel transport predictions based on (a) form drag corrections and (b) higher reference stresses.For the most part, they found that both approaches could provide similar accuracy.They also noted that "uncertainties in predicted transport rates remain huge (up to roughly 3 orders of magnitude)" (Schneider et al., 2015), and suggested that factors including supply effects may account for remaining discrepancies.Although beyond the scope of the present analysis, form drag effects could be separated from best-fit τ * rs50 by using a calculated skin friction stress rather than total stress.However, doing so would add extra uncertainty to the shear stresses, while still not directly accounting for effects of sediment supply.Implicitly subsuming form drag into τ * c(q s ) arguably provides a simpler and more parsimonious approach for modeling transport and morphodynamics.

Conclusions
I propose a new model in which feedback causes τ * c(q s ) , the nondimensional critical shear stress for gravel transport, to evolve through time as a function of sediment transport disequilibrium (Eqs. 8,10).Net erosion tends to increase local τ * c(q s ) (reducing transport rates), while net deposition tends to decrease τ * c(q s ) (increasing transport rates).Laboratory flume experiments described by Johnson et al. (2015) are used to evaluate the proposed τ * c(q s ) model.The experiments intentionally explored disequilibrium bedload transport and morphodynamic adjustment.Thresholds of motion were backcalculated from the experimental data using the Wilcock and Crowe (2003) model for mixed grain size transport.I also show that the Wilcock and Crowe (2003) hiding function is consistent with our experimental data, supporting its applicability to steep channels.
After empirically calibrating three model parameters, the τ * c(q s ) model -a differential equation -can explain nearly 70 % of the variability in experimental thresholds of motion.I then incorporate τ * c(q s ) into a simple morphodynamic model for channel profile evolution.Changes in τ * c(q s ) are negative feedbacks on morphodynamic response, because not only slope but also τ * c(q s ) evolve when perturbed.Finally, τ * c(q s ) is redefined to be a state variable for fluvial channels.State and path variables are fundamental to many disciplines such as thermodynamics, because they allow the evolution of systems to be calculated.The same should be true for morphodynamics.Conceptualizing landscape evolution models in terms of feedbacks among evolving state variables and path functions may improve our ability to predict landscape responses to land use, climate change and tectonic forcing.

Data availability
Data sources for Fig. 1 are given in its caption.An additional data table (data set "grl53586-sup-0002-supplementary.xls") that is provided with Johnson et al. (2015) contains data on Fig. 2. Data and matlab code to make figures in the paper based on experimental calibrations are available through Figshare (doi:10.6084/m9.figshare.3635298.v1).

*Figure 3
Figure 3. τ *rs50 fits to the experimental data with the W&CM."W&CM fit" uses Eq. (16) to calculate hiding function exponent b, while "Power-law fit" calculates a best-fit b along with τ * rs50 .Error bars give 95 % confidence intervals on τ * rs50 based on the regressions; although uncertainty can be broad the trends are clear and consistent.Shaded area indicates times of fine gravel addition (sediment feed) in each experiment.Data in these plots are available inJohnson (2016).

Figure 4 .
Figure 4. Data points are based on power-law fits to exponent b.The W&CM hiding function (Eq.16) does a good job matching the data, although it was not fit to these points.The first six measurements of each experiment (roughly the first 10 min) were excluded because of large scatter associated with the greatest bed instability.The axes reflect the left and right hand sides of Eq. (15), but uses dimensional stresses to be consistent with plots shown inWilcock and Crowe (2003).

Figure 5 .
Figure 5. Best-fit models (Eqs.4, 8 and 10) compared to experimental constraints.The periods of upstream sediment supply (Q sfeed ) are indicated by the grey boxes for each experiment.
Figure6.Profile evolution, comparing the morphodynamic responses of models with and without threshold evolution.The initial condition is an equilibrium channel with τ * c(q s ) = 0.06, upstream sediment supply q s = 1 × 10 −3 m 2 s −1 , and an initial equilibrium slope of 0.0147.Sediment supply is increased 5× at t = 0. Lines are each 5 model days apart, and indicate the evolution to a new transport equilibrium.

Figure 7 .
Figure 7. Slope and critical shear stress evolution, for sediment supply increases (which correspond to Fig. 6 models) and decreases by factors of 5.As in Fig.6, t = 0 corresponds to an equilibrium condition where the initial slope and initial threshold are consistent with the initial upstream sediment supply.Slope and τ * c(q s ) were averaged over nodes 3-10, leaving out the first and last two nodes because of minor model boundary effects.

Figure 8 .
Figure 8. Morphodynamic model sensitivity to sediment supply perturbations and k.All models at the same equilibrium condition as shown in Figs. 6 and 7. (a) Slope adjustment, normalized by the initial equilibrium slope.The correspondence of Eq. (20) and the morphodynamic model calculations demonstrate that the models did asymptotically attain equilibrium slopes.(b) S eq ratio is the ratio of equilibrium slopes of the Exner +τ * c(q s ) model divided by S eq for the Exner-only model, to show the relative affect that τ * c evolution has on equilibrium slopes.(c) Equilibrium timescales for model adjustment.(d) t eq ratio is the ratio of t eq for the Exner +τ * c(q s ) model divided by t eq for the Exner-only model.Values are lower than 1, indicating that the τ * c evolution has a large influence on equilibrium timescales.(e) Evolution of τ * c(q s ) .

Figure 9 .
Figure 9. Spatial and temporal evolution of morphodynamic slopes, for the same models shown in Fig. 6.Slope is initially at equilibrium and responds to the 5× increase in upstream sediment supply at t = 0. (a) The Exner-only model initially has spatial slope variability, but evolves to a uniform new equilibrium slope.(b, c) In the Exner +τ * c(q s ) model, spatial variability in both slope and τ * c(q s )

Figure 10 .
Figure 10.Comparison of experimental and best-fit model constraints on τ * rs50 , compared to proposed constraints for D 84 reference stress bounds for low and high sediment supply from Recking (2012).

Table 1 .
Best-fit threshold evolution models.