Mass balance, grade, and adjustment timescales in bedrock channels

Rivers are dynamical systems that are thought to evolve towards a steady-state configuration. Then, geomorphic parameters, such as channel width and slope, are constant over time. In the mathematical description of the system, the steady state corresponds to a fixed point in the dynamic equations in which all time derivatives are equal to zero. In alluvial rivers, steady state is characterized by grade. This can be expressed as a so-called order principle: an alluvial river evolves to achieve a state in which sediment transport is constant along the river channel and is equal to transport capacity everywhere. In bedrock rivers, steady state is thought to be achieved with a balance between channel incision and uplift. The corresponding order principle is the following: a bedrock river evolves to achieve a vertical bedrock incision rate that is equal to the uplift rate or base-level lowering rate. In the present work, considerations of process physics and of the mass balance of a bedrock channel are used to argue that bedrock rivers evolve to achieve both grade and a balance between channel incision and uplift. As such, bedrock channels are governed by two order principles. As a consequence, the recognition of a steady state with respect to one of them does not necessarily imply an overall steady state. For further discussion of the bedrock channel evolution towards a steady state, expressions for adjustment timescales are sought. For this, a mechanistic model for lateral erosion of bedrock channels is developed, which allows one to obtain analytical solutions for the adjustment timescales for the morphological variables of channel width, channel bed slope, and alluvial bed cover. The adjustment timescale to achieve steady cover is of the order of minutes to days, while the adjustment timescales for width and slope are of the order of thousands of years. Thus, cover is adjusted quickly in response to a change in boundary conditions to achieve a graded state. The resulting change in vertical and lateral incision rates triggers a slow adjustment of width and slope, which in turn affects bed cover. As a result of these feedbacks, it can be expected that a bedrock channel is close to a graded state most of the time, even when it is transiently adjusting its bedrock channel morphology.

1 Introduction

Bedrock rivers are important geomorphic landforms in mountain regions. They set the base level for hillslope response and evacuate the products of erosion, weathering, and hillslope mass wasting (e.g. Hovius and Stark, 2006). As such, they integrate the upstream erosional signal of the landscape, and the material transported in rivers can be used to estimate catchment-averaged denudation rates on various timescales (e.g. Turowski and Cook, 2017). Further, their morphology is thought to be indicative of past climate and tectonic conditions (e.g. Stark et al., 2010; Wobus et al., 2006). Consequently, they provide archives that can be exploited to unravel the Earth's history.

River channels are dynamical systems. Their state variables – for example, slope, cross-sectional shape, and bed roughness – evolve over time under the influence of externally imposed driving variables including water discharge, sediment supply, and tectonic uplift (e.g. Heimann et al., 2015; Lague, 2010; Parker, 1979; Wickert and Schildgen, 2019). Like in many other dynamical systems, there exists a fixed point in the descriptions of river dynamics at which all state variables are constant over time. In an alluvial river, at this fixed point, entrainment and deposition of sediment are in balance along the river profile, implying that sediment transport rate is constant and that sediment transport capacity matches sediment supply. A river that exhibits these features is said to be “in grade” or “graded”, because it is neither aggrading nor degrading (Mackin, 1948). Since its introduction, the graded stream concept has become a central paradigm in river morphodynamics (e.g. Blom et al., 2017; Church, 2006). There are several reasons for this importance. Chiefly, rivers are physically complicated systems, and the description of their steady-state forms is a problem that is considerably simpler than the full description of their dynamics. Further, many variables of natural rivers are challenging to measure. Yet, comparatively simple scaling relations have been observed between variables such as discharge or drainage area, on the one hand, and channel width or channel slope, on the other hand (e.g. Gleason, 2015; Leopold and Maddock Jr., 1953; Whitbread et al., 2015). These scaling relationships are thought to be explainable using steady-state models (e.g. Eaton and Church, 2004; Smith, 1974; Turowski, 2018; Wobus et al., 2006).

The condition of grade in a stream is tightly connected to the description of its sediment mass balance. For alluvial rivers, this mass balance is typically described by one of two approaches: the Exner equation or the entrainment–deposition framework (e.g. An et al., 2018). In the Exner equation (e.g. Chen et al., 2014; Paola and Voller, 2005), the rate of change of the sediment bed elevation h_s is related to the along-stream divergence of sediment supply per unit width, q_s.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial h_{\mathrm{s}}}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{s}}(\mathrm{1}-p)}}{\displaystyle \frac{\partial q_{\mathrm{s}}}{\partial x}}.\end{array}$$

Here, p is the porosity and ρ_s is the density of the sediment, t is the time, and x is the distance in the downstream direction. In steady state, for a graded stream, the time derivative on the left-hand side is zero, which implies that the spatial derivative on the right-hand side is zero also. As a result, the sediment flux is constant along the stream – the stream is in grade. Any bed elevation change leads to an adjustment of slope. The condition of grade thus implies that transport capacity is also constant and equal to sediment supply along the stream. In the entrainment–deposition framework (e.g. Charru et al., 2004; Davy and Lague, 2009; Shobe et al., 2017), the entrainment rate E and deposition rate D of sediment mass per unit area are tracked explicitly, giving the mass balance for the mobile sediment mass per unit area M_m (e.g. Turowski and Hodge, 2017):

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial M_{\mathrm{m}}}{\partial t}}=-{\displaystyle \frac{\partial q_{\mathrm{s}}}{\partial x}}+E-D.\end{array}$$

The sediment bed elevation change is then described by a second equation:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial h_{\mathrm{s}}}{\partial t}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{s}}(\mathrm{1}-p)}}(D-E).\end{array}$$

Within this framework, in steady state, time derivatives are set to zero, implying that entrainment needs to equal deposition (Eq. 3) and sediment flux along the stream needs to be constant (Eq. 2). Again, this means that the stream is in grade. The main advantage of the erosion-deposition framework is that it separately keeps track of stationary and moving sediment mass. This allows one to predict a lagged response of bed elevation to changes in sediment supply, due to the interplay of entrainment, deposition, and lateral sediment movement (e.g. An et al., 2018). Its main disadvantage is that both entrainment and deposition (E and D in Eqs. 2 and 3) need to be quantified in terms of hydraulic drivers. In contrast, to use the Exner equation, only transport capacity or transport rate needs to be quantified, which is considerably easier to measure than deposition and entrainment rate, and therefore the relevant relationships are better constrained. Nevertheless, both approaches are related and the entrainment–deposition equations (Eqs. 2 and 3) can be transformed into the Exner equation (Eq. 1) when combining mobile and stationary mass into a single total mass term (Appendix A).

In bedrock channels, the concept of grade has not been widely applied. One of the main reasons for this is that bedrock channels are usually viewed as detachment-limited systems, where sediment supply is much smaller than transport capacity (e.g. Tinkler and Wohl, 1998; Whipple et al., 2013), which is in direct contrast to the assumption of grade. As a result, the system is assumed to be driven by its potency for erosion (e.g. Whipple, 2004). The evolution of bedrock channel bed elevation h_b is described by the equation (e.g. Howard, 1994)

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial h_{\mathrm{b}}}{\partial t}}=T_{\mathrm{U}}-I.\end{array}$$

Here, T_U is the uplift rate or relative base-level fall rate, and I the bedrock incision rate. According to Eq. (4), bedrock channels adjust to a steady state in which incision rate I equals uplift rate T_U.

Over the last 2 decades, evidence has been mounting that fluvial bedrock erosion is driven by the impacts of sediment particles in many settings (e.g. Cook et al., 2013; Johnson et al., 2010; Sklar and Dietrich, 2004). The amount of sediment in the channel affects erosion rates by two main effects. First, an increase in the number of moving particles leads to an increase in the number of impacts on the bed, increasing erosion rates. This is known as the tools effect. Second, sediment residing on the bed may protect the rock surface from impacts, reducing erosion rates. This is known as the cover effect. Evidence for both tools and cover effects has been described in laboratory and field studies (e.g. Beer et al., 2016; Cook et al., 2013; Johnson and Whipple, 2010; Sklar and Dietrich, 2001; Turowski et al., 2008a). In addition, large sediment bodies are common in many mountain regions, and they can reside in and around stream channels for potentially a long time (e.g. Korup et al., 2006; Schoch et al., 2018). All of these observations imply that a description of the mass balance of sediment should be an essential part of any theoretical description of bedrock channels. In addition, recent observations have been interpreted such that bedrock channels are in a graded state, similarly to alluvial channels (Phillips and Jerolmack, 2016). Thus, it seems that the view that bedrock channels are in a detachment-limited state, in which long-term sediment supply is smaller than transport capacity (e.g. Whipple et al., 2013), is insufficient to account for all observations made in natural streams.

In this paper, I have three separate, yet related, aims. First, I develop a description of the mass balance of bedrock channels, based on previous work by Turowski and Hodge (2017) and Turowski (2018). The mass balance is used to derive and discuss the concept of the graded stream for bedrock channels. Second, I derive expressions for response timescales for bedrock channels to adjust to a graded state. Third, for this it is necessary to develop a description of bedrock channel wall erosion by impacting particles. The concepts are used to discuss the current notion of bedrock channels, their possible routes to a graded state, and the relevant response timescales.

2 Theoretical considerations

2.1 Mass balance equations for sediment

Landscapes form by the interplay of bedrock erosion and the entrainment, transport and deposition of sediment, as determined by various drivers such as climate, tectonics, and biological activity. Each erosion process has a minimum of two phases: the breakdown of rock mass by chemical or physical weathering and the entrainment and evacuation of loose pieces of rock that are produced in this way (Gilbert, 1877). From this, it is clear that a minimum description of any eroding landscape needs to include a mass balance equation each for bedrock and for loose sediment. Consider a control volume within a river (Fig. 1), with width W, length L, and a height ranging from the surface, i.e. the interface between bedrock or sediment and the water, to a fixed reference level somewhere in the bedrock below. The loose material, sediment, overlays the bedrock. Uplift pushes new bedrock into the control volume at a rate T_U, while incision converts it into sediment at a rate I. We assume that the erosion products are small enough so that they are subsequently transported in suspension. Then the rate of change of bedrock mass per unit area M_b is given by

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial M_{\mathrm{b}}}{\partial t}}={\mathit{\rho }}_{\mathrm{r}}\left(T_{\mathrm{U}}-I\right).\end{array}$$

Here, ρ_r represents the density of the bedrock. Dividing Eq. (5) by ρ_r, and realizing that $h_{\mathrm{b}}=M_{\mathrm{b}}/{\mathit{\rho }}_{\mathrm{r}}$, we retrieve the usual form of the bedrock mass balance (Eq., 4). Details of the derivation of the mass balance for sediment have been given by Turowski and Hodge (2017). Note that working with mass instead of a deposit thickness is advantageous for bedrock channels, because sediment may not be equally distributed on the bed. The entrainment–deposition framework is preferable, because it makes it possible to distinguish between moving and stationary sediment, which is necessary to treat the cover and the tools effects. This is not possible when using the Exner approach (Appendix A). The mass balance for the mobile sediment per unit area M_m is given by Eq. (2):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial M_{\mathrm{m}}}{\partial t}}=-{\displaystyle \frac{\partial q_{\mathrm{s}}}{\partial x}}+E-D.\end{array}$$

The mass balance for the stationary sediment per unit area M_s is given by

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial M_{\mathrm{s}}}{\partial t}}=D-E.\end{array}$$

Finally, sediment flux per unit width q_s and the mobile sediment mass per unit area M_m are connected via the downstream particle speed U:

$$\begin{array}{}\text{(8)}& q_{\mathrm{s}}=U{M}_{\mathrm{m}}.\end{array}$$

Figure 1 Schematic side view of a control volume within a bedrock channel. The bedrock (bottom) is overlain by stationary sediment (centre), which exchanges particles via entrainment E and deposition D with the mobile sediment in the water column (top). The bedrock surface h_b lowers at the incision rate I, while the sediment surface h_s evolves according to the balance of entrainment and deposition (Eq. 6).

2.2 Lateral erosion in bedrock channels by impacting particles

Considering impact erosion to be the dominant erosion process, the lateral erosion rate E_L of bedrock channels is driven by particle impacts. It can therefore, similarly to the formulation of the saltation–abrasion model (Sklar and Dietrich, 2004), be written as the product of two terms: (i) the average volume eroded by a single impact, V_i, and (ii) the impact rate per area and time I_r. The latter term can be subdivided into two terms. The first of these quantifies the number of available particles per unit time and area, F_T, which describes the tools effect. The second term, F_C, describes the effect of bed cover, which captures the effects of the distribution of sediment in the channel on lateral erosion. The need for this term arises because bedload particles generally travel parallel to the channel walls. Sideward deflection is controlled by the interaction of moving particles with the bed (Beer et al., 2017; Fuller et al., 2016) and specifically with stationary sediment, i.e. bed cover (cf. Turowski, 2018).

$$\begin{array}{}\text{(9)}& E_{\mathrm{L}}=V_{\mathrm{i}}{I}_{\mathrm{r}}=V_{\mathrm{i}}{F}_{\mathrm{T}}{F}_{\mathrm{C}}\end{array}$$

The volume eroded per impact for lateral erosion should be the same as for vertical erosion and has been quantified by Sklar and Dietrich (2004) as the energy of the impact divided by a material constant. It can be evaluated by

$$\begin{array}{}\text{(10)}& V_{\mathrm{i}}={\displaystyle \frac{\mathrm{2}Y}{k_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}{\displaystyle \frac{M_{\mathrm{p}}{w}_{\mathrm{i}}^{\mathrm{2}}}{\mathrm{2}}}.\end{array}$$

Here, the first term is related to material properties, where Y and σ_T are Young's modulus of the bedrock and its tensile strength, respectively, and k_v is the rock resistance coefficient. The second term gives the kinetic energy of the impacting grain. Here, M_p is the mass of a single particle and w_i is the impact speed normal to the wall.

As in vertical bedrock erosion (Beer and Turowski, 2015; Inoue et al., 2014; Sklar and Dietrich, 2004), the tools effect can be modelled as a linear function of bedload supply Q_s (Mishra et al., 2018) multiplied by a dimensionless factor κ_T with values between 0 and 1 that describes the fraction of bedload available for lateral erosion. To obtain the number of impacting particles per unit area, this product needs to be divided by the mass of a single particle and the total area of the wall, A_w, that is eroded to give

$$\begin{array}{}\text{(11)}& F_{\mathrm{T}}={\displaystyle \frac{{\mathit{\kappa }}_{\mathrm{T}}{Q}_{\mathrm{s}}}{A_{\mathrm{w}}{M}_{\mathrm{p}}}}.\end{array}$$

Substituting Eqs. (10) and (11) into Eq. (9), the lateral erosion rate of a bedrock channel can thus be written as

$$\begin{array}{}\text{(12)}& E_{\mathrm{L}}={\mathit{\kappa }}_{\mathrm{T}}{\displaystyle \frac{Y}{k_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}{\displaystyle \frac{Q_{\mathrm{s}}{w}_{\mathrm{i}}^{\mathrm{2}}}{A_{\mathrm{w}}}}{F}_{\mathrm{CD}}.\end{array}$$

In Eq. (12), there are three parameters that require further discussion: the impact speed w_i, the eroded area A_w, and the cover-dependent term F_CD. In a previous paper (Turowski, 2018), I argued that lateral erosion and channel width development are intimately related to bed cover. The quantification of all three parameters springs from the physical–conceptual model developed in this previous paper. For this reason the cover-dependent term, F_CD, will be discussed first, leading to a quantification of the other two terms, w_i and A_w.

In a straight bedrock channel the motion of water and sediment is generally parallel to the walls. Lateral erosion occurs when sediment particles are deflected sideways such that they impact the walls with sufficient force to cause damage. For a given reach, we can define a sideward deflection length scale d, which is relevant for reach-scale lateral erosion (Turowski, 2018). The relevant cross section for setting reach-scale channel width is assumed to be located where the sinuous bedload particle stream crosses from the gravel bar onto the smooth bedrock at the apex of the bar (Fig. 2). Only there, several conditions come together that are favourable to achieve the maximal sideward deflection distances (Turowski, 2018). These are (i) the high particle concentration, (ii) a vector of motion of the particle stream that is already pointing towards the walls, (iii) the existence of roughness necessary for sideward deflection provided by the alluvium, and (iv) the smooth bedrock that does not hinder sideward motion. We expect that the wall is eroded if the uncovered width W_uncovered in the cross section is smaller than d (Fig. 3). As a result, we can quantify the cover-dependent term F_CD as

$$\begin{array}{}\text{(13)}& F_{\mathrm{CD}}=\left\{\begin{array}{ll}\mathrm{1}& \mathrm{if}d>W_{\mathrm{uncovered}}\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right..\end{array}$$

We can write the eroded area on the wall A_w as the product of a length scale and a height. From the argument above, the same particle can attack the wall once when passing each gravel bar. Therefore, the relevant length scale for lateral erosion is the distance between bars on a given side of the channel, i.e. the wavelength of bar spacing, λ. Fuller et al. (2016) observed that for sideward-deflected particles, the erosion height on the wall is larger than the typical saltation hop height. Beer et al. (2017) observed a similar increase in wall erosion rates near boulder obstacles in the channel. When the roughness elements that cause deflection are related to stationary alluvium, we can expect that the height scale is the maximum saltation height of bedload particles at the wall, H_w.

$$\begin{array}{}\text{(14)}& A_{\mathrm{w}}=\mathit{\lambda }{H}_{\mathrm{w}}\end{array}$$

Note that the entire area is not eroded at the same time. Rather, particles are deflected towards the wall near the apex of the bars (cf. Turowski, 2018). Consequently, only a small area is eroded at a given time, and the locus of erosion slowly moves downstream as the bars migrate. Likewise, the impact speed w_i and the sideward deflection distance d are related to saltation properties. It is, of course, possible that a particle undergoes several saltation cycles until it impacts the wall. However, in this case, in each additional saltation hop, the sideward component of motion would reduce due to downstream hydraulic forces and frictional loss of momentum. Here, I assume that only during the first hop do particles have sufficient lateral momentum to cause erosion upon impact on the wall. This assumption needs to be verified experimentally.

Figure 2 Schematic top view of a straight bedrock channel, with alternating submerged gravel bars (dark grey) on a bedrock bed (white). The sinuous thalweg (light grey) and bedload path (transparent dark grey) are indicated. The dashed black line indicates the cross section that is ideal for sideward deflection of particles; here, the bedload particle stream crosses the boundary between gravel and smooth bedload. The wavelength of the alternating bars and therefore of the bedload path should scale with channel width. Adapted from Turowski (2018).

Figure 3 The sideward deflection length scale d interacts with bed cover and channel width to determine whether the lateral erosion occurs (a) or not (b, c). Adapted from Turowski (2018).

Since, within the model, sideward deflection is caused by stationary alluvium, particle trajectories should follow those observed for saltation over alluvium (e.g. Abbot and Francis, 1977; Niño et al., 1994), rather than those over bedrock (e.g. Chatanantavet et al., 2013; Auel et al., 2017a). Because the wall-normal component of the motion is relevant for impact erosion (e.g. Sklar and Dietrich, 2004), the particle trajectory needs to be corrected for its angle γ of motion with respect to the wall. Then, the sideward deflection distance d is related to the saltation hop length L_s by

$$\begin{array}{}\text{(15)}& d=L_{\mathrm{s}}\sin \left(\mathit{\gamma }\right).\end{array}$$

Likewise, the impact speed w_i is related to the particle speed U by

$$\begin{array}{}\text{(16)}& w_{\mathrm{i}}=U\sin \left(\mathit{\gamma }\right).\end{array}$$

Auel et al. (2017a) proposed empirical equations to describe saltation properties over a sediment bed as a function of hydraulics, based on their own experiments and a data compilation from various sources. They give the saltation hop length L_s by

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{L_{\mathrm{s}}}{D}}=\mathrm{1.17}\left({\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}}-\mathrm{1}\right).\end{array}$$

Here, D is grain diameter, θ is the Shields stress, and θ_c is the critical Shields stress for the onset of bedload motion. Similarly, hop height H_s is given by

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{H_{\mathrm{s}}}{D}}=\mathrm{0.025}\left[\left({\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}}-\mathrm{1}\right)+\mathrm{24}\right]\end{array}$$

and downstream particle speed U is given by

$$\begin{array}{}\text{(19)}& U=\mathrm{1.46}{\left[\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right)gD\right]}^{\mathrm{0.5}}{\left({\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}}-\mathrm{1}\right)}^{\mathrm{0.5}}.\end{array}$$

Here, ρ and ρ_s are the densities of the water and sediment, respectively, and g is the acceleration due to gravity. Finally, to close the system of equations, we need some relations describing the geometry of the gravel bars. Alternating bars in bedrock channels have been little studied (e.g. Nelson and Seminara, 2012), and the necessary relations are not available. From a large data compilation of bar width and length in braided channels, Kelly (2006) found that bar length L_bar is related to bar width W_bar by

$$\begin{array}{}\text{(20)}& L_{\mathrm{bar}}=\mathrm{4.95}{W}_{\mathrm{bar}}^{\mathrm{0.97}}.\end{array}$$

Based on this observation, I assume that in bedrock channels the wavelength of the bars scales with their width, such that

$$\begin{array}{}\text{(21)}& \mathit{\lambda }={\mathit{\kappa }}_{\mathrm{bar}}{W}_{\mathrm{bar}}={\mathit{\kappa }}_{\mathrm{bar}}{W}_{\mathrm{covered}}.\end{array}$$

Here, the bar width has been identified with the covered width W_covered (Figs. 2 and 3), and κ_bar is a dimensionless constant with a value of 2–10, by analogy with bar shapes in alluvial rivers (e.g. Kelly, 2006). Bed cover C is the ratio of covered bed area A_c to total bed area A_tot, which can be related to the covered width W_covered as follows:

$$\begin{array}{}\text{(22)}& C={\displaystyle \frac{A_{\mathrm{c}}}{A_{\mathrm{tot}}}}={\displaystyle \frac{W_{\mathrm{covered}}}{W}}={\displaystyle \frac{W_{\mathrm{bar}}}{W}}.\end{array}$$

Here, W is the channel width. As a result, the bar length can be written as

$$\begin{array}{}\text{(23)}& \mathit{\lambda }={\mathit{\kappa }}_{\mathrm{bar}}WC.\end{array}$$

Assuming that the maximum saltation height at the wall corresponds to the maximum saltation hop height, H_w=H_s, and substituting Eqs. (13) to (16) and Eqs. (18) to (23) into Eq. (12), we obtain

$$\begin{array}{}\text{(24)}& \begin{array}{rl}& E_{\mathrm{L}}=\\ & \left\{\begin{array}{ll}\frac{\mathit{\kappa }Yg}{k_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}\left(\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}-\mathrm{1}\right){\sin }^{\mathrm{2}}\left(\mathit{\gamma }\right)\frac{Q_{\mathrm{s}}}{WC}\frac{\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}-\mathrm{1}\right)}{\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}-\mathrm{1}\right)+\mathrm{24}}& \mathrm{if}d>W_{\mathrm{uncovered}}\mathrm{and}\mathit{\theta }>{\mathit{\theta }}_{\mathrm{c}}\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right..\end{array}\end{array}$$

Here, $\mathit{\kappa }=\mathrm{85}{\mathit{\kappa }}_{\mathrm{T}}/{\mathit{\kappa }}_{\mathrm{bar}}$ is a dimensionless constant. The sideward deflection length scale d can be estimated by the hop length L_s (Eq. 24):

$$\begin{array}{}\text{(25)}& d=\mathrm{1.17}D\sin \left(\mathit{\gamma }\right)\left({\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}}-\mathrm{1}\right).\end{array}$$

Finally, the uncovered width can be related to bed cover using Eq. (22).

$$\begin{array}{}\text{(26)}& W_{\mathrm{uncovered}}=W-W_{\mathrm{covered}}=W(\mathrm{1}-C)\end{array}$$

The rate of change of channel width, in the case of a widening channel, should be twice the lateral erosion rate given in Eq. (24), since both sides are eroded at the same time.

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\mathrm{d}W}{\mathrm{d}t}}=\mathrm{2}{E}_{\mathrm{L}}\end{array}$$

Note that, when d=W_uncovered, the model gives a steady-state channel width consistent with the model of Turowski (2018), with the sideward deflection distance given by Eq. (25).

2.3 Timescales of morphological adjustment in bedrock channels

I will now derive analytical expressions for the response time of the channel to perturbations in the boundary conditions, such as changes in discharge, sediment supply, or uplift rate. This will be done for three key parameters: channel bed slope, channel width, and cover. For the derivation, it is necessary to assume that, on the timescale of adjustment of one variable, the other variables stay essentially constant. This assumption is reasonable if a particular variable adjusts much slower than another variable. For example, slope adjustment takes much longer time than the adjustment of bed cover.

2.3.1 Response time of channel bed slope

Taking the spatial derivative of Eq. (4) and assuming spatially constant uplift rate T_U, we obtain

$$\begin{array}{}\text{(28)}& {\displaystyle \frac{\partial }{\partial x}}{\displaystyle \frac{\partial h_{\mathrm{b}}}{\partial t}}=-{\displaystyle \frac{\partial I}{\partial x}}.\end{array}$$

Channel bed slope S is defined as the topographic gradient in the downstream direction:

$$\begin{array}{}\text{(29)}& S=-{\displaystyle \frac{\partial h_{\mathrm{b}}}{\partial x}}.\end{array}$$

Equation (28) can thus be rewritten as

$$\begin{array}{}\text{(30)}& {\displaystyle \frac{\partial I}{\partial x}}=-{\displaystyle \frac{\partial }{\partial x}}{\displaystyle \frac{\partial h_{\mathrm{b}}}{\partial t}}=-{\displaystyle \frac{\partial }{\partial t}}{\displaystyle \frac{\partial h_{\mathrm{b}}}{\partial x}}={\displaystyle \frac{\partial S}{\partial t}}.\end{array}$$

According to the revised saltation–abrasion equation by Auel et al. (2017b), the vertical erosion equation takes the form

$$\begin{array}{}\text{(31)}& I={\displaystyle \frac{gY}{\mathrm{230}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right){\displaystyle \frac{Q_{\mathrm{s}}}{W}}(\mathrm{1}-C).\end{array}$$

Steady-state cover can be described with the equation by Turowski and Hodge (2017):

$$\begin{array}{}\text{(32)}& C=\left(\mathrm{1}-e^{-\frac{Q_{\mathrm{s}}}{M_{\mathrm{0}}UW}}\right){\displaystyle \frac{Q_{\mathrm{s}}}{Q_{\mathrm{t}}}}.\end{array}$$

The bedload transport capacity can be written as

$$\begin{array}{}\text{(33)}& {\displaystyle \frac{Q_{\mathrm{t}}}{W}}=K_{\mathrm{bl}}{Q}^{\mathrm{m}}{S}^{\mathrm{n}}.\end{array}$$

Substituting Eqs. (29) to (33) into Eq. (28), and assuming that all variables apart from slope are constant, the slope evolution equation takes the form

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{\partial S}{\partial t}}+nB{S}^{-n-\mathrm{1}}{\displaystyle \frac{\partial S}{\partial x}}=\mathrm{0}.\end{array}$$

Here, B is assumed to be constant.

$$\begin{array}{}\text{(35)}& B={\displaystyle \frac{gY}{\mathrm{230}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right)\left(\mathrm{1}-e^{-\frac{Q_{\mathrm{s}}}{M_{\mathrm{0}}WU}}\right){\displaystyle \frac{Q_{\mathrm{s}}^{\mathrm{2}}}{K_{\mathrm{bl}}{W}^{\mathrm{2}}{Q}^{\mathrm{m}}}}.\end{array}$$

Equation (34) is a non-linear wave equation with celerity c_S

$$\begin{array}{}\text{(36)}& c_{\mathrm{S}}=nB{S}^{-n-\mathrm{1}}.\end{array}$$

The timescale of slope adjustment T_S can therefore be written as

$$\begin{array}{}\text{(37)}& \begin{array}{rl}{T}_{\mathrm{S}}& ={\displaystyle \frac{L}{c_{\mathrm{S}}}}={\displaystyle \frac{L{S}^{n+\mathrm{1}}}{nB}}={\displaystyle \frac{k_{\mathrm{bl}}{Q}^{\mathrm{m}}L{S}^{n+\mathrm{1}}}{nk\left(\mathrm{1}-e^{-\frac{\stackrel{\mathrm{‾}}{q_{\mathrm{s}}}}{M_{\mathrm{0}}U}}\right)q_{\mathrm{s}}^{\mathrm{2}}}}\\ & ={\displaystyle \frac{q_{\mathrm{t}}LS}{nk\left(\mathrm{1}-e^{-\frac{\stackrel{\mathrm{‾}}{q_{\mathrm{s}}}}{M_{\mathrm{0}}U}}\right)q_{\mathrm{s}}^{\mathrm{2}}}}={\displaystyle \frac{LSW}{nk{Q}_{\mathrm{s}}C}}.\end{array}\end{array}$$

Here, L is the length of the reach in question, and k is the erodibility, which, according to the revised saltation–abrasion equation by Auel et al. (2017b), takes the form

$$\begin{array}{}\text{(38)}& k={\displaystyle \frac{gY}{\mathrm{230}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right).\end{array}$$

2.3.2 Response time of channel width

For the adjustment of channel width, it is necessary to distinguish between narrowing and widening channels. While channel widening is controlled by the lateral erosion of bedrock walls (see Sect. 2.2, Eq. 24), a bedrock channel can only narrow when incising vertically. Therefore, the response timescale of narrowing is related to the vertical incision rate. The timescale of narrowing can be estimated by the time necessary to incise the flow depth H. After this time, the wetted channel cross section has been completely replaced. Thus, using the continuity equation (Eq. D3) and the expression for flow velocity (Eq. D4), the timescale of channel narrowing is

$$\begin{array}{}\text{(39)}& T_{\mathrm{N}}={\displaystyle \frac{H}{I}}={\displaystyle \frac{(gS{)}^{\frac{\mathit{\alpha }-\mathrm{1}}{\mathrm{2}}}{R}^{\frac{\mathrm{3}\mathit{\alpha }-\mathrm{1}}{\mathrm{2}}}}{k_{\mathrm{V}}I}}{\left({\displaystyle \frac{Q}{W}}\right)}^{\mathrm{1}-\mathit{\alpha }}.\end{array}$$

The technique of perturbation analysis can be used to obtain an analytical solution for the width response time in the case of a widening channel (e.g. Braun et al., 2015; Turowski and Hodge, 2017). The mathematical details are given in Appendix C, leading to the equation

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{T}_{\mathrm{W}}& ={\displaystyle \frac{\mathrm{18}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}{\mathit{\kappa }Y\left(\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}-\mathrm{1}\right)g}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{c}}}{\mathit{\theta }}}{\displaystyle \frac{W^{\mathrm{2}}}{Q_{\mathrm{t}}}}\left({\displaystyle \frac{\mathrm{3}}{\mathrm{2}C}}{\displaystyle \frac{Q_{\mathrm{s}}}{M_{\mathrm{0}}UW}}\left({\displaystyle \frac{\mathrm{1}}{C}}{\displaystyle \frac{Q_{\mathrm{s}}}{Q_{\mathrm{t}}}}-\mathrm{1}\right)\right.\\ & {\left.-{\displaystyle \frac{(\mathit{\alpha }-\mathrm{1})}{C}}+(\mathit{\alpha }-\mathrm{2})\right)}^{-\mathrm{1}}.\end{array}\end{array}$$

Here, M₀ is the minimum mass necessary to cover the bed per unit area, and α≈0.6 is a dimensionless exponent that appears in the flow velocity equation (see Eq. D4; Nitsche et al., 2012). The minimum mass M₀ can be evaluated by assuming that a single layer of close-packed spherical grains resides on the bed (Turowski, 2009; Turowski and Hodge, 2017)

$$\begin{array}{}\text{(41)}& M_{\mathrm{0}}={\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\mathrm{s}}D}{\mathrm{3}\sqrt{\mathrm{3}}}}.\end{array}$$

2.3.3 Response time of bed cover

The response time for the adjustment of bed cover T_C was previously derived by Turowski and Hodge (2017) and is given by

$$\begin{array}{}\text{(42)}& T_{\mathrm{C}}={\displaystyle \frac{L{M}_{\mathrm{0}}W}{Q_{\mathrm{t}}C}}.\end{array}$$

2.3.4 Response time ratios

The dynamics of the channel during adjustment is to some extent determined by the relative magnitude of the response times. For example, if the response time for the adjustment of bed slope is always much longer than the response time for bed cover, on the timescale of slope adjustment, it can be assumed that bed cover is always at a steady state. The ratio of the response time for slope and width (widening channel) is given by

$$\begin{array}{}\text{(43)}& \begin{array}{rl}{\displaystyle \frac{T_{\mathrm{S}}}{T_{\mathrm{W}}}}=& {\displaystyle \frac{\mathrm{115}\mathit{\kappa }}{\mathrm{9}n}}{\displaystyle \frac{SL}{WC}}{\displaystyle \frac{Q_{\mathrm{t}}}{Q_{\mathrm{s}}}}{\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}}\left({\displaystyle \frac{(\mathit{\alpha }-\mathrm{1})}{C}}-{\displaystyle \frac{\mathrm{3}}{\mathrm{2}C}}{\displaystyle \frac{Q_{\mathrm{s}}}{M_{\mathrm{0}}UW}}\right.\\ & \left.\left({\displaystyle \frac{\mathrm{1}}{C}}{\displaystyle \frac{Q_{\mathrm{s}}}{Q_{\mathrm{t}}}}-\mathrm{1}\right)-(\mathit{\alpha }-\mathrm{2})\right).\end{array}\end{array}$$

Similarly, for a narrowing channel

$$\begin{array}{}\text{(44)}& {\displaystyle \frac{T_{\mathrm{S}}}{T_{\mathrm{N}}}}={\displaystyle \frac{k_{\mathrm{V}}ILSW}{nk{Q}_{\mathrm{s}}C}}(gS{)}^{\frac{\mathrm{1}-\mathit{\alpha }}{\mathrm{2}}}{R}^{\frac{\mathrm{1}-\mathrm{3}\mathit{\alpha }}{\mathrm{2}}}{\left({\displaystyle \frac{Q}{W}}\right)}^{\mathit{\alpha }-\mathrm{1}}.\end{array}$$

The ratio of the response time for cover and slope is given by

$$\begin{array}{}\text{(45)}& {\displaystyle \frac{T_{\mathrm{C}}}{T_{\mathrm{S}}}}={\displaystyle \frac{gY}{\mathrm{230}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right){\displaystyle \frac{n{M}_{\mathrm{0}}}{S}}{\displaystyle \frac{Q_{\mathrm{s}}}{Q_{\mathrm{t}}}}.\end{array}$$

The ratio of the response time for cover and width is given by

$$\begin{array}{}\text{(46)}& \begin{array}{rl}{\displaystyle \frac{T_{\mathrm{C}}}{T_{\mathrm{W}}}}& ={\displaystyle \frac{\mathit{\kappa }Y\left(\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}-\mathrm{1}\right)g{M}_{\mathrm{0}}}{\mathrm{18}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}{\displaystyle \frac{L}{WC}}{\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}}\left({\displaystyle \frac{(\mathit{\alpha }-\mathrm{1})}{C}}\right.\\ & \left.-{\displaystyle \frac{\mathrm{3}}{\mathrm{2}C}}{\displaystyle \frac{Q_{\mathrm{s}}}{M_{\mathrm{0}}UW}}\left({\displaystyle \frac{\mathrm{1}}{C}}{\displaystyle \frac{Q_{\mathrm{s}}}{Q_{\mathrm{t}}}}-\mathrm{1}\right)-(\mathit{\alpha }-\mathrm{2})\right).\end{array}\end{array}$$

Similarly, for a narrowing channel

$$\begin{array}{}\text{(47)}& {\displaystyle \frac{T_{\mathrm{S}}}{T_{\mathrm{N}}}}={\displaystyle \frac{k_{\mathrm{V}}IL{M}_{\mathrm{0}}W}{Q_{\mathrm{t}}C}}(gS{)}^{\frac{\mathrm{1}-\mathit{\alpha }}{\mathrm{2}}}{R}^{\frac{\mathrm{1}-\mathrm{3}\mathit{\alpha }}{\mathrm{2}}}{\left({\displaystyle \frac{Q}{W}}\right)}^{\mathit{\alpha }-\mathrm{1}}.\end{array}$$

3 Results

To illustrate the dependence of channel morphology and of the adjustment timescales on control and channel morphology parameters, I used parameter values for Lushui at the Liwu River, Taiwan (Table 1; see Turowski et al., 2007). The values of reach parameters were either measured in the field or estimated using literature data. The value for discharge is representative of bedload-carrying flows, using the partitioning method proposed by Sklar and Dietrich (2006). The values of the exponent and prefactor of the flow velocity equation (Eq. D4) were selected using data by Nitsche et al. (2012).

Table 1 Parameter values used for the example calculations, following estimates by Turowski et al. (2007) for the Liwu River, at Lushui, Taiwan.

3.1 Steady-state channel morphology

The sideward deflection length scale d is an important parameter setting channel morphology in steady state, in particular the channel width, which depends on the square root of d (Turowski, 2018).

$$\begin{array}{}\text{(48)}& W=\sqrt{{\displaystyle \frac{k{Q}_{\mathrm{s}}d}{I}}}=\sqrt{{\displaystyle \frac{\mathrm{1.17}D\sin \left(\mathit{\gamma }\right)\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}-\mathrm{1}\right)k{Q}_{\mathrm{s}}}{I}}}\end{array}$$

Here, d is estimated using saltation hop length of bedload particles over bare bedrock (Eq. 25). Saltation hop length is dependent on the Shields stress, and the new formulation can consequently alter steady-state scaling of channel width and slope. Unfortunately, Eq. (48) cannot be solved analytically, since Shields stress θ is non-linearly dependent on channel width and slope (see Eq. D6), and a numerical solution is necessary (Fig. 4). As in the model by Turowski (2018), channel width is independent of discharge (Fig. 4a), and the observed scaling between width and discharge arises from a co-dependence of discharge and sediment supply (see Fig. 4b).

Figure 4 Steady-state channel width (solid line), channel bed slope (dashed line), and bed cover (dotted line) against forcing variables discharge (a), sediment supply (b), and uplift rate (c). For the calculations, all other parameters have been kept constant (Table 1).

3.2 Controls on adjustment timescales

For the calculation of adjustment timescales, the dependence of width and slope on discharge, sediment supply, and uplift rate, and on each other, needs to be explicitly taken into account. From the derivation (Appendix C), the relevant width and slope in the timescale equations are those of the steady-state morphology corresponding to the relevant control variables. As such, they are not independent of sediment supply, discharge, or other control variables. Within the model, steady-state channel width and slope cannot be evaluated analytically or written in a closed-form equation. Thus, a numerical solution is necessary. Adjustment timescales of width are generally longer than those for slope and for cover (Fig. 5), at least for the parameter values used in the example calculations (Table 1).

Figure 5 Timescales (a, c, e) and timescale ratios (b, d, f) for channel adjustment, using appropriate steady-state values corresponding to imposed discharge, sediment supply, and uplift rate for slope, width, and cover against forcing variables discharge (a, b), sediment supply (c, d), and uplift rate (e, f). For the calculations, all other parameters have been kept constant (Table 1). For the timescale ratios (b, d, f), only the timescale for widening channels was used, due to its similarity with the timescale for narrowing channels (a, c, e). The solid red line in (b, d, f) indicates a ratio of one.

4 Discussion

4.1 Lateral erosion equation

Equation (24) is a mechanistic description of lateral fluvial bedrock erosion by impacting particles. Field and laboratory data that can be used to test the model are scarce, and the few data sets that exist do not include information on all necessary parameters to test it (e.g. Cook et al., 2014; Fuller et al., 2016; Suzuki, 1982; Mishra et al., 2018). The minimum parameters needed for a meaningful test are the lateral erosion rates measured in parallel with relevant driving variables, including water discharge and bedload transport rate, in a channel with self-formed sediment cover and alternating gravel bars. Nevertheless, the model provides a starting point for future investigations, providing a clear mechanistic description and a host of testable assumptions and predictions.

Due to a lack of direct relevant data and to keep the complexity of the model reasonable, it was necessary to make some simplifications and assumptions on relevant processes and geomorphic response. For example, bedrock channels at high slopes tend to adjust their bed into a step-pool morphology (Duckson and Duckson, 1995; Scheingross et al., 2019). The feedbacks necessary to develop these bedforms, and how they may affect the flow hydraulics and erosion rates, have not been considered in the present model (e.g. Scheingross and Lamb, 2017; Yager et al., 2012). In addition, it was necessary to quantify the wavelength of alternating bars. For the considerations on timescales presented here, the assumption of steady-state cover had to be made, implying fully developed bars and ignoring a potential braiding instability at large channel widths. Nelson and Seminara (2012) provided a linear stability analysis of bar formation over an initially bare bed. They stated explicitly that their considerations do not apply to the geometry of fully formed bars. However, their results and numerical model predictions by Inoue et al. (2016) could be interpreted to suggest that during the transient adjustment to fully formed bars from an initially empty bed under constant forcing conditions, bar wavelength varies little over time. Experimental evidence is rare. Some circumstantial observations can be found in the paper of Chatanantavet and Parker (2008), but these authors do not provide a systematic investigation or conclusive evidence for any type of scaling. In summary, none of the available studies was set up to investigate the controls of fully formed alternating bars, and a full understanding of the controls of their geometry is currently lacking. In absence of a full theory of alternating bars in bedrock channels, I have chosen to keep bar aspect ratio constant (Eq. 21), by analogy with observations in alluvial channels (e.g. Kelly, 2006). Yet, due to the coupling with bed cover (Eq. 23), this assumption leads to unphysical behaviour in the limit of small degrees of cover. In this case, the bar wavelength is small, implying small bar width in comparison to channel width. As a consequence, the meandering bedload path has a large amplitude in comparison to its wavelength, and the deflection angle γ approaches 90^∘. The assumption about bar wavelength is a minor piece in the model, affecting only the response time channel widening, which is linearly dependent on bar aspect ratio. For a full treatment of bar wavelength, we can speculate on the behaviour at two limits. First, at low values of cover, bar wavelength should be independent of cover, and it is likely controlled by channel width or depth. Second, at high cover values, neighbouring bars start to overlap and the relationship to cover likely becomes more complicated. For a fully covered bed, bar dynamics should resemble those in alluvial channels. Further theoretical and experimental investigations are necessary to resolve this issue.

The lateral erosion equation (Eq. 24) generally aligns with expected relations. Lateral erosion rates increase with increasing shear stress, sediment supply, and erodibility. However, they are inversely proportional to bed cover. This negative relationship arises because gravel bars increase their length as cover increases, due to their constant aspect ratio (Eq. 23). This leads to less frequent impacts on the wall by travelling bedload. Fuller et al. (2016) observed that bedrock wall erosion is positively correlated with bed roughness in laboratory experiments. Similarly, Beer et al. (2017) observed higher wall erosion rates next roughness elements in a field study. The data from both of these papers are not sufficient for constraining a functional relationship between roughness and lateral erosion rates. In the model, lateral erosion rate (Eq. 24) depends implicitly on roughness, with a positive relationship, via the dependence on shear stress (see Eq. D6). A similar implicit dependence can be found for the sideward deflection distance d (Eq. 25). Nevertheless, dedicated data on sideward deflection distances are needed to test the current equations and to guide future theoretical developments. Another aspect that is lacking in the current formulation is the dependence of lateral erosion rate on channel curvature. Recent work has attempted to address this within the stream-power framework of bedrock erosion (e.g. Langston and Tucker, 2018; Limaye and Lamb, 2014). Including channel curvature into the present model needs further work on bar deposition and bedload paths within curved channels (cf. Bunte et al., 2006; Fernández et al., 2019; Mishra et al., 2018; Turowski, 2018).

4.2 Steady-state channel morphology

In comparison to the model by Turowski (2018), the sideward deflection length scale d has been explicitly quantified in terms of hydraulics (Eq. 25), which may alter steady-state relationships in comparison to the previously published model. In general, the updated model's predictions align with the results of Turowski (2018). It is somewhat surprising that channel width, like in the model by Turowski (2018), is explicitly independent of discharge (Fig. 4a), and, instead, it is set by sediment supply (Fig. 4b). This implies that channel bed slope adjusts to changes in discharge without an effect on channel width, as long as sediment supply stays constant. The results arise because slope and discharge only feature in the same two equations, in that for Shields stress (Eq. D6) and that for the bedload transport capacity (Eq. C7). Using common parameter values for the relevant exponents m, n, and α, the relationship between slope and discharge is the same in these two equations, allowing the two parameters to co-vary without affecting other parameters. Considering all other parameters constant, the first of these (Eq. D6) gives the relation

$$\begin{array}{}\text{(49)}& S\sim Q^{\frac{\mathrm{2}(\mathit{\alpha }-\mathrm{1})}{\mathit{\alpha }+\mathrm{1}}},\end{array}$$

while the second one (Eq. C7) gives the relation

$$\begin{array}{}\text{(50)}& S\sim Q^{-\frac{m}{n}}.\end{array}$$

With the common parameter choice of α=0.6 (see Nitsche et al., 2012), m=1, and n=2 (see Turowski, 2018), we find that the two exponents are equal:

$$\begin{array}{}\text{(51)}& {\displaystyle \frac{\mathrm{2}(\mathit{\alpha }-\mathrm{1})}{\mathit{\alpha }+\mathrm{1}}}=-{\displaystyle \frac{m}{n}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

Thus, a change in discharge can be offset by a change in slope, without the need to vary any of the other parameters. Mathematically, this means that by substitution the number of parameters and equations can be reduced by one, and slope can be eliminated. A different choice of the m∕n ratio or of α would yield a direct dependence of width on discharge, and a dynamic co-evolution of slope and width.

4.3 Order principles and grade in bedrock channels

The condition of grade can be stated as what I call an order principle, which is a principle after which a dynamic system adjusts state variables to comply with forcing variables. Considering a stream without tributaries or hillslope sediment supply, the order principle for the condition of grade can be stated as follows: a river adjusts such that sediment flux is constant along the stream. The order principle is a direct consequence of the description of the sediment mass balance of the stream (see Sect. 1).

Unlike alluvial channels, which feature a single type of material (the alluvium), in bedrock channels we need to also consider bedrock. This necessitates a second mass balance equation for bedrock (Eq. 4), in addition to that for alluvium (see Sect. 2.1). Accepting that a sediment mass balance cannot be neglected for a mechanistic description of bedrock channel dynamics, a bedrock river thus adjusts to two order principles, rather than one. The first of these is related to the mass balance of sediment (Sect. 2.1) and leads to a condition of grade, as discussed above. The second of these is related to the mass balance of bedrock (Eq. 4) and can be stated as follows: the river adjusts such that the vertical erosion rate is equal to the uplift or base-level lowering rate.

When control variables change, the river responds by adjusting its morphology – slope, width, and bed cover – to comply with both of the order principles. However, due to the different adjustment timescales, the path to a new steady-state morphology may be complex. As an example, consider a river at steady state when sediment supply increases. The river responds by depositing sediment and increasing stationary sediment mass (Eqs. 6 and 7). The increase in available stationary sediment increases entrainment rates (cf. Turowski and Hodge, 2017). Deposition continues until the river reaches a graded state in which sediment outflux from the considered reach is equal to sediment supply (Eqs. 6 and 7). At the same time, any change in stationary sediment directly affects bed cover (Eq. 32), and the immediate response of the stationary sediment mass is reflected in the short response times of bed cover (Eq. 42; Fig. 5). Changes in cover, in turn, affect both vertical and lateral incision rates, initiating slope and width adjustment. These adjust much more slowly than bed cover (Fig. 5) until the vertical erosion rate matches the uplift rate. Yet, adjustments in width and slope feed back into the sediment dynamics, e.g. by affecting transport capacity. Again, the river responds by depositing or entraining material to maintain grade. The mutual feedback continues until both order principles – grade and the erosional balance with matching incision and uplift rates – are satisfied.

With two order principles controlling bedrock channel adjustment, the river may be in a steady state with respect to one of them but not with the other. Because the adjustment timescale for cover is shortest (Fig. 5), with values that range from minutes to days, it can be expected that bedrock rivers are close to a graded state most of the times (cf. Phillips and Jerolmack, 2016). Given the long adjustment times for width and slope, this does not necessarily mean, however, that they are in a steady state with respect to bedrock elevation, where incision rate matches uplift rate.

4.4 What is a bedrock channel?

The considerations and arguments presented in this paper affect the conceptual view of a bedrock channel, and the use of relevant terminology. We can distinguish detachment-limited and transport-limited channels, which are identified with the two endmember descriptions focusing on the mass balance description of bedrock (detachment limited) and sediment (transport limited), respectively (cf. Shobe et al., 2017). For detachment-limited channels, we assume that the transport of sediment (Eq. 6) can be neglected; i.e. sediment transport does not significantly impact channel dynamics and morphology. Formally, this assumption is valid if sediment supply is very much smaller than transport capacity, implying that stationary sediment mass M_s∼0 or that sediment transport does not affect bedrock erosion. For transport-limited channels, we assume that bedrock incision can be neglected (Eq. 4). Formally, this assumption is valid if deposition or erosion has a negligible effect on the stationary sediment mass, in the mathematical limit as M_s goes to infinity. The latter point implies that entrainment or deposition of sediment does not significantly affect stationary sediment mass.

A formal definition of bedrock channels should fulfil a number of criteria (cf. Turowski et al., 2008b). First, the definition should comply with the intuition of field workers. Alluvial and bedrock channels are endmembers on a continuum of channel types, and, therefore, there will always be debated cases. But generally, most geomorphologists would agree on whether the particular river should be classified as an alluvial or bedrock river when seeing it in the field. Second, it should not rely on observations of field parameters that can change quickly, e.g. over a single flood. Third, a useful definition should not rely on parameters that cannot be measured. Fourth, it should not rely on theoretical concepts that are untested, untestable, or debated. Fifth, a definition rooted in the understanding of relevant processes or dynamics is preferable to one that relies solely on descriptions of morphology.

Bedrock channels, in general, have often been classified as detachment-limited channels, in which long-term sediment-supply is (much) smaller than long-term sediment transport capacity (e.g. Whipple, 2004; Whipple et al., 2013). Further, this condition is generally assumed to result in partial sediment cover and exposed bedrock on channel bed and banks. Bedrock exposure in the channel can easily be observed in the field, and it is therefore often used for channel classification (e.g. Montgomery et al., 1996; Tinkler and Wohl, 1998). A number of formal definitions of bedrock channels have been put forward based on these considerations. Exemplary, I will quote and discuss the most recent definition by Whipple et al. (2013):

Bedrock rivers may satisfy either or both of the following conditions: (1) the long-term capacity of the river to transport bedload (Q_c) exceeds the long-term supply of bedload (Q_s), resulting in generally sediment-starved conditions, significant rock exposure in bed and banks, and only thin, patchy, and temporary alluvial cover; or (2) the river is, over the long term (millennial to geologic timescales), actively incising through in-place rock.

Few geomorphologists would argue against the second part of the definition, although it may be difficult to assess this aspect in the field. Nevertheless, it is the first part of the definition that is relevant to the points made here, and which I reject based on the following general arguments and on the concepts developed in the present paper. First, the definition is theoretically laden in the sense that a theoretical concept is imposed and equated to a field observation. To my knowledge, no methods currently exist that allow one to reliably measure either long-term sediment supply or transport capacity. Even the inaccurate estimates that are currently possible need extensive field and modelling work, partly require strong assumptions, and are subject to large errors (e.g. Schneider et al., 2015). As such, the statement is not useful for the identification of bedrock channels in the field. Second, using mass balance arguments, I have demonstrated that bedrock channels adjust to a graded state. Unlike alluvial rivers, this does not imply that sediment supply is equal to transport capacity. Rather, the relationship between cover and the ratio of supply and capacity is modulated by the deposition and entrainment of stationary sediment mass, i.e. bed cover. The model by Turowski and Hodge (2017) predicts partial cover for sediment supply values larger than transport capacity in some parameter configurations. Similarly, the simulations by Inoue et al. (2016) predict partially covered bed for conditions where sediment supply equals transport capacity. This shows that – depending on the theoretical formulation and the relevant concepts, assumptions, and definitions – sediment supply values equal to or larger than transport capacity may be possible for bedrock channels. Third, even if the long-term sediment supply is lower than transport capacity, alluvial cover is not necessarily thin, patchy, or temporary, as is assumed in the definition. Rather, there can be thick, substantial, widespread, or persistent cover in the channel. For example, Shepherd (1972) and Fernández et al. (2019) documented persistent gravel bars in experimental meandering bedrock channels. Theoretical cover models (e.g. Hodge and Hoey, 2012; Sklar and Dietrich, 2004; Turowski and Hodge, 2017) predict substantial cover for certain sediment supply values that are smaller than the transport capacity. Experimental observations of runaway alluviation (e.g. Chatanantavet and Parker, 2008) provide evidence for this. Fourth, in a natural channel, sediment supply and discharge vary over timescales that are short in comparison to the adjustment timescales of channel width and slope. Upscaling discharge variability and sediment supply with a numerical model (Lague, 2010) showed that the channel bed is either fully covered or free of sediment for the majority of the time. Long-term mean cover values in his simulations exceeded a value of 0.5 in all cases, prohibiting the use of a detachment- or transport-limited approximation. Taken together, the arguments suggest that the connection between patchy, thin, and temporary alluvial cover and a ratio of sediment supply to transport capacity smaller than one is not tenable. As such, the definition, as proposed, is neither useful nor does it reflect current knowledge of bedrock channel dynamics. Turowski et al. (2008b) proposed an alternative definition, stating that “a bedrock channel cannot substantially widen, lower, or shift its bed without eroding bedrock”. This definition has been discussed and slightly altered by Meshkova et al. (2012). It does not stand in contradiction to field observations, current process knowledge, or newly emerged concepts, and it can be readily applied in the field (see Turowski et al., 2008b, for a discussion of relevant field criteria).

5 Conclusions

Bedrock channel dynamics are controlled by two dominant order principles. They adjust their morphology both to achieve grade, in which the sediment transport rate is constant along the stream, and to match incision rate to uplift or base-level lowering. The recognition of a steady state corresponding to one of these principles does not necessarily imply that the other has also been achieved. With minutes to days, the adjustment timescale for bed cover is short relative to the timescales for channel width and slope, and cover may be adjusted by changing supply conditions even over the duration of a single flood event. Thus, it can be expected that bedrock channels are close to a graded state most of the time. In the example calculations (Fig. 5), adjustment timescales for slope and width are of the order of thousands of years. This is shorter than the major cyclic variations of Earth's climate (e.g. Roe, 2006) or the typical timescales of mountain building. The results therefore suggest that many bedrock channels are also close to an erosional steady state, in which erosion rate is equal to uplift rate.

Appendix A: Deriving the Exner equation from the entrainment–deposition framework

Substituting Eq. (7) into Eq. (6) to eliminate entrainment and deposition rates, we obtain

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial M_{\mathrm{m}}}{\partial t}}=-{\displaystyle \frac{\partial q_{\mathrm{s}}}{\partial x}}-{\displaystyle \frac{\partial M_{\mathrm{s}}}{\partial t}}.\end{array}$$

We can rearrange to get

$$\begin{array}{}\text{(A2)}& {\displaystyle \frac{\partial \left(M_{\mathrm{m}}+M_{\mathrm{s}}\right)}{\partial t}}=-{\displaystyle \frac{\partial q_{\mathrm{s}}}{\partial x}}.\end{array}$$

We can define a total sediment mass per unit area, $M_{\mathrm{tot}}=M_{\mathrm{m}}+M_{\mathrm{s}}$, and divide by the sediment density, ρ_r(1−p), to obtain the Exner equation:

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{r}}(\mathrm{1}-p)}}{\displaystyle \frac{\partial M_{\mathrm{tot}}}{\partial t}}={\displaystyle \frac{\partial h_{\mathrm{s}}}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{r}}(\mathrm{1}-p)}}{\displaystyle \frac{\partial q_{\mathrm{s}}}{\partial x}}.\end{array}$$

Appendix B: Estimating the deflection angle

Assume that the bedload particle path through the channel follows a sinusoidal path with a wavelength equal to the gravel bar spacing and an amplitude A_bar:

$$\begin{array}{}\text{(B1)}& y=A_{\mathrm{bar}}\sin \left(\mathrm{2}\mathit{\pi }{\displaystyle \frac{x}{\mathit{\lambda }}}\right).\end{array}$$

Here, y denotes the distance in the cross-channel direction, with the channel centre line located at y=0, and x denotes the distance in the along-channel direction. The tangent of the angle γ is given by the derivative of Eq. (B1):

$$\begin{array}{}\text{(B2)}& \tan \left(\mathit{\gamma }\right)={\displaystyle \frac{\mathrm{d}y}{\mathrm{d}x}}=\mathrm{2}\mathit{\pi }{\displaystyle \frac{A_{\mathrm{bar}}}{\mathit{\lambda }}}\cos \left(\mathrm{2}\mathit{\pi }{\displaystyle \frac{x}{\mathit{\lambda }}}\right).\end{array}$$

We are interested in the deflection angle γ at the edge of the gravel bar, a distance W_covered (the covered part of the channel width) from the channel boundary, which corresponds to $y=W/\mathrm{2}-W_{\mathrm{covered}}$. Hence, at the corresponding x position x_edge

$$\begin{array}{}\text{(B3)}& \mathrm{2}\mathit{\pi }{\displaystyle \frac{x_{\mathrm{edge}}}{\mathit{\lambda }}}={\sin }^{-\mathrm{1}}.\end{array}$$

Here, sin ⁻¹ denotes the inverse sinus function. Combining Eqs. (B1)–(B3) and writing the path amplitude as a fraction, $f=\mathrm{2}A/W$, of the half channel width, we obtain

$$\begin{array}{}\text{(B4)}& \sin \left(\mathit{\gamma }\right)=\sin \left\{{\tan }^{-\mathrm{1}}\left[\mathit{\pi }{\displaystyle \frac{fW}{\mathit{\lambda }}}\cos \left({\sin }^{-\mathrm{1}}\left({\displaystyle \frac{\mathrm{1}}{f}}-{\displaystyle \frac{\mathrm{2}{W}_{\mathrm{covered}}}{fW}}\right)\right)\right]\right\}.\end{array}$$

Here, tan ⁻¹ denotes the inverse tangent function. Substituting Eq. (23) for λ and C for W_covered∕W (Eq. 22), we obtain

$$\begin{array}{}\text{(B5)}& \sin \left(\mathit{\gamma }\right)=\sin \left\{{\tan }^{-\mathrm{1}}\left[{\displaystyle \frac{f\mathit{\pi }}{k_{\mathrm{bar}}C}}\cos \left({\sin }^{-\mathrm{1}}\left({\displaystyle \frac{\mathrm{1}}{f}}-{\displaystyle \frac{\mathrm{2}C}{f}}\right)\right)\right]\right\}.\end{array}$$

Assuming f=1, a reasonable approximation for the square of Eq. (B5) (as it appears in all equations) is

$$\begin{array}{}\text{(B6)}& {\sin }^{\mathrm{2}}\left(\mathit{\gamma }\right)\approx \mathrm{1}-C.\end{array}$$

Appendix C: Deriving the response timescale of width adjustment using perturbation analysis

For the following analysis we assume that all parameters are kept constant apart from sediment supply, which varies sinusoidally over time. This choice allows one to obtain an analytical solution for the problem, and it does not affect the result for the timescale of transient adjustment. Sediment supply can then be written as the sum of the average supply $\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}$ and a perturbation term δQ_s. The variation of the latter is described with a sinusoidal oscillation around zero.

$$\begin{array}{}\text{(C1)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{s}}=\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}+\mathit{\delta }{Q}_{\mathrm{s}},\text{(C2)}& {\displaystyle }& {\displaystyle }\mathit{\delta }{Q}_{\mathrm{s}}=K\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{P}}\right).\end{array}$$

Here, K is a constant and P the period of the perturbation. Using linearized approximations to the differential equations (i.e. using first-order Taylor series to approximate non-linear functions), we then derive the width response to this perturbation, which can also be written as the sum of a time-independent term $\stackrel{\mathrm{‾}}{W}$ and a time-dependent term δW.

$$\begin{array}{}\text{(C3)}& W=\stackrel{\mathrm{‾}}{W}+\mathit{\delta }W\end{array}$$

To obtain an equation describing the time evolution of channel width, we combine Eqs. (24) and (27) to obtain

$$\begin{array}{}\text{(C4)}& {\displaystyle \frac{\mathrm{d}W}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{2}\mathit{\kappa }Yg}{k_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right){\sin }^{\mathrm{2}}\left(\mathit{\gamma }\right){\displaystyle \frac{Q_{\mathrm{s}}}{WC}}{\displaystyle \frac{\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}-\mathrm{1}\right)}{\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}-\mathrm{1}\right)+\mathrm{24}}}.\end{array}$$

We substitute the squared sine of the angle by Eq. (B6) to obtain

$$\begin{array}{}\text{(C5)}& {\displaystyle \frac{\mathrm{d}W}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{2}\mathit{\kappa }Yg}{k_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right){\displaystyle \frac{\mathrm{1}-C}{C}}{\displaystyle \frac{Q_{\mathrm{s}}}{W}}{\displaystyle \frac{\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}-\mathrm{1}\right)}{\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}-\mathrm{1}\right)+\mathrm{24}}}.\end{array}$$

To simplify the equation further, I make the assumption that excess transport stage θ∕θ_c rarely exceeds the value of 10. Then, we can approximate (cf. Auel et al., 2017a):

$$\begin{array}{}\text{(C6)}& {\displaystyle \frac{\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}-\mathrm{1}\right)}{\left(\frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}-\mathrm{1}\right)+\mathrm{24}}}\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{36}}}{\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}}.\end{array}$$

The width evolution equation is then

$$\begin{array}{}\text{(C7)}& {\displaystyle \frac{\mathrm{d}W}{\mathrm{d}t}}={\displaystyle \frac{\mathit{\kappa }Yg}{\mathrm{18}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right){\displaystyle \frac{\mathrm{1}-C}{C}}{\displaystyle \frac{Q_{\mathrm{s}}}{W}}{\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}}.\end{array}$$

Steady-state cover can be described with Eq. (32) (Turowski and Hodge, 2017).

$$\begin{array}{}\text{(C8)}& C=\left(\mathrm{1}-e^{-\frac{Q_{\mathrm{s}}}{M_{\mathrm{0}}UW}}\right){\displaystyle \frac{Q_{\mathrm{s}}}{Q_{\mathrm{t}}}}\end{array}$$

The bedload transport capacity can be written using Eq. (33) (see Turowski, 2018):

$$\begin{array}{}\text{(C9)}& {\displaystyle }{\displaystyle \frac{Q_{\mathrm{t}}}{W}}=K_{\mathrm{bl}}{Q}^{\mathrm{m}}{S}^{\mathrm{n}},\text{(C10)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}W}{\mathrm{d}t}}& =A{K}_{\mathrm{bl}}{Q}^{\mathrm{m}}{S}^{\mathrm{n}}{\left(\mathrm{1}-e^{-\frac{Q_{\mathrm{s}}}{M_{\mathrm{0}}UW}}\right)}^{-\mathrm{1}}(W{)}^{\mathit{\alpha }-\mathrm{1}}\\ & -A{Q}_{\mathrm{s}}(W{)}^{\mathit{\alpha }-\mathrm{2}},\end{array}\end{array}$$

with

$$\begin{array}{}\text{(C11)}& A={\displaystyle \frac{\mathit{\kappa }Yg}{\mathrm{18}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}}\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right){\displaystyle \frac{(gS{)}^{\frac{\mathit{\alpha }+\mathrm{1}}{\mathrm{2}}}{R}^{\frac{\mathrm{3}\mathit{\alpha }-\mathrm{1}}{\mathrm{2}}}}{k_{\mathrm{V}}\left(\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}-\mathrm{1}\right)D{\mathit{\theta }}_{\mathrm{c}}}}(Q{)}^{\mathrm{1}-\mathit{\alpha }}.\end{array}$$

Here, to reduce the number of parameters and reveal implicit dependencies, the Shields stress has been substituted using standard hydraulic scaling relations (Appendix D). The parameter α is a dimensionless constant that typically takes a value of 0.6, and it is a measure of roughness with the dimensions of length (see Nitsche et al., 2012). Next, Eqs. (C1) and (C3) are substituted into Eq. (C4), and expanded using first-order Taylor approximations of the form

$$\begin{array}{}\text{(C12)}& {\displaystyle }\begin{array}{rl}\left(B(\stackrel{\mathrm{‾}}{W}+\mathit{\delta }W{)}^{\mathit{\alpha }-\mathrm{1}}-\mathrm{1}\right)& \approx \left(B{\stackrel{\mathrm{‾}}{W}}^{\mathit{\alpha }-\mathrm{1}}-\mathrm{1}\right)\\ & +(\mathit{\alpha }-\mathrm{1})B{\stackrel{\mathrm{‾}}{W}}^{\mathit{\alpha }-\mathrm{2}}\mathit{\delta }W,\end{array}\text{(C13)}& {\displaystyle }\begin{array}{rl}\left(\mathrm{1}\right.& {\left.-e^{-\frac{Q_{\mathrm{s}}}{M_{\mathrm{0}}UW}}\right)}^{-\mathrm{1}}\approx {\left(\mathrm{1}-e^{-\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\right)}^{-\mathrm{1}}\\ & -{\left(\mathrm{1}-e^{-\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\right)}^{-\mathrm{2}}\left[{\displaystyle \frac{\mathrm{1}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}{e}^{-\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\mathit{\delta }{Q}_{\mathrm{s}}\right.\\ & \left.+{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}}}\left({\displaystyle \frac{\mathrm{1}}{U^{\mathrm{2}}\stackrel{\mathrm{‾}}{W}}}{\displaystyle \frac{\partial U}{\partial W}}+{\displaystyle \frac{\mathrm{1}}{U{\stackrel{\mathrm{‾}}{W}}^{\mathrm{2}}}}\right)e^{-\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\mathit{\delta }W\right].\end{array}\end{array}$$

After some algebra and dropping terms that are quadratic or cubic in the delta terms δQ_s and δW, we obtain

$$\begin{array}{}\text{(C14)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}\mathit{\delta }W}{\mathrm{d}t}}=& A{K}_{\mathrm{bl}}{\left(\mathrm{1}-e^{-\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\right)}^{-\mathrm{1}}(\stackrel{\mathrm{‾}}{W}{)}^{\mathit{\alpha }-\mathrm{1}}-A\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}(\stackrel{\mathrm{‾}}{W}{)}^{\mathit{\alpha }-\mathrm{2}}\\ & +\left[A{K}_{\mathrm{bl}}{Q}^{\mathrm{m}}{S}^{\mathrm{n}}{\left(\mathrm{1}-e^{-\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\right)}^{-\mathrm{1}}(\mathit{\alpha }-\mathrm{1}){\stackrel{\mathrm{‾}}{W}}^{\mathit{\alpha }-\mathrm{2}}\right.\\ & -A{K}_{\mathrm{bl}}{Q}^{\mathrm{m}}{S}^{\mathrm{n}}{\left(\mathrm{1}-e^{-\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\right)}^{-\mathrm{2}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}}}\\ & \left({\displaystyle \frac{\mathrm{1}}{U^{\mathrm{2}}\stackrel{\mathrm{‾}}{W}}}{\displaystyle \frac{\partial U}{\partial W}}+{\displaystyle \frac{\mathrm{1}}{U{\stackrel{\mathrm{‾}}{W}}^{\mathrm{2}}}}\right)e^{-\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}{\stackrel{\mathrm{‾}}{W}}^{\mathit{\alpha }-\mathrm{1}}\\ & \left.-A(\mathit{\alpha }-\mathrm{2})\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}{\stackrel{\mathrm{‾}}{W}}^{\mathit{\alpha }-\mathrm{3}}\right]\mathit{\delta }W+\left[A{K}_{\mathrm{bl}}{Q}^{\mathrm{m}}{S}^{\mathrm{n}}\right.\\ & {\left(\mathrm{1}-e^{-\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\right)}^{-\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{M_{\mathrm{0}}U}}{e}^{-\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}{\stackrel{\mathrm{‾}}{W}}^{\mathit{\alpha }-\mathrm{2}}\\ & \left.-A{\stackrel{\mathrm{‾}}{W}}^{\mathit{\alpha }-\mathrm{2}}\right]\mathit{\delta }{Q}_{\mathrm{s}}.\end{array}\end{array}$$

Resubstituting for cover, particle speed, and so on, we obtain

$$\begin{array}{}\text{(C15)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}\mathit{\delta }W}{\mathrm{d}t}}& =A{K}_{\mathrm{bl}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{\stackrel{\mathrm{‾}}{C}{Q}_{\mathrm{t}}}}(\stackrel{\mathrm{‾}}{W}{)}^{\mathit{\alpha }-\mathrm{1}}-A\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}(\stackrel{\mathrm{‾}}{W}{)}^{\mathit{\alpha }-\mathrm{2}}+A\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}{\stackrel{\mathrm{‾}}{W}}^{\mathit{\alpha }-\mathrm{3}}\\ & \left[{\displaystyle \frac{(\mathit{\alpha }-\mathrm{1})}{\stackrel{\mathrm{‾}}{C}}}-{\displaystyle \frac{\mathrm{3}}{\mathrm{2}\stackrel{\mathrm{‾}}{C}}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\left({\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{C}}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{\stackrel{\mathrm{‾}}{Q_{\mathrm{t}}}}}-\mathrm{1}\right)-(\mathit{\alpha }-\mathrm{2})\right]\\ & \mathit{\delta }W+A{\stackrel{\mathrm{‾}}{W}}^{\mathit{\alpha }-\mathrm{2}}\left[{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\left({\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{C}}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{\stackrel{\mathrm{‾}}{Q_{\mathrm{t}}}}}-\mathrm{1}\right)-\mathrm{1}\right]\mathit{\delta }{Q}_{\mathrm{s}}.\end{array}\end{array}$$

Next, Eq. (C2) is substituted into Eq. (C13) to obtain a differential equation of the form

$$\begin{array}{}\text{(C16)}& {\displaystyle \frac{\mathrm{d}\mathit{\delta }W}{\mathrm{d}t}}=K_{\mathrm{1}}\left[K_{\mathrm{2}}+K\sin \left\{{\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{P}}\right\}+K_{\mathrm{3}}\mathit{\delta }W\right].\end{array}$$

The general solution to Eq. (C14) is

$$\begin{array}{}\text{(C17)}& \begin{array}{rl}\mathit{\delta }W=& {\displaystyle \frac{K_{\mathrm{1}}{K}_{\mathrm{2}}\left(\frac{P}{\mathrm{2}\mathit{\pi }}\right)}{K_{\mathrm{1}}^{\mathrm{2}}{K}_{\mathrm{3}}^{\mathrm{2}}{\left(\frac{P}{\mathrm{2}\mathit{\pi }}\right)}^{\mathrm{2}}+\mathrm{1}}}\sqrt{\left(K_{\mathrm{1}}^{\mathrm{2}}{K}_{\mathrm{3}}^{\mathrm{2}}{\left({\displaystyle \frac{P}{\mathrm{2}\mathit{\pi }}}\right)}^{\mathrm{2}}+\mathrm{1}\right)}\sin \\ & \left\{{\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{P}}+\mathit{\phi }\right\}+{\displaystyle \frac{K{K}_{\mathrm{1}}^{\mathrm{2}}{K}_{\mathrm{3}}^{\mathrm{2}}{\left(\frac{P}{\mathrm{2}\mathit{\pi }}\right)}^{\mathrm{2}}+K}{K_{\mathrm{1}}^{\mathrm{2}}{K}_{\mathrm{3}}^{\mathrm{2}}{\left(\frac{P}{\mathrm{2}\mathit{\pi }}\right)}^{\mathrm{2}}+\mathrm{1}}}\\ & +c_{\mathrm{1}}\exp \left\{K_{\mathrm{1}}{K}_{\mathrm{3}}t\right\}.\end{array}\end{array}$$

Here, c₁ is the integrative constant and φ is a phase shift of the width response to the perturbation in sediment supply. The exponential term describes transient adjustment to the steady state and can be used to obtain the response time.

$$\begin{array}{}\text{(C18)}& T_{\mathrm{W}}=-{\displaystyle \frac{\mathrm{1}}{K_{\mathrm{1}}{K}_{\mathrm{3}}}}.\end{array}$$

Collecting the terms, we obtain

$$\begin{array}{}\text{(C19)}& \begin{array}{rl}{T}_{\mathrm{W}}& ={\displaystyle \frac{\mathrm{18}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}{\mathit{\kappa }Yg}}{\displaystyle \frac{k_{\mathrm{V}}gD{\mathit{\theta }}_{\mathrm{c}}}{{\left(gS\right)}^{\frac{\mathit{\alpha }+\mathrm{1}}{\mathrm{2}}}{R}^{\frac{\mathrm{3}\mathit{\alpha }-\mathrm{1}}{\mathrm{2}}}}}{\displaystyle \frac{Q^{\mathit{\alpha }-\mathrm{1}}}{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}}{\stackrel{\mathrm{‾}}{W}}^{\mathrm{3}-\mathit{\alpha }}\\ & {\left({\displaystyle \frac{\mathrm{3}}{\mathrm{2}\stackrel{\mathrm{‾}}{C}}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\left({\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{C}}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{Q_{\mathrm{t}}}}-\mathrm{1}\right)-{\displaystyle \frac{(\mathit{\alpha }-\mathrm{1})}{\stackrel{\mathrm{‾}}{C}}}+(\mathit{\alpha }-\mathrm{2})\right)}^{-\mathrm{1}}.\end{array}\end{array}$$

Equation (C17) is considerably simpler in terms of shear stress

$$\begin{array}{}\text{(C20)}& \begin{array}{rl}{T}_{\mathrm{W}}& ={\displaystyle \frac{\mathrm{18}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}{\mathit{\kappa }Y\left(\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}-\mathrm{1}\right)g}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{c}}}{\mathit{\theta }}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{W}}^{\mathrm{2}}}{Q_{\mathrm{t}}}}\left({\displaystyle \frac{\mathrm{3}}{\mathrm{2}\stackrel{\mathrm{‾}}{C}}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{M_{\mathrm{0}}U\stackrel{\mathrm{‾}}{W}}}\left({\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{C}}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{Q_{\mathrm{t}}}}-\mathrm{1}\right)\right.\\ & {\left.-{\displaystyle \frac{(\mathit{\alpha }-\mathrm{1})}{\stackrel{\mathrm{‾}}{C}}}+(\mathit{\alpha }-\mathrm{2})\right)}^{-\mathrm{1}}.\end{array}\end{array}$$

In the linear cover approximation (cover-dominated limit; see Turowski, 2018), we have

$$\begin{array}{}\text{(C21)}& \stackrel{\mathrm{‾}}{C}={\displaystyle \frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{Q_{\mathrm{t}}}}.\end{array}$$

Thus, Eq. (C18) becomes

$$\begin{array}{}\text{(C22)}& T_{\mathrm{W},\mathrm{cover}}={\displaystyle \frac{\mathrm{18}{k}_{\mathrm{v}}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}}{\mathit{\kappa }Y\left(\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}-\mathrm{1}\right)g}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{c}}}{\mathit{\theta }}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{W}}^{\mathrm{2}}}{Q_{\mathrm{t}}}}\left({\displaystyle \frac{\mathrm{1}}{(\mathit{\alpha }-\mathrm{2})\frac{\stackrel{\mathrm{‾}}{Q_{\mathrm{s}}}}{Q_{\mathrm{t}}}-(\mathit{\alpha }-\mathrm{1})}}\right).\end{array}$$

Appendix D: Writing shear stress and bedload speed in terms of discharge

The development here is equivalent to that given by Turowski (2018).

The reach-averaged Shields stress θ is defined by

$$\begin{array}{}\text{(D1)}& \mathit{\theta }={\displaystyle \frac{\mathit{\tau }}{\left({\mathit{\rho }}_{\mathrm{s}}-\mathit{\rho }\right)gD}}.\end{array}$$

Here, τ is the shear stress, given by the DuBoys equation

$$\begin{array}{}\text{(D2)}& \mathit{\tau }=\mathit{\rho }gHS.\end{array}$$

The continuity equation for water flow is

$$\begin{array}{}\text{(D3)}& Q=WHV.\end{array}$$

There are a number of different equations available to compute water flow velocity V. For mountain streams, a discharge-based variable power flow resistance equation has been found to be a good description of available data (Ferguson, 2007; Nitsche et al., 2012):

$$\begin{array}{}\text{(D4)}& V=k_{\mathrm{V}}(gS{)}^{\frac{\mathrm{1}-\mathit{\alpha }}{\mathrm{2}}}{R}^{\frac{\mathrm{1}-\mathrm{3}\mathit{\alpha }}{\mathrm{2}}}{\left({\displaystyle \frac{Q}{W}}\right)}^{\mathit{\alpha }}.\end{array}$$

Here, R is a measure of bed roughness with dimensions of length, e.g. the standard deviation of the bed surface (e.g. Nitsche et al., 2012), and k_V≈1 and α≈0.6 are constants. Combining Eqs. (C2)–(C4), shear stress can be written as

$$\begin{array}{}\text{(D5)}& \mathit{\tau }={\displaystyle \frac{\mathit{\rho }}{k_{\mathrm{V}}}}(gS{)}^{\frac{\mathit{\alpha }+\mathrm{1}}{\mathrm{2}}}{R}^{\frac{\mathrm{3}\mathit{\alpha }-\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{Q}{W}}\right)}^{\mathrm{1}-\mathit{\alpha }}.\end{array}$$

The Shields stress is thus given by

$$\begin{array}{}\text{(D6)}& \mathit{\theta }={\displaystyle \frac{(gS{)}^{\frac{\mathit{\alpha }+\mathrm{1}}{\mathrm{2}}}{R}^{\frac{\mathrm{3}\mathit{\alpha }-\mathrm{1}}{\mathrm{2}}}}{k_{\mathrm{V}}\left(\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}-\mathrm{1}\right)gD}}{\left({\displaystyle \frac{Q}{W}}\right)}^{\mathrm{1}-\mathit{\alpha }}.\end{array}$$

The downstream bedload velocity arises in the cover relation (Eq. 32), and it can be written as

$$\begin{array}{}\text{(D7)}& U=\mathrm{1.46}{\left(\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right)gD\right)}^{\mathrm{1}/\mathrm{2}}{\left({\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{c}}}}-\mathrm{1}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

In terms of discharge, this evaluates to

$$\begin{array}{}\text{(D8)}& \begin{array}{rl}U=& \mathrm{1.46}{\left(\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}}-\mathrm{1}\right)gD\right)}^{\mathrm{1}/\mathrm{2}}\\ & {\left({\displaystyle \frac{(gS{)}^{\frac{\mathit{\alpha }+\mathrm{1}}{\mathrm{2}}}{R}^{\frac{\mathrm{3}\mathit{\alpha }-\mathrm{1}}{\mathrm{2}}}}{k_{\mathrm{V}}{\mathit{\theta }}_{\mathrm{c}}\left(\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathit{\rho }}-\mathrm{1}\right)gD}}{\left({\displaystyle \frac{Q}{W}}\right)}^{\mathrm{1}-\mathit{\alpha }}-\mathrm{1}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}\end{array}$$

Appendix E: Notation

Abar	Bar amplitude (m)
Acover	Covered bed area (m2)
Atot	Total bed area (m2)
Aw	Actively eroding channel wall area (m2)
a	Scaling exponent, d−A
B	Constant in non-linear wave equation, describing slope development (m s−1)
b	Scaling exponent, β−A
C	Fraction of covered bed
CSS	Steady-state cover
c	Scaling exponent, Q−A
cs	Celerity of non-linear wave equation, describing slope development (m s−1)
d	Sideward deflection length scale, reach (m)
D	Sediment deposition rate per bed area (kg m−2 s−1)
D50	Median grain size (m)
e	Base of the natural logarithm
E	Sediment entrainment rate per bed area (kg m−2 s−1)
f	Bedload path amplitude as fraction of channel width
FCD	Cover-dependent term in the lateral erosion equation
FT	Tools factor in the lateral erosion equation (kg s−1)
g	Acceleration due to gravity (m s−2)
hb	Bedrock elevation (m)
hs	Sediment elevation (m)
H	Water depth (m)
I	Vertical erosion rate (m s−1)
Kbl	Bedload transport efficiency (kg m−3m s−m)
ke	Erodibility in stream power model (m1−3m s1−m)
kh	Hydrology coefficient (m3−2c s−1)
kLE	Lateral erosion coefficient (m kg−1)
ks	Steepness index (m2θ)
ktools	Lumped constant, tools-dominated channel slope
kv	Rock erodibility coefficient
kV	Velocity coefficient (m2α)
L	Reach length (m)
M0	Minimum mass per area necessary to cover the bed (kg m−2)
Mm	Mobile sediment mass (kg m−2)
Ms	Stationary sediment mass (kg m−2)
m	Discharge exponent in bedload equation
m′	Discharge exponent in the stream power model
n	Slope exponent in bedload equation
n′	Slope exponent in the stream power model
qs	Mass sediment transport rate per unit width (kg ms−1)
qt	Mass sediment transport capacity per unit width (kg ms−1)
Q	Water discharge (m3 s−1)
$Q_{\mathrm{c}}^{*}$	Relative sediment supply at the critical cover
Qs	Upstream sediment mass supply (kg s−1)
$Q_{\mathrm{s}}^{*}$	Relative sediment supply; sediment transport rate over transport capacity
Qt	Mass sediment transport capacity (kg s−1)
R	Bed roughness length scale (m)
S	Channel bed slope
Scover	Channel bed slope predicted in the cover-dominated approximation
Stools	Channel bed slope predicted in the tools-dominated approximation
SV	Valley slope

TC	Timescale of cover adjustment (s)
TN	Timescale of width adjustment, for a narrowing channel (s)
TS	Timescale of slope adjustment (s)
TW	Timescale of width adjustment, for a widening channel (s)
U	Bedload speed (m s−1)
V	Water flow velocity (m s−1)
W	Channel width (m)
Wcover	Covered length within the channel width (m)
Wuncover	Uncovered length within the channel width (m)
Wss	Steady-state channel width (m)
x	Dimensional streamwise spatial coordinate (m)
Y	Young's modulus of the bedrock (kg m−1 s−2)
α	Scaling exponent, V−Q
β	Fraction of sediment transported as bedload
θ	Shields stress
θc	Critical Shields stress for the onset of sediment motion
λ	Bar wavelength (m)
ρ	Density of water (kg m−3)
ρbulk	Bulk density of sediment (kg m−3)
Pr	Density of bedrock (kg m−3)
ρs	Density of sediment (kg m−3)
σT	Rock tensile strength (kg m−1 s−2)
κ	Lumped constant, width evolution equation (m−2)
κbar	Bar aspect ratio
κC	Coefficient in the cover term of width evolution
κT	Coefficient in the tools term of width evolution
ξ	Average sediment thickness above the bedrock (m)
τ	Bed shear stress (N m−2)
τc	Critical bed shear stress at the onset of bedload motion (N m−2)