Upscaling Sediment-Flux-Dependent Fluvial Bedrock Incision to Long Timescales

Fluvial bedrock incision is driven by the impact of moving bedload particles. Mechanistic, sediment-flux-dependent incision models have been proposed, but the stream power incision model (SPIM) is frequently used to model landscape evolution over large spatial and temporal scales. This disconnect between the mechanistic understanding of fluvial bedrock incision on the process scale, and the way it is modeled on long time scales presents one of the current challenges in quantitative geomorphology. Here, a mechanistic model of fluvial bedrock incision that is rooted in current process understanding is explicitly upscaled to long time scales by integrating over the distribution of discharge. The model predicts a channel long profile form equivalent to the one yielded by the SPIM

observations from laboratory experiments and natural streams have by now been reported, demonstrating that bedload transport exerts a dominant control on the patterns and rates of erosion (e.g., Beer et al., 2017;Finnegan et al., 2007;Mishra et al., 2018;Shepherd, 1972;Turowski, Hovius, Hsieh, et al., 2008;Wohl & Ikeda, 1997).A number of sediment-related effects have been identified (e.g., Sklar & Dietrich, 2001;Turowski et al., 2013), two of which seem to be most important.The tools effect arises because fluvial bedrock erosion is driven by the impacts of moving bedload particles, implying that an increasing number of moving particles leads to an increasing number of impacts and therefore higher erosion rates (e.g., Beer & Turowski, 2015;Cook et al., 2013;Foley, 1980;Inoue et al., 2014).The cover effect arises because sediment residing on the bed can shield the bedrock from impacts, thereby decreasing erosion rates (e.g., Chatanantavet & Parker, 2008;Mishra & Inoue, 2020;Turowski, Hovius, Hsieh, et al., 2008).Yet, in landscape evolution models designed for long timescales, fluvial bedrock erosion is commonly described by the stream power incision model (SPIM), in which incision rate is a power function of water discharge and channel bed slope (e.g., Barnhart et al., 2020;Seidl & Dietrich, 1992).The SPIM is unable to account for the tools and cover effects on the process to decadal time scales, but it remains popular because of its simple form.In addition, it explains the widely observed power law scaling between channel bed slope and drainage area, and spatial patterns of knickpoint migration speeds (see Lague, 2014, for a summary of field evidence in the context of the SPIM).
The gap between mechanistic processes understanding on short timescales, and the popularity of the SPIM on long timescales currently represents a central challenge in the study of bedrock channel morphodynamics (Venditti et al., 2020).Diverse temporal scales can be theoretically connected by explicit upscaling, integrating instantaneous process descriptions over the distributions of forcing variables (e.g., Blom et al., 2017;Lague et al., 2005).Attempts to upscale sediment-flux-dependent incision models in this way have so far been scarce, because multiple interacting variables make analytical solutions challenging.Turowski et al. (2007) partitioned sediment-carrying and clean flows using a method suggested by Sklar and Dietrich (2006) in an analytical model of bedrock channel morphology including both tools and cover effects.Lague (2010) included the cover effect, but not the tools effect, into a numerical model of bedrock channel evolution, forced by random time series of daily discharge following an inverse gamma distribution (Crave & Davy, 2001).None of these attempts captures the entire range of conditions and dynamic behavior that can be expected for natural bedrock rivers.
Here, I present analytical solutions for the long-term incision rate and steady state channel morphology using a mechanistic incision law including both tools and cover effects.The solutions demonstrate that the steady state channel long profile is set by bedload transport rather than bedrock incision processes, and offers insights into the role of thresholds and channel width, and the river's adjustment to variable discharge.

Theoretical Treatment
In this section, I develop a description of a steady state bedrock channel, upscaling from a sediment-flux-dependent erosion law including both tools and cover effects.Several stochastically varying forcing variables, including water discharge and bedload transport rate, and dependent variables such as bed cover that may exhibit a strong history dependence, are addressed in turn to explain the assumptions made to make an analytical solution possible.By the end, solutions for the long-term mean sediment transport rate, bed cover and incision rate are obtained.

General Consideration
In previous treatments, water discharge was assumed to be the only stochastically fluctuating parameter.In this case, to upscale instantaneous incision laws to long time scales, we need to integrate over the distribution of discharge, assumed to follow the inverse gamma distribution (e.g., Crave & Davy, 2001;Lague et al., 2005;Molnar et al., 2006) Here,   / Q Q Q is the instantaneous discharge Q normalized by the long-term mean discharge Q, the constant k is a measure of the variability of discharge.For natural streams, k takes typical values between 0.1 and 5 (Molnar et al., 2006).Note that k decreases with increasing discharge variability, that is, k = 0.5 signifies a channel with highly variable discharge and k = 3 signifies a channel with little discharge variability.exp{x} denotes the natural exponential function, and Γ(x) denotes the gamma function, defined by Here, z is a dummy variable.The long-term mean of a particular discharge-dependent quantity X of interest can be obtained by integrating over the distribution Here, Q min and Q max denote the minimum and maximum discharge considered to be relevant for setting X, and the overbar denotes the long-term mean of a parameter, as obtained by the integral in Equation 3.For the analysis in the paper, I assume that Q max is sufficiently high such that the distribution of discharge is adequately captured by integrating the right-hand power law tail of the discharge distribution to infinity (see Lague et al., 2005, for a detailed discussion of the effects of this assumption).Then, Equation 3becomes When dealing with sediment-flux-dependent incision laws, the bedload transport rate is another driving variable affecting incision rates directly via the tools effect and indirectly via the cover effect.Bedload transport rates can fluctuate strongly, and measured rates can scatter over several orders of magnitude for a given discharge (e.g., Turowski, 2010).In addition, the amount of sediment residing on the bed, which determines bed cover, is a history-dependent state variable.Integrating explicitly over the temporal variation of these variables would prevent an analytical treatment and necessitate a numerical solution.To deal with this problem, I introduce an intermediate timescale.At this timescale, the short-term fluctuations of bedload transport rates and sediment cover are averaged out, and the average can be treated as a deterministic function of discharge.To clearly distinguish the quantities at the two timescales, I use the term "average" and angle brackets <> for the intermediate timescale, and the term "mean" and an overbar for the geological timescale.

Treatment of Channel Width
A fully dynamic model of bedrock channel width in a sediment-flux-dependent setting is currently not available.Commonly, channel width is assumed to depend on discharge according to a power law, using standard downstream and at-a-station hydraulic geometry relationships of the form (e.g., Lague et al., 2005 Here, W is the channel width corresponding to the long-term mean discharge at a particular station, W is the instantaneous width varying locally with discharge, and ω d and ω a are dimensionless exponents.Within the present treatment, I replace Equation 5 with the steady state width equation obtained from the model of Turowski (2018) TUROWSKI 10.1029/2020JF005880 3 of 20 Here, k e is a measure of the bedrock erodibility, Q s is the bedload transport rate and I the incision rate, and the sideward deflection distance d is the distance by which bedload particles can be deflected in the cross-channel direction (Turowski, 2018).Here, d is treated as a constant, which can be viewed as a general scaling factor with unit of length within the context of long-term channel morphology.

Upscaling Bedload Transport
Bedload transport rate can fluctuate strongly even if hydraulic conditions stay constant over time (e.g., Turowski, 2010).However, at a given discharge, there exists a well-defined mean transport rate that scales with discharge.At the intermediate timescale, I assume short-term fluctuations of transport rates can be neglected, and average bedload supply at a given discharge is a function of discharge and slope of the form (e.g., Rickenmann, 2001;Smith & Bretherton, 1972) Here, k BL is a dimensional constant, S is the channel bed slope, Q ct is the critical discharge for the onset of bedload transport m, n, and q are dimensionless exponents, and the angle brackets denote the average quantity at a given discharge at the intermediate time scale.Many standard sediment transport formulas can be expressed in the form of Equation 8.The power function of channel width W is included to make possible the modeling of the varying scaling of sediment flux with channel width (Carson & Griffiths, 1987;Cook et al., 2020), depending on the sediment transport and flow velocity equations that are used (Figure 1, Appendix A).Note that, depending on the choice of bedload transport equation, m, n, and q are not independent of each other (Appendix A).Using the normalized discharge Q*, the bedload transport rate can be rewritten as Combining Equation 9with Equations 6 and 7, we obtain Using Equation 4, the long-term sediment flux is then given by The integral evaluates to Here,  Appendix A).Note that the threshold of motion cuts off the relationship both for large width, because flow depth becomes too small for transport to occur, and for small width, because shear stress is partitioned from the bed to the channel walls (e.g., Turowski, Hovius, Hsieh, et al., 2008).The scaling bracketing the relationship for small width, giving q = 5/2 (black line), for intermediate width using the Darcy-Weissbach friction equation, giving q = 0 (dashed red line), and the Manning equation, giving q = 1/10 (solid red line), are indicated (see Appendix A).The plot was generated using Equations A1, A8, A9, and A10 (Appendix A), with    0.045 The upper incomplete gamma function is defined by x c z z z (14)

Upscaling Bed Cover
Bed cover C can vary over short timescales, and is dependent on the history of sediment supply and hydraulic forcing (e.g., Fernández et al., 2019;Lague, 2010;Turowski & Hodge, 2017).However, response timescales of bed cover to varying flow conditions are orders of magnitudes smaller than those of the adjustment of channel width and slope (Turowski, 2020).As a result, similar to bedload transport, cover can be treated to be independent of discharge at the intermediate timescale, following a distribution with a well-defined average for a given discharge.This implies that instantaneous cover C can be viewed as independent of discharge, and the intermediate-term average cover <C > systematically varies with discharge.Here, the relationship between the average cover and discharge is modeled by a power law function with a scaling exponent α (Turowski et al., 2013), from hereon called the cover exponent (Figure 2).The bed changes from fully to partially covered at a characteristic dimensionless discharge Q .When α > 0, bed cover increases with increasing discharge, and the bedrock is exposed during small flows, for . This is the flood-depositing case, for which the cover function is given by When α < 0, bed cover decreases with increasing discharge, and bedrock is exposed during large flows, for . This is the flood-cleaning case, for which the cover function is given by For convenience, the cover threshold can be written as a multiple b of the threshold discharge  ct Q for the onset of bedload transport.The ratio b is defined as With the assumptions made so far, the system features two discharge thresholds, one for the onset of bedload motion and therefore the activity of the tools effect, the other for the change of fully to partially covered TUROWSKI 10.1029/2020JF005880 5 of 20  Turowski et al., 2013).For the flood-depositing case, where cover exponent α > 0, the covered fraction of the bed increases with increasing discharge, implying that the bed is partially covered at discharge smaller than the characteristic discharge Q and fully covered at discharges above it.Bedrock erosion occurs at low and intermediate discharges.For the flood-cleaning case, where the cover exponent α < 0, the covered fraction of the bed decreases with increasing discharge, implying that the bed is partially covered at discharge larger than the characteristic discharge  cc Q and fully covered at discharges below it.
Bedrock erosion occurs during floods.
bed.Together with the two types of cover behavior of the channel (Figure 2), flood-depositing (α > 0, Equation 15) and flood-cleaning (α < 0, Equation 16), we can distinguish four cases, yielding different integrative limits and solutions for the long-term results (Table 1).In the flood-cleaning case, when , erosion occurs for all discharges greater than ), erosion occurs for all discharges greater than  ct Q .In the flood-depositing case, when , erosion occurs for all discharges greater than  ct Q and smaller than ), no erosion occurs, because bedload moves and tools are only available at discharges when the bed is fully covered.For the three cases in which erosion occurs at some discharges, Equation 4 can be applied to calculate the long-term mean cover.The solutions have the general form Here, F C is a dimensionless function.Full solutions for F C are given in Appendix B.

Upscaling Sediment-Flux-Dependent Incision
A sediment-flux-dependent erosion law, including tools and cover effects, can be given by (Auel et al., 2017;Sklar & Dietrich, 2004;Turowski, 2018;Zhang et al., 2015): The dimensional constant k e depends on rock, sediment, and fluid properties, given by Auel et al. (2017) as Here, g is the acceleration due to gravity, Y is Young's modulus of the bedrock, σ T its tensile strength, k v is the dimensionless rock resistance coefficient, and ρ s and ρ are the sediment and fluid density, respectively.At the intermediate timescale, Equation 19 can be rewritten as Combining Equations 4,6,7,10,12,13,17,18,and 21, the long-term incision rate can be evaluated by the integral Here, the limits of integration Q depend on the values of α and b (see Table 1).The full solutions for all three cases are given in Appendix C, and take the general form , , , , , , .
, , , , , , Threshold of motion larger than cover threshold a Bedload transport occurs and tools are available only at large discharges, when the bed is fully covered.

Table 1
The Ranges of Q* for Which Erosion is Possible in the Four Cases

Results
In general, there are four unknown variables, channel bed slope S, the long-term cover fraction C (Equation 18; see also Appendix B, Equations B3, B6, and B9), the long-term bedload sediment supply s Q (Equation 12), and the ratio between the cover threshold and the threshold of bedload motion b (Equation 16).The solutions provide three equations.The long-term incision rate I (Equation 23) can be treated as an independent variable that is determined by the long-term uplift or baselevel lowering rate.Another equation can be obtained from the conditions for steady state, when the long-term bedload supply is related to the longterm incision rate by Here,  is the long-term mean fraction of sediment that is transported as bedload, and A is the drainage area.I further substitute discharge with a simple hydrologic relation Here, R is the long-term mean runoff.To illustrate the dependence of channel morphology and of the adjustment time scales on control and channel morphology parameters, I used parameter values oriented on Lushui at the Liwu River, Taiwan (Table 2; see Turowski et al., 2007, andTurowski, 2020).The values of reach parameters were either measured in the field or estimated using literature data.

Steady State Channel Long Profile
Both long-term bedload supply (Equation 12) and long-term incision rate (Equation 23) show the same dependence on channel bed slope, while long-term mean cover is independent of slope (Appendix B).As a result, channel bed slope S can be calculated from the equation for long-term bedload supply (Equation 12), which is independent of the cover threshold and of long-term cover.Inverting Equation 12for S and substituting Equations 24 and 25 yields For certain combinations of the parameter values, the function F Qs may be negative or not give a solution at all.In these cases, Equation 26does not yield a valid solution for the channel bed slope.Parameter combinations without solutions occur mainly for 1 < k < 2.5 (Figure 3).

Scaling With Discharge Variability
The controls on channel morphology by discharge variability k are complicated (Figure 4), and depend both on the cover scaling exponent α and the width scaling exponent q.For the same discharge variability k, multiple possible solutions are available for most of the parameter space.The solution for the channel bed slope S is independent of α, but strongly dependent on q (Figure 4a).For small values of k, S increases with increasing k, for intermediate values of k, no solutions are available (see also Figure 3), and for large values of k, S decreases with increasing k (Figure 4a).The ratio of cover threshold to threshold of motion b (Equation 17), and the long-term mean cover C can increase or decrease with TUROWSKI 10.1029/2020JF005880  (Turowski, 2020;Turowski et al., 2007) increasing discharge variability k, depending on the values of q and α (Figures 4b and 4c).The threshold ratio b is weakly dependent on k for k < 2 and k > 3.5 (Figure 4b).For 2 < k < 3.5, it decreases strongly with increasing k when α = 2, and increases weakly when α = −2.The longterm mean cover is constant for k < 1.5, regardless of the values of q and α (Figure 4c).For k > 2.5, it increases when α = −2, but it decreases when α = 2.

Control of Reach-Scale Cover Behavior
The ratio of cover threshold to threshold of motion b (Equation 17), and the long-term mean cover, depends on reach-scale cover behavior, that is, whether the channel behaves as flood-cleaning (when the cover scaling exponent α < 0) or flood-depositing (α > 0).For flood-cleaning channels, b is close to zero or to one, for flood-depositing channels, it is larger than one (Figure 5a & 5c).For streams with low discharge variability (high k) the long-term mean cover increases for increasing α for flood-cleaning channels, and decreases for increasing α for flood-depositing channels (Figure 5b).For streams with high discharge variability, the function is complicated (Figure 5d)

Steady State Channel Long Profile
Empirically, the channel long profile of bedrock rivers is often described by a power law function, as has been observed in many natural settings (e.g., Whipple, 2004;Whitbread et al., 2015): Here, k s is known as the steepness index and θ is known as the concavity index.While the concavity index θ typically falls into a narrow range between 0.4 and 0.7 in natural bedrock channels (e.g., Lague, 2014;Whipple, 2004), the value of the steepness index k s can vary over several orders of magnitude (e.g., Barnhart et al., 2020).The upscaled model yields a similar slope-area scaling (Equation 26; Figure 6), in which the steady state channel long-profile is controlled by the mechanics of bedload transport, rather than the mechanics of bedrock incision.This notion is consistent with field observations of Johnson et al. (2009).Assuming that the long-term average bedload fraction  scales with drainage area    17), and (c) of the long-term mean cover C with the discharge variability parameter k.Lines for the cover exponent α = −2 (flood-cleaning), crosses for α = 2 (flood-depositing).Black for the width exponent q = 0, blue for q = 1, and red for q = 2.5.Note that channel bed slope is independent of α (Equation 26).For many conditions, there are two solutions available, corresponding to the solutions for b smaller (dashed) or larger than one (solid) (see Table 1; Appendices B and C).

A) B ) C )
The concavity index θ is then given by    26) and (b) corresponding elevation profiles.The concavity index θ was calculated using Equation 29.Relationships are shown for m = 2 and n = 1, and B = q = 0 (black line), B = 0 and q = 1 (blue line), B = 0 and q = 2.5 (red line), and B = 1 (black dashed).The steepness index was calibrated to the conditions at the Liwu (red circle, Table 2), with A = 435 km 2 , S = 0.02 ( The bedload fraction typically decreases with increasing drainage when different river catchments are compared (Turowski et al., 2010).Based on field observations in the Himalayas, Dingle et al. (2017) suggested that the bedload transport rate in actively eroding rivers is constant along a river, despite increasing drainage area, implying B = 1 for a steady state catchment (cf.Equation 24).In the following, I will discuss the endmember cases of B = 1 (bedload transport rate independent of drainage area) and B = 0 (bedload fraction independent of drainage area).In the former case, for B = 1, the concavity index is equal to the ratio of the discharge and slope exponents in the bedload transport Equation 8, m/n, independent of the value of q.In the latter case, for B = 0, a wide range of values can be obtained for concavity index, depending on the choices for m, n, and q.In natural rivers, θ is usual within the range between 0.4 and 0.7 (Lague, 2014;Whipple, 2004).A value of θ = 5/8 = 0.625 is obtained for m = 1 and n = 2, as in the bedload transport equation of Rickenmann (2001) (see also the discussion of Turowski, 2018), and q = 5/2, corresponding to the limit behavior for narrow channels (Figure 1; Appendix A).Smaller values of q decrease the concavity, with θ = 0 for q = 0, and θ = 1/4 = 0.25 for q = 1.A value of θ = 0.5, a standard choice in many modeling exercises, is obtained for q = 2.
The channel long-profile solution of the upscaled model (Equation 26) gives an explicit expression for the steepness index in terms of long-term incision rate I , erodibility k e (Equation 20), long-term mean bedload fraction  , long-term mean runoff , R as well as channel geometry parameters (e.g., the width exponent q) and discharge variability k: Here, F Qs is given by Equation 13.The steepness index k s can in principle be calibrated to measurable parameters (Appendix D).The formulation helps explaining the large variability observed for the steepness index in nature.For example, erodibility depends on the inverse square of rock tensile strength (Equation 20), which varies over more than 2 orders of magnitude for natural rock (e.g., Sklar & Dietrich, 2001).
In addition, Equation 30 allows stream profile inversions for incision rate (e.g., Wobus et al., 2006) with considerably more detail than the commonly used SPIM.However, many of the necessary parameters, such as the width exponent q, are not usually known for natural channels, and would require extensive field measurements.As a result, no data sets are currently available that include all necessary parameter values for natural channels.Therefore, a full test of the model is currently not possible.
The lack of valid solutions for channel bed slope for certain parameter combinations (Figure 4) occurs when the sum of four terms including the gamma or incomplete gamma functions in Equation 13 are equal to or smaller than zero.The values of the four terms depend on discharge variability k, the discharge exponent m in the bedload transport equation (Equation 8), and the product of the width exponent q (Equation 8) and the at-a-station hydraulic geometry exponent for width, ω a (Equation 6).This suggests that the river needs to adjust its absolute width (which changes the width exponent q) and its at-a-station hydraulic geometry for width (i.e., the cross-sectional shape, which changes ω a ) to achieve a channel long-profile that is consistent with the condition of grade, in which sediment deposition and entrainment are balanced.The results underline the importance of channel width for understanding bedrock channel dynamics.

Cover and Thresholds
According to the model, long-term channel dynamics are controlled by at least two discharge thresholds, the critical discharge for the onset of bedload motion, and the discharge at which the channel switches from a fully to a partially covered bed.The relationship between these two threshold discharges, quantified in their ratio b, depends strongly both on the reach-scale cover behavior (quantified by the cover scaling exponent α; see Figure 5a) and discharge variability k (Figure 5b).Often multiple solutions for the channel geometry are available for the same set of boundary conditions (Figure 4).Turowski et al. (2013) showed that both flood-cleaning (α < 0) and flood-depositing (α > 0) bedrock channels exist in nature, and sometimes reaches behaving one way or the other alternate in a single stream (e.g., Heritage et al., 2004).It is unclear what controls the cover scaling exponent α, and correspondingly, why a particular reach or stream behaves flood-cleaning or flood-depositing.It can be expected that both hillslope processes and in-channel processes contribute to this control (cf.Turowski et al., 2013).For example, precipitation amount and intensity control hillslope sediment supply to the channel by landsliding or surface wash, but also the in-channel sediment transport capacity via their relationship to discharge.Water discharge, in turn, affects upstream sediment supply to a given reach, as well as bank and bed erosion rates.The cover exponent α could be related to or depend on other parameters such as discharge variability k, and the cross-sectional shape, quantified by the at-a-station hydraulic geometry exponent for width, ω a , and the width exponent q.All four of these parameters are treated as independent variables in the model, but may adjust their values through yet unknown feedback mechanisms.These topics provide starting points for future research.The parameters ω a , k, and α can be obtained from field observations.The hydraulic geometry exponent ω a can be constrained from parallel measurements of width and discharge (e.g., Turowski, Hovius, Wilson, & Horng, 2008), and the cover exponent α from parallel observations of cover and discharge.To estimate the discharge variability parameter k, a discharge time series is necessary spanning at least a few years, better decades (e.g., Molnar et al., 2006).The deflection length scale d is likely difficult to measure in natural settings, but should be well suited for study in laboratory experiments.
The combination of a flood-depositing channel with a threshold discharge for the onset of bedload motion that is higher than the discharge at which the bed becomes fully alluviated allows for a solution in which no incision occurs.Bedrock channels may evolve to this state when tectonic activity of a mountain belt ceases and uplift stops, rather than turning into an alluvial channel.

Upscaling Instantaneous and Intermediate Timescale Process Descriptions
In the description of long-term evolution of river channels or entire landscapes, it is common to upscale instantaneous or intermediate timescale process descriptions by substituting a representative long-term discharge for the instantaneous or intermediate discharge.The upscaled version for bedload transport (Equation 12) resembles the process description at intermediate timescales (Equation 8), similar to the upscaled version of the SPIM (Lague et al., 2005).There are three key differences: (i) instantaneous and intermediate timescale variables are replaced by their long-term equivalents, (ii) the threshold vanishes in the long-term equation, and (iii) a dimensionless function, dependent on discharge variability k, channel morphology parameters, ω a , and q, and thresholds, such as the threshold of sediment motion Q ct , is multiplied on (e.g., Equation 13).In contrast, the upscaled equations for bed cover (Equation 18) does not resemble the intermediate timescale version (Equations 15 and 16).In comparison to the instantaneous incision rate (Equation 19), in the equation for the upscaled incision rate (Equation 23) the dependence on cover is replaced by a different dimensionless function (Appendix C).As such, while equations similar to the instantaneous process law can be used in the long term, the determination of the dimensionless scaling factor may be complicated, and depend, for example, on the discharge variability parameter k.Both channel bed slope S (Figure 4a) and the threshold ratio b (Figure 4b) become independent of the discharge variability parameter k for large k, corresponding to small discharge variability.Consistent with previous notions (cf.Deal et al., 2018;Lague et al., 2005), this may imply that discharge variability does not need to be precisely known for environments with high k (low variability), when channel slope S is the primary interest.However, the long-term cover C behaves in the opposite way (Figure 4c), with constant values for small k and high variability for large k.Further, the relationship between the threshold ratio b and the cover exponent α are strongly modulated by k (Figures 5a and 5c), and slope depends on q even for large k.Again, questions about the feedbacks between α, the discharge variability parameter k, the cross-sectional shape, quantified by the at-a-station hydraulic geometry exponent for width, ω a , and the width exponent q arise and provide starting points for future research.

Comparison to the Upscaled Stream-Power Incision Model
The upscaled SPIM (Appendix E; Lague et al., 2005) is able to capture the nonlinear dependence of incision rates on discharge variability observed in the Himalaya (DiBiase & Whipple, 2011;Scherler et al., 2017).The results obtained from sediment-flux-dependent incision models give similar relationships for flood-cleaning channels and a ratio of cover threshold to threshold of motion b that is smaller than one (Figure 7).26).The SPIM has been claimed to provide a description of bedrock channel dynamics on long time scales (e.g., Venditti et al., 2020), even though its mechanistic assumption-a direct scaling between erosion rate and stream power (Seidl & Dietrich, 1992)-has been falsified on the process time scale (e.g., Beer & Turowski, 2015;Sklar & Dietrich, 2001).Whipple and Tucker (2002) already recognized that a wide range of incision models yield similar or even identical predictions for the channel long profile, and concluded that observations of transient dynamics need to be used to assess model efficacy.However, studies that have attempted this arrived at conflicting results.For example, van der Beek and Bishop (2003) found that all of the tested models could be parameterized to explain their observations, while Tomkin et al. (2003) concluded that none of the tested models could be fit to their data with physically meaningful parameter values.Valla et al. (2010) found that a transport-limited model better described their data, while Attal et al. ( 2011) argued that a SPIM yielded the best fit, provided a threshold of erosion was included.Even though it is limited to steady state channels, the model developed herein yields multiple possible solution for a given set of boundary conditions.It has previously been shown that sediment-flux-dependent incision models can yield transient behavior that mimics either transport-or detachment-limited conditions or a mixture of both (e.g., Davy & Lague, 2009;Gasparini et al., 2007;Whipple & Tucker, 2002).This suggests that sediment-flux-dependent incision models as used here can yield the rich transient behavior inferred from observations in natural bedrock channels (cf.Lague, 2014).

Conclusions
Bedrock channels are thought to evolve toward a steady state in which the long-term incision rate is equal to the long-term baselevel-lowering rate (e.g., Lague et al., 2005;Whipple, 2004), but also have the need to transport the supplied sediment load (e.g., Turowski, 2020).Upscaling a sediment-flux-dependent incision model rooted in current mechanistic understanding of fluvial bedrock incision processes yields several solutions that are consistent with the requirements for the long-term steady state, but differ in their short time dynamics, for example in the relationship of bed cover and discharge, or transport capacity and channel width.Some of these solutions are similar to the behavior expected in the stream power paradigm, but other solutions are also possible.In the stream-power incision model, both in its instantaneous and in its upscaled form, the effects of channel width, of sediment supply and transport, and of bedrock erodibility TUROWSKI 10.1029/2020JF005880 12 of 20 are lumped together in a single calibration parameter, usually called the erodibility.In contrast, the model presented here makes the relationship of the long-term incision rate and of channel geometry with these effects explicit.It thus advances a more detailed picture into the long-term behavior of bedrock channels and offers a wider range of testable predictions and assumptions.
For certain parameter combinations, the model does not yield a solution for channel bed slope.However, it is in most cases possible to find a solution by adjusting the dependence of the bedload transport capacity on channel width (via the exponent q), the cross-sectional geometry (via the at-a-station hydraulic geometry exponent for width, ω a ), or the reach-scale cover dynamics (via the exponent α).It is unclear what controls these parameters in nature and whether they interdepend on each other.Still, the model results suggest that the long-term dynamic behavior of bedrock channels is richer than previously thought, and that there are controls and feedbacks that have been little explored so far.
The long-term incision rate (Equation 26) is explicitly dependent on bedrock erodibility, channel width, and discharge variability, as well as the free parameters determining reach-scale cover behavior and cross-sectional shape.All of these effects have previously been argued to be important factors in setting incision rates (e.g., Bursztyn et al., 2015;Cook et al., 2013;Lague et al., 2005;Whitbread et al., 2015).However, within the SPIM, they are lumped into a single calibration parameter.The new formulation makes it possible to separate all of these effects.This offers rich new possibilities for testing the model using data from natural streams.

Appendix A: Width Dependence of Bedload Transport Rate
A commonly used equation for bedload transport has the form (e.g., Fernandez Luque & van Beek, 1976;Meyer-Peter & Müller, 1948) Here, γ is a dimensionless coefficient and D is a representative grain size.The reach-averaged Shields stress   is defined by Here, τ is the shear stress, given by the DuBoys equation Here, R h is the hydraulic radius.The continuity equation for water flow is Here, A c is the cross-sectional area.For a rectangular channel, hydraulic radius and cross-sectional area are given by Here, H is the flow depth.A generic form for the cross-section averaged flow velocity can be written as Here, k V is a friction coefficient and δ takes the value of 1/2 for the Darcy-Weissbach equation and 2/3 for the Manning equation.Eliminating H, R h , A c , and V by combining Equations A3-A7, we obtain Equation A8 does not permit a closed-form solution for τ (see Figure 1 for a numerical solution).However, we can make some statements on scaling.If width is small, the quadratic term in width can be neglected.Then For intermediate width, the term independent of width on the left-hand side can be neglected, and scales as Applying Equations A10 or A11 to the bedload transport Equation A1, we can obtain sensible values for the width exponent q in Equation 8: Using Equation A10, shear stress τ is independent of discharge (m = 0), and proportional to both slope and width, and n = 3/2 and q = 5/2.Using Equation A11, we need to distinguish for the two friction equations.
For the Darcy-Weissbach equation, when δ = 1/2, this results in m = n = 1 and q = 0, while for the Manning equation, when δ = 2/3, this results in m = 9/10, n = 21/20 and q = 1/10.Note that Rickenmann (2001) showed that n = 2 gives a better description of field data than the theoretical values given above.He suggested that this deviation results from the decreasing occurrence and importance of macro-roughness elements such as stationary boulders or step-pool sequences when moving downstream.

Appendix B: Solutions for Long-Term Average Cover
For the flood-cleaning case (α < 0) with Combining with Equation 4, the integral becomes The integral evaluates to For the flood-cleaning case (α < 0) with The integral becomes For the flood-depositing case (α > 0) with The integral becomes

Appendix C: Solutions for the Long-Term Incision Rate
For the flood-cleaning case (α < 0) with For the flood-cleaning case (α < 0) with For the flood-depositing case (α > 0) with , we obtain (C3)

Appendix D: Calibrating the Model Using Field Observations of Slope and Width
Channel width and slope as well as drainage area and discharge can be measured in the field with commonly available methods.Under the assumption that the channel is at a long-term steady state, such measurements can be used to calibrate various free parameters in the equations.From Equations 7 and 24, channel width is given by With measurements of slope, width, drainage area, mean discharge, and long-term incision rate, and an estimate of the bedload fraction, Equation D2 can be used to calculate k BL .Subsequently, erodibility k e can be calculated using Equations D4 and D1 can be used to calculate the deflection length scale d.Lague et al. (2005) gave a comprehensive discussion of the upscaled stream power incision model (SPIM).In a stream with variable discharge, the bedrock incision rate according to the SPIM is given by (cf.Lague et al., 2005)  

Appendix E: Upscaling the Stream Power Incision Model
The integral is given by The long-term incision rate is then given by

Figure 2 .
Figure2.Illustration of the scaling relationship of cover with discharge and the definitions of flood-cleaning and flood-depositing channels (adapted fromTurowski et al., 2013).For the flood-depositing case, where cover exponent α > 0, the covered fraction of the bed increases with increasing discharge, implying that the bed is partially covered at discharge smaller than the characteristic discharge function depending on the values of α and b (Appendix C).

Figure 3 .
Figure 3. Space of valid solutions of slope (Equation 26) are shown in gray, for a thresholds of motion   ct 1 Q .Other choices for  ct Q only slightly

Figure 4 .
Figure 4. (a) Scaling of slope S (Equation 26) (b) the ratio of cover threshold to threshold of motion b (Equation17), and (c) of the long-term mean cover C with the discharge variability parameter k.Lines for the cover exponent α = −2 (flood-cleaning), crosses for α = 2 (flood-depositing).Black for the width exponent q = 0, blue for q = 1, and red for q = 2.5.Note that channel bed slope is independent of α (Equation26).For many conditions, there are two solutions available, corresponding to the solutions for b smaller (dashed) or larger than one (solid) (see Table1; Appendices B and C).

Figure 5 .
Figure 5. (a and c) Scaling of the ratio of cover threshold to threshold of motion b (Equation17) and (b and d) of the long-term mean cover with the cover scaling exponent α (see Equations 15 and 16), for a (a and b) low variability (k = 3); and (c and d) high variability climate (k = 3).For flood-cleaning streams (α < 0) (cf.Figure2, Table1), separate solutions are shown for the threshold ratio b (Equation17), for b > 1 (dashed line) and b < 1 (solid line).

Figure 6 .
Figure 6.(a) Slope-area relationship predicted by the model (Equation26) and (b) corresponding elevation profiles.The concavity index θ was calculated using Equation29.Relationships are shown for m = 2 and n = 1, and B = q = 0 (black line), B = 0 and q = 1 (blue line), B = 0 and q = 2.5 (red line), and B = 1 (black dashed).The steepness index was calibrated to the conditions at the Liwu (red circle, Table2), with A = 435 km 2 , S = 0.02 (Table2), and elevation 410 m.Hack's law A = k H x H was used to convert from drainage area to channel length, with k H = 1.2 and H = 2.
yield several other potential solutions, indicating a larger flexibility of the channel to deal with different climatic situations.The model presented here uses a sediment-flux-dependent bedrock incision models rooted in current mechanistic understanding of fluvial bedrock incision.As such, it connects reach-scale to landscape-scale approaches in modeling bedrock river dynamics, addressing what Venditti et al. (2020) called a current grand challenge of geomorphology.The channel long profile predicted by the model yields a power-law dependence of channel bed slope on incision rate and drainage area similar to the one obtained from the stream-power incision model (SPIM) (Equation

Figure 7 .
Figure 7.Comparison of the dependence of incision rate on discharge variability k obtained from the stream power incision model (SPIM, solid line; see Appendix E) and the sediment-flux-dependent model used herein (see Appendix C).Calculations were made for values of the width exponent of (a) q = 0 and (b) q = 5/2 (see Appendix A), and a threshold of motion   1 ct Q.
23 gives a further constraint on erodibility k e in case the long-term incision rate is known.
Discharge variability-dependent function for the long-term mean bed cover.F I Discharge variability-dependent function for the long-term mean incision rate.F Qs Discharge variability-dependent function for the long-term mean sediment transport rate.pdf(Q*) Probability density function of the dimensionless discharge Q* (Equation 1).Γ(x) Gamma function of x (Equation 2).Γ(x,c) Upper incomplete gamma function of x (Equation 14).Variables A Drainage area [m 2 ].A c Cross-sectional area of the flow [m 2 ]. a Scaling exponent, C-Q*.B Scaling exponent, β-A.b Coefficient of proportionality, Downstream hydraulic geometry exponent for width, scaling exponent <W> − Q ω a At-a-station hydraulic geometry exponent for width, scaling exponent W − Q*

Table 2
Parameter Values Used for the Example Calculations, Following Estimates for the Liwu River, at Lushui, Taiwan

Table 2
), and elevation 410 m.Hack's law A = k H x H was used to convert from drainage area to channel length, with k H = 1.2 and H = 2.