A probabilistic framework for the cover effect in bedrock erosion

Abstract

The cover effect in fluvial bedrock erosion is a major control on bedrock channel morphology and long-term channel dynamics. Here, we suggest a probabilistic framework for the description of the cover effect that can be applied to field, laboratory, and modelling data and thus allows the comparison of results from different sources. The framework describes the formation of sediment cover as a function of the probability of sediment being deposited on already alluviated areas of the bed. We define benchmark cases and suggest physical interpretations of deviations from these benchmarks. Furthermore, we develop a reach-scale model for sediment transfer in a bedrock channel and use it to clarify the relations between the sediment mass residing on the bed, the exposed bedrock fraction, and the transport stage. We derive system timescales and investigate cover response to cyclic perturbations. The model predicts that bedrock channels can achieve grade in steady state by adjusting bed cover. Thus, bedrock channels have at least two characteristic timescales of response. Over short timescales, the degree of bed cover is adjusted such that the supplied sediment load can just be transported, while over long timescales, channel morphology evolves such that the bedrock incision rate matches the tectonic uplift or base-level lowering rate.

Introduction

Bedrock channels are shaped by erosion caused by countless impacts of the sediment particles they carry along their bed (Beer and Turowski, 2015; Cook et al., 2013; Sklar and Dietrich, 2004). There are feedbacks between the evolving channel morphology, the bedload transport, and the hydraulics (e.g. Finnegan et al., 2007; Johnson and Whipple, 2007; Wohl and Ikeda, 1997). Impacting bedload particles driven forward by the fluid forces erode and therefore shape the bedrock bed. In turn, the morphology of the channel determines the pathways of both sediment and water, and the forces which the latter exerts on the former, and thus sets the stage for the entrainment and deposition of the sediment (Hodge and Hoey, 2016). Sediment particles play a key role in this erosion process; they provide the tools for erosion and also determine where bedrock is exposed such that it can be worn away by impacting particles (Gilbert, 1877; Sklar and Dietrich, 2004).

The importance of the cover effect – that a stationary layer of gravel can shield the bedrock from bedload impacts – has by now been firmly established in a number of field and laboratory studies (e.g. Chatanantavet and Parker, 2008; Finnegan et al., 2007; Hobley et al., 2011; Johnson and Whipple, 2007; Turowski and Rickenmann, 2009; Turowski et al., 2008; Yanites et al., 2011). Sediment cover is generally modelled with generic relationships that predict the decrease in the fraction of exposed bedrock area $A^{*}$ with the increase in the relative sediment supply $Q_{\text{s}}^{*}$, usually defined as the ratio of sediment supply to transport capacity. Based on laboratory experiments and simple modelling, Turowski and Bloem (2016) argued that the focus on covered area is generally justified on the reach scale and that the erosion of bedrock under a thin sediment cover can be neglected. However, the behaviour of sediment cover under flood conditions is currently unknown and the assumption that the cover distribution at low flow is representative of that at high flow may not be justified (cf. Beer et al., 2016; Turowski et al., 2008).

The most commonly used function to describe the cover effect is the linear decline (Sklar and Dietrich, 1998), which is the simplest function connecting the steady-state endmembers of an empty bed when relative sediment supply $Q_{\text{s}}^{*}=\mathrm{0}$ and full cover when $Q_{\text{s}}^{*}=\mathrm{1}$: $${\displaystyle }{A}^{*}=\left\{\begin{array}{cc}\mathrm{1}-Q_{\text{s}}^{*}& \begin{array}{cc}\mathrm{for}& Q_{\text{s}}^{*}<\mathrm{1}\end{array}\\ \mathrm{0}& \mathrm{otherwise}.\end{array}\right.$$

In contrast, the exponential cover function arises under the assumption that particle deposition is equally likely for each part of the bed, whether it is covered or not (Turowski et al., 2007). $${\displaystyle }{A}^{*}=\left\{\begin{array}{cc}{e}^{-Q_{\text{s}}^{*}}& \begin{array}{cc}\mathrm{for}& Q_{\text{s}}^{*}<\mathrm{1}\end{array}\\ \mathrm{0}& \mathrm{otherwise}.\end{array}\right.$$ Here, $e$ is the base of the natural logarithm.

Hodge and Hoey (2012) obtained both the linear and the exponential functions using a cellular automaton (CA) model that modulated grain entrainment probabilities by the number of neighbouring grains. However, consistent with laboratory flume data, the same model also produced other behaviours under different parameterizations. One alternative behaviour is runaway alluviation, which was attributed by Chatanantavet and Parker (2008) to the differing roughness of bedrock and alluvial patches. Due to a decrease in flow velocity, an increase in surface roughness, and differing grain geometry, the likelihood of deposition is higher over bed sections covered by alluvium compared to smooth, bare bedrock sections (Hodge et al., 2011). This can lead to rapid alluviation of the entire bed once a minimum fraction has been covered. The relationship between sediment flux and cover is also affected by the bedrock morphology; flume experiments have demonstrated that on a non-planar bed, the location of sediment cover is driven by bed topography and hydraulics (e.g. Finnegan et al., 2007; Inoue et al., 2014). Johnson and Whipple (2007) observed that stable patches of alluvium tend to form in topographic lows such as potholes and at the bottom of slot canyons, whereas Hodge and Hoey (2016) found that local flow velocity also controls sediment cover location.

The relationship between roughness, bed cover, and incision was explored in a number of recent numerical modelling studies. Nelson and Seminara (2011, 2012) were one of the first to model the impact that the differing roughness of bedrock and alluvial areas has on sediment patch stability. Zhang et al. (2015) formulated a macro-roughness cover model, in which sediment cover is related to the ratio of sediment thickness to bedrock macro-roughness. Aubert et al. (2016) directly simulated the dynamics of particles in a turbulent flow and obtained both linear and exponential cover functions. Johnson (2014) linked sediment transport and cover to bed roughness in a reach-scale model. Using a model formulation similar to that of Nelson and Seminara (2011), Inoue et al. (2016) reproduced bar formation and sediment dynamics in bedrock channels. All of these studies used slightly different approaches and mathematical formulations to describe alluvial cover, making a direct comparison difficult.

Over timescales including multiple floods, the variability in sediment supply is also important (e.g. Turowski et al., 2013). Lague (2010) used a model formulation in which cover was written as a function of the average sediment depth to upscale daily incision processes to long timescales. He found that over the long term, cover dynamics are largely independent of the precise formulation on the process scale and are rather controlled by the magnitude–frequency distribution of discharge and sediment supply. Using the CA model of Hodge and Hoey (2012), Hodge (2017) found that, when sediment supply was very variable (alternating large pulses with no sediment supply), the amount of sediment cover was primarily determined by the recent supply history, rather than by the relationships identified under constant sediment supply.

So far, it has been somewhat difficult to compare and discuss the different cover functions obtained from theoretical considerations, numerical models, and experiments, since a unifying framework and clear benchmark cases have been missing. Here, we propose such a framework and develop type cases linked to physical considerations of the flow hydraulics and sediment erosion and deposition. We show how this framework can be applied to data from a published model (Hodge and Hoey, 2012). Furthermore, we develop a reach-scale erosion–deposition model that allows the dynamic modelling of cover and prediction of steady states. Thus, we clarify the relationship between cover, deposited mass, and relative sediment supply. As part of this model framework, we investigate the response time of a channel to a change in sediment input, which we illustrate using data from a natural channel.

A probabilistic framework

Development

Here we build on the arguments put forward by Turowski et al. (2007) and Turowski (2009). Consider a bedrock bed on which sediment particles are distributed. We can view the deposition of each particle as a random process, and each area element on the bed surface can be assigned a probability for the deposition of a particle. When assuming that a given number of particles are distributed on the bed, the mean behaviour of the exposed area $A^{*}$ can be calculated from the following equation (Fig. 1): $${\displaystyle }\text{d}{A}^{*}=-P\left(A^{*},,,M_{\text{s}}^{*},,,\mathrm{\dots }\right)\text{d}{M}_{\text{s}}^{*}.$$ $P$ is the probability that a given particle is deposited on the exposed part of the bed, which here is a function of the fraction of exposed area ($A^{*})$ and a dimensionless mass of particles on the bed per area ($M_{\text{s}}^{*}$, explained below), but which can be expected to also be a function of the relative sediment supply, the bed topography and roughness, the particle size, the local hydraulics, or other control variables. $M_{\text{s}}^{*}$ is a dimensionless mass equal to the total mass of the particles residing on the bed per area, which is suitably normalized. A suitable mass for normalization is the minimum mass required to cover a unit area, $M_{\mathrm{0}}$, as will become clear later. The minus sign is introduced because the fraction of the exposed area reduces as $M_{\text{s}}^{*}$ increases. As most previous relationships are expressed in terms of relative sediment supply $Q_{\text{s}}^{*}$, the relation of $M_{\text{s}}^{*}$ to $Q_{\text{s}}^{*}$ will be discussed later.

Cartoon illustration of a bed partially covered by sediment. For purpose of illustration, the bed is divided into a square raster, with each pixel of the size of a single grain. For a given number of particles in the area of the bed of interest, the exposed area fraction of the bed is dependent on the distribution of particles. Grains that sit on top of other grains do not contribute to cover. The probability that a new grain is deposited on uncovered bed is given by $P$ (Eq. 3).

[Figure omitted. See PDF]

We can make some general statements about $P$. First, $P$ is defined for the range $\mathrm{0}\le A^{*}\le \mathrm{1}$ and undefined elsewhere. Second, $P$ takes values between 0 and 1 for $\mathrm{0}\le A^{*}\le \mathrm{1}$. Third, $P(A^{*}=\mathrm{0})=\mathrm{0}$ and $P(A^{*}=\mathrm{1})=\mathrm{1}$. Note that $P$ is not a distribution function and therefore does not need to integrate to 1. Neither does it have to be continuous and differentiable everywhere.

For purpose of illustration, we will next discuss two simple forms of the probability function $P$ that lead to the linear and exponential forms of the cover effect. First, consider the case that all particles are always deposited on exposed bedrock. In this case, formally, to keep with the conditions stated above, we define $P=\mathrm{1}$ for $\mathrm{0}<A^{*}\le \mathrm{1}$ and $P=\mathrm{0}$ for $A^{*}=\mathrm{0}$. Thus, we can write $${\displaystyle }\begin{array}{ll}\text{d}{A}^{*}=-\text{d}{M}_{\text{s}}^{*}& \begin{array}{ll}\mathrm{for}& \mathrm{0}<A^{*}\le \mathrm{1}\end{array}\\ \text{d}{A}^{*}=\mathrm{0}& \begin{array}{cc}\mathrm{for}& A^{*}=\mathrm{0}\end{array}\end{array}.$$

Integrating, we obtain $${\displaystyle }{A}^{*}=-M_{\text{s}}^{*}+C,$$ where the constant of integration $C$ is found to equal 1 by using the condition $A^{*}(M_{\text{s}}^{*}=\mathrm{0})=\mathrm{1}$. Thus, we obtain a linear cover function. Note that the linear cover function gives a theoretical lower bound for the amount of cover: it arises when all available sediment always falls on uncovered ground, and thus no additional sediment is available that could facilitate quicker alluviation. In essence, this is a mass conservation argument. Now it is obvious why $M_{\mathrm{0}}$ is a convenient way to normalize: in plots of $A^{*}$ against $M_{\text{s}}^{*}$, we obtain a triangular region bounded by the points [0,1], [0,0], and [1,0] in which the cover function cannot exist (Fig. 2).

Similarly to above, if we set $P$ to a constant value, $k$, smaller than 1 for $\mathrm{0}<A^{*}\le \mathrm{1}$, we obtain $${\displaystyle }{A}^{*}=\mathrm{1}-k{M}_{\text{s}}^{*}.$$

It is clear that the assumption of $P=k$ is physically unrealistic because it implies that the probability of deposition on exposed ground is independent of the amount of uncovered bedrock. Especially when $A^{*}$ is close to 0, it seems unlikely that, say, 90 % of the sediment always falls on uncovered ground. A more realistic assumption is that the probability of deposition on uncovered ground is independent of location and other possible controls but is equal to the fraction of exposed area; i.e. $P=A^{*}$. In a probabilistic sense, this is also the simplest plausible assumption one can make. Then $${\displaystyle }\text{d}{A}^{*}=-A^{*}\text{d}{M}_{\text{s}}^{*},$$ giving upon integration $${\displaystyle }{A}^{*}=e^{-M_{\text{s}}^{*}}.$$

The argument used here to obtain the exponential cover effect in Eq. (8) essentially corresponds to the one given by Turowski et al. (2007). Since this case presents the simplest plausible assumption, we will use it as a benchmark case, to which we will compare other possible functional forms of $P$.

In principle, the probability function $P$ can be varied to account for various processes that make deposition more likely either on already covered ground by decreasing $P$ for the appropriate range of $A^{*}$ from the benchmark case $P=A^{*}$ or on uncovered ground by increasing $P$ from the benchmark case $P=A^{*}$. As has been identified previously (Chatanantavet and Parker, 2008; Hodge and Hoey 2012), roughness feedbacks to the flow can cause either case, depending on whether subsequent deposition is adjacent to or on top of existing sediment patches. In the former case, particles residing on an otherwise bare bedrock bed act as obstacles for moving particles and create a low-velocity wake zone in the downstream direction. In addition, particles residing on other single particles are unstable and stacks of particles are unlikely. Hence, newly arriving particles tend to deposit either upstream or downstream of stationary particles and the probability is generally higher for deposition on uncovered ground than in the benchmark case. In the latter case, larger patches of stationary particles increase the surface roughness of the bed, thus decreasing the local flow velocity and stresses, making deposition on the patch more likely. In this way, the probability of deposition on already covered bed is increased in comparison to the benchmark case.

A simple functional form that can be used to take into account either one of these two effects is a power law dependence of $P$ on $A^{*}$, taking the form $P=A^{{*}^{\mathit{\alpha }}}$ (Fig. 2a). Then, the cover function becomes (Fig. 2b) $${\displaystyle }{A}^{*}={\left(\mathrm{1},-,\left(\mathrm{1},-,\mathit{\alpha }\right),M_{\text{s}}^{*}\right)}^{\frac{\mathrm{1}}{\mathrm{1}-\mathit{\alpha }}}.$$

Here, the probability of deposition on uncovered ground is increased in comparison to the benchmark exponential case if $\mathrm{0}<\mathit{\alpha }<\mathrm{1}$ and decreased if $\mathit{\alpha }>\mathrm{1}$.

A convenient and flexible way to parameterize $P\left(A^{*}\right)$ in general is the cumulative version of the Beta distribution, given by $${\displaystyle }P\left(A^{*}\right)=B\left(A^{*},;,a,,,b\right).$$

Here, $B(A^{*};a,b)$ is the regularized incomplete Beta function with two shape parameters $a$ and $b$, which are both real positive numbers, defined by $${\displaystyle }B\left(A^{*},;,a,,,b\right)={\displaystyle \frac{{\int }_{\mathrm{0}}^{A^{*}}{y}^{a-\mathrm{1}}{\left(\mathrm{1},-,y\right)}^{b-\mathrm{1}}\text{d}y}{{\int }_{\mathrm{0}}^{\mathrm{1}}{y}^{a-\mathrm{1}}{\left(\mathrm{1},-,y\right)}^{b-\mathrm{1}}\text{d}y}}.$$

Here, $y$ is a dummy variable. With suitable choices for $a$ and $b$, cover functions resembling the exponential ($a=b=$1), the linear form ($a=\mathrm{0}$, $b>\mathrm{0}$), and the power law form ($a\gg b$ or $a\ll b)$ can be retrieved. Wavy functions are also a possibility (Fig. 3); thus, both of the roughness effects described above can be modelled in a single scenario. Unfortunately, the integral necessary to obtain $A^{*}\left(M_{\text{s}}^{*}\right)$ does not give a closed-form analytical solution and needs to be computed numerically.

In principle, a suitable function $P$ could also be defined to account for the influence of bed topography on sediment deposition. Such a function is likely dependent on the details of the particular bed, hydraulics, and sediment flow paths in a complex way and needs to be mapped out experimentally.

(a) Various examples for the probability function $P$ as a function of bedrock exposure $A^{*}$. (b) Corresponding analytical solutions for the cover function between $A^{*}$ and dimensionless sediment mass $M_{\text{s}}^{*}$ using Eqs. (6), (7), and (9). Grey shading depicts the area where the cover function cannot run due to conservation of mass.

[Figure omitted. See PDF]

Examples for the use of the regularized incomplete Beta function (Eq. 11) to parameterize $P$, using various values for the shape parameters $a$ and $b$. The choice $a=b=\mathrm{1}$ gives a dependence that is equivalent to the exponential cover function. Grey shading depicts the area where the cover function cannot run due to conservation of mass.

[Figure omitted. See PDF]

Example of application using model data

To illustrate how the framework can be used, we apply it to data obtained from the CA model developed by Hodge and Hoey (2012). The CA model reproduces the transport of individual sediment grains over a smooth bedrock surface. In each time step, the probability of a grain being entrained is a function of the number of neighbouring grains. If five or more of the eight neighbouring cells contain grains, then the grain has a probability of entrainment $p_{\text{c}}$; otherwise it has the probability $p_{\text{i}}$. In most model runs $p_{\text{c}}$ was set to a value less than that of $p_{\text{i}}$, thus accounting for the impact of sediment cover in decreasing local shear stress (through increased flow resistance) and increasing the critical entrainment shear stress for grains (via lower grain exposure and increased pivot angles). Thus, in the model, grain scale dynamics of entrainment are varied by adjusting the values of $p_{\text{i}}$ and $p_{\text{c}}$. This has a direct effect on the reach-scale distribution of cover, which is captured by our $P$ function (Eq. 3).

The model is run with a domain that is 100 cells wide by 1000 cells long, with each cell having the same area as a grain. Up to four grains can potentially be entrained from each cell in a time step, limiting the maximum sediment flux. In each time step random numbers and the probabilities are used to select the grains that are entrained, which are then moved a step length of 10 cells downstream and deposited. Model results are insensitive to the step length. A fixed number of grains are also supplied to the upstream end of the model domain. A smoothing algorithm is applied to prevent unrealistically tall piles of grains developing in cells if there are far fewer grains in adjacent cells. After around 500 time steps the model typically reaches a steady-state condition in which the number of grains supplied to and leaving the model domain are equal. Sediment cover is measured in a downstream area of the model domain and is defined as grains that are not entrained in a given time step. Consequently grains that are deposited in one time step and entrained in the following one do not contribute to the sediment cover, and so the model implicitly incorporates the effect of local sediment cover on grain deposition.

Model runs were completed with six different combinations of $p_{\text{i}}$ and $p_{\text{c}}$: 0.95/0.95, 0.95/0.75, 0.75/0.10, 0.75/0.30, 0.30/0.30, and 0.95/0.05. These combinations were selected to cover the range of relationships between relative sediment supply $Q_{\text{s}}^{*}$ and the exposed bed fraction $A^{*}$ observed by Hodge and Hoey (2012). For each pair of $p_{\text{i}}$ and $p_{\text{c}}$, model runs were completed with at least 20 different values of $Q_{\text{s}}^{*}$ in order to quantify the model behaviour.

Cover bed fraction and total mass on the bed produced by the model were converted using Eq. (3) into the new probabilistic framework (Fig. 4). The derivative was approximated by simple linear finite differences, which, in the case of run-away alluviation, resulted in a non-continuous curve due to large gradients. The exponential benchmark (Eq. 8) is also shown for comparison. The different model parameterizations produce results in which the probability of deposition on bedrock is both more and less likely than in the baseline case, with some runs showing both behaviours. Cases where the probability is more than the baseline case (i.e. grains are more likely to fall on uncovered areas) are associated with runs in which grains in clusters are relatively immobile. These runs are likely to be particularly affected by the smoothing algorithm that acts to move sediment from alluviated to bedrock areas. All model parameterizations predict greater bed exposure for a given normalized mass than is predicted by a linear cover relationship (Fig. 3b). Runs with relatively more immobile cluster grains have a lower exposed fraction for the same normalized mass. Runs with low values of $p_{\text{i}}$ and $p_{\text{c}}$ seem to lead to behaviour in which cover is more likely than in the exponential benchmark, while for high values, it is less likely. However, these are complex interactions and it is difficult to generalize the model behaviour.

Probability functions $P$ and cover function derived from data obtained from the model of Hodge and Hoey (2012). The grey dashed line shows the exponential benchmark behaviour. Grey shading depicts the area where the cover function cannot run due to conservation of mass. The legend gives values of the probabilities of entrainment $p_{\text{i}}$ and $p_{\text{c}}$ used for the runs (see text).

[Figure omitted. See PDF]

Cover development in time and space

Model derivation

Previous descriptions of the cover effect relate the exposed fraction of the bed to the relative sediment supply $Q_{\text{s}}^{*}$ (see Eqs. 1 and 2). In this section, we derive a model to clarify the relationship between the exposed fraction, $Q_{\text{s}}^{*}$, and $M_{\text{s}}$ and put it on a sound physical basis. To this end, the probabilistic formulation introduced previously is extended to allow the calculation of the temporal and spatial evolution of sediment cover in a stream. Here, we will derive the equations for the one-dimensional case (linear flume), but extensions to higher dimensions are possible in principle. The derivation is inspired by the erosion–deposition framework (e.g. Charru et al., 2004), with some necessary adaptions to make it suitable for channels with partial sediment cover (e.g. Turowski, 2009). In our system, we consider two separate mass reservoirs within a control volume. The first reservoir contains all particles in motion, the total mass per bed area of which is denoted by $M_{\text{m}}$, while the second reservoir contains all particles that are stationary on the bed, the total mass per bed area of which is denoted by $M_{\text{s}}$. The reservoirs exchange mass by entrainment and deposition; i.e. when a stationary particle is entrained it becomes mobile and when a mobile particle is deposited, it becomes stationary. In addition to Eq. (3), we then need three further equations: one to connect the rate of change in mobile mass to the sediment flux in the flume and one each to describe mass conservation in the two reservoirs. Instead of the common approach tracking the height of the sediment over a reference level, as is done in the classic mass conservation in fluvial systems, the Exner equation (e.g. Paola and Voller, 2005), we use the total sediment mass on the bed as a variable. Mobile sediment mass is supplied from upstream ($\mathrm{\Delta }$in), leaves in the downstream direction ($\mathrm{\Delta }$out), and can be exchanged between the stationary and the mobile mass reservoirs by entrainment ($E_{\text{tot}})$ and deposition ($D_{\text{tot}})$ (Fig. 5). The latter two parameters describe the exchange of particles between reservoirs; in the single-reservoir Exner equation these terms are not needed. It is clear that for the problem at hand the choice of total mass or volume as a variable to track the amount of sediment in the reach of interest is preferable to the height of the alluvial cover, since necessarily, when cover is patchy, the height of the alluvium varies across the bed.

Sediment dynamics at the bed are modelled by two reservoirs for stationary and mobile mass, which can exchange material by entrainment ($E_{\text{tot}})$ and deposition ($D_{\text{tot}})$. Sediment mass can be supplied from upstream ($\mathrm{\Delta }$in) and can leave in the downstream direction ($\mathrm{\Delta }$out).

[Figure omitted. See PDF]

The difference form of the mass balance for the mobile sediment is then given by (cf. Fig. 5) $${\displaystyle }\mathrm{\Delta }{M}_{\text{m}}=\left(\mathrm{\Delta },\mathrm{in},-,\mathrm{\Delta },\mathrm{out},+,E_{\text{tot}},-,D_{\text{tot}}\right)\mathrm{\Delta }t.$$

Here, $\mathrm{\Delta }{M}_{\text{m}}$ is the change in mobile sediment mass and $\mathrm{\Delta }t$ is a change in time. As the length of a time step is reduced to 0, a continuous version of Eq. (12) is obtained, which reads $${\displaystyle \frac{\partial M_{\text{m}}}{\partial t}}=-{\displaystyle \frac{\partial q_{\text{s}}}{\partial x}}+E-D.$$

Here, $x$ is the coordinate in the streamwise direction, $t$ the time, and $q_{\text{s}}$ the sediment mass transport rate per unit width, while $E$ is the mass entrainment rate per bed area and $D$ is the mass deposition rate per bed area. Similarly, in the mass balance for the stationary mass reservoir, the rate of change of the stationary sediment mass $M_{\text{s}}$ in time is the difference of the deposition rate $D$ and the entrainment rate $E$: $${\displaystyle \frac{\partial M_{\text{s}}}{\partial t}}=D-E.$$

It is useful to work with dimensionless variables by defining $t^{\ast }=t/T$ and $x^{*}=x/L$, where $T$ and $L$ are suitable time and length scales, respectively. The dimensionless mobile mass per bed area $M_{\text{m}}^{*}$ is equal to $M_{\text{m}}/M_{\mathrm{0}}$, and Eq. (13) becomes $${\displaystyle \frac{\partial M_{\text{m}}^{*}}{\partial t^{*}}}=-{\displaystyle \frac{\partial q_{\text{s}}^{*}}{\partial x^{*}}}+E^{*}-D^{*}.$$

Here, $${\displaystyle }{q}_{\text{s}}^{*}={\displaystyle \frac{T}{L{M}_{\mathrm{0}}}}{q}_{\text{s}}.$$

The dimensionless entrainment and deposition rates, $E^{*}$ and $D^{*}$, are equal to $TE/M_{\mathrm{0}}$ and $TD/M_{\mathrm{0}}$, respectively. Similarly, the balance for the stationary mass (Eq. 14) can be written as $${\displaystyle \frac{\partial M_{\text{s}}^{*}}{\partial t^{*}}}=D^{*}-E^{*}.$$

We also need sediment entrainment and deposition functions. The entrainment rate needs to be modulated by the availability of sediment on the bed. If $M_{\text{s}}^{*}$ is equal to 0, no material can be entrained. A plausible assumption is that the maximal entrainment rate, $E_{\text{max}}^{*}$, is equal to the transport capacity. $${\displaystyle }{E}_{\text{max}}^{*}=q_{\text{t}}^{*}$$

Here, $q_{\text{t}}^{*}$ is the dimensionless mass transport capacity, which is related to the transport capacity per unit width $q_{\text{t}}$ by a relation similar to Eq. (16). To first order, the rate of change in entrainment rate, d$E$, is proportional to the difference of $E_{\text{max}}$ and $E$ and to the rate of change in mass on the bed. $${\displaystyle }\text{d}{E}^{*}=\left(E_{\text{max}}^{*},-,E^{*}\right)\text{d}{M}_{\text{s}}^{*}=\left(q_{\text{t}}^{*},-,E^{*}\right)\text{d}{M}_{\text{s}}^{*}$$

Integrating, we obtain $${\displaystyle }{E}^{*}=E_{\text{max}}^{*}\left(\mathrm{1},-,e^{-M_{\text{s}}^{*}}\right)=\left(\mathrm{1},-,e^{-M_{\text{s}}^{*}}\right)q_{\text{t}}^{*}.$$

Here, we used the condition $E^{*}(M_{\text{s}}^{*}=\mathrm{0})=\mathrm{0}$ to fix the integration constant to $E_{\text{max}}^{*}$. As required, Eq. (20) approaches $E_{\text{max}}^{*}$ as $M_{\text{s}}^{*}$ goes to infinity and is equal to 0 when $M_{\text{s}}^{*}$ is equal to 0. Using a similar line of argument and by assuming the maximum deposition rate to be equal to $q_{\text{s}}^{*}$, we arrive at an equation for the deposition rate $D^{*}$. $${\displaystyle }{D}^{*}=\left(\mathrm{1},-,e^{-M_{\text{m}}^{*}}\right)q_{\text{s}}^{*}$$

When $M_{\text{m}}^{*}$ is small, then the amount that can be deposited is limited by $M_{\text{m}}^{*}$. If $M_{\text{m}}*$ is large, then deposition is limited by sediment supply. Substituting Eqs. (20) and (21) into Eq. (17), we obtain $$\begin{array}{cc}& {\displaystyle }& {\displaystyle \frac{\partial M_{\text{s}}^{*}\left(x^{*},,,t^{*}\right)}{\partial t^{*}}}=D^{*}-E^{*}={\displaystyle }& {\displaystyle }\left(\mathrm{1},-,e^{-M_{\text{m}}^{*}\left(x^{*},,,t^{*}\right)}\right)q_{\text{s}}^{*}\left(x^{*},,,t^{*}\right)-\left(\mathrm{1},-,e^{-M_{\text{s}}^{*}\left(x^{*},,,t^{*}\right)}\right)q_{\text{t}}^{*}\left(x^{*},,,t^{*}\right).\end{array}$$

Note that $q_{\text{s}}^{*}/q_{\text{t}}^{*}=Q_{\text{s}}^{*}$. The equation for the mobile mass (Eq. 14) becomes $$\begin{array}{cc}& {\displaystyle \frac{\partial M_{\text{m}}^{*}\left(x^{*},,,t^{*}\right)}{\partial t^{*}}}=& {\displaystyle }-{\displaystyle \frac{\partial q_{\text{s}}^{*}}{\partial x^{*}}}-\left(\mathrm{1},-,e^{-M_{\text{m}}^{*}\left(x^{*},,,t^{*}\right)}\right)q_{\text{s}}^{*}\left(x^{*},,,t^{*}\right){\displaystyle }& {\displaystyle }+\left(\mathrm{1},-,e^{-M_{\text{s}}^{*}\left(x^{*},,,t^{*}\right)}\right)q_{\text{t}}^{*}\left(x^{*},,,t^{*}\right).\end{array}$$

Finally, the sediment transport rate needs to be proportional to the mobile sediment mass times the downstream sediment speed $U$, and we can write $${\displaystyle }{q}_{\text{s}}^{*}\left(x^{*},,,t^{*}\right)=U^{*}\left(x^{*},,,t^{*}\right)M_{\text{m}}^{*}\left(x^{*},,,t^{*}\right).$$ Here $${\displaystyle }{U}^{*}={\displaystyle \frac{T}{L}}U.$$

After incorporating the original equation between $A^{*}$ and $M_{\text{s}}^{*}$ (Eq. 3), the system of four differential Eqs. (3), (22), (23) and (24) contains four unknowns: the downstream gradient in the sediment transport rate $\partial q_{\text{s}}^{*}/\partial x^{*}$, the exposed fraction of the bed $A^{*}$, the non-dimensional stationary mass $M_{\text{s}}^{*}$, and the non-dimensional mobile mass $M_{\text{m}}^{*}$, while the non-dimensional transport capacity $q_{\text{t}}^{*}$ and the non-dimensional downstream sediment speed $U^{*}$ are input variables and $P$ is a externally specified function. In addition, sediment input $q_{\text{s}}^{*}$ needs to be specified as an upstream boundary condition and initial values for the mobile mass $M_{\text{m}}^{*}$, and the stationary mass $M_{\text{s}}^{*}$ need to be specified everywhere.

Time-independent solution

In this chapter, we discuss the steady solution to the system of equations and thus clarify the relationship between cover, stationary sediment mass, sediment supply, and transport capacity. Setting the time derivatives to 0, we obtain a time-independent solution, which links the exposed area directly to the ratio of sediment transport rate to transport capacity. From Eq. (23) it follows that in this case, the entrainment rate is equal to the deposition rate, and we obtain $${\displaystyle }\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{M_{\text{m}}^{*}}}\right)\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}=\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}}\right)q_{\text{t}}^{*}.$$

Here, the bar over the variables denotes their steady-state value. Substituting Eq. (24) to eliminate $\stackrel{\mathrm{‾}}{M_{\text{m}}^{*}}$ and solving for $\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}$ gives $$\begin{array}{cc}& {\displaystyle }\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}& {\displaystyle }=-\ln \left\{\mathrm{1},-,\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}}\right),{\displaystyle \frac{\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}}{q_{\text{t}}^{*}}}\right\}{\displaystyle }& {\displaystyle }=-\ln \left\{\mathrm{1},-,\left(\mathrm{1},-,e^{-\frac{q_{\text{t}}^{*}}{U^{*}}\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}}\right),\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}\right\}.\end{array}$$

Note that we assume here that sediment cover is only dependent on the stationary sediment mass on the bed, and we thus neglect grain–grain interactions known as the dynamic cover (Turowski et al., 2007). In analogy to Eq. (24), we can write $${\displaystyle }{q}_{\text{t}}^{*}=U^{*}{M}_{\mathrm{0}}^{*}.$$

Here, $M_{\mathrm{0}}^{*}$ is a characteristic dimensionless mass that depends on hydraulics and therefore implicitly on transport capacity (which should not be confused with the minimum mass necessary to fully cover the bed $M_{\mathrm{0}})$. When sediment transport rate equals transport capacity, then $M_{\mathrm{0}}^{*}$ is equal to the mobile mass of sediment normalized by the reference mass $M_{\mathrm{0}}$. It can be viewed as a proxy for the transport capacity and is a convenient parameter to simplify the equations. The mobile mass can then, in general, be written as follows (cf. Turowski et al., 2007), remembering that the relative sediment supply $Q_{\text{s}}^{*}=\mathrm{1}$ when supply is equal to capacity: $${\displaystyle }{M}_{\text{m}}^{*}=M_{\mathrm{0}}^{*}{Q}_{\text{s}}^{*}.$$

If we use the exponential cover function (Eq. 8) with Eqs. (27), (28), and (29), we obtain $$\begin{array}{cc}{\displaystyle }\stackrel{\mathrm{‾}}{A^{*}}& {\displaystyle }=\mathrm{1}-\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}}\right){\displaystyle \frac{\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}}{q_{\text{t}}^{*}}}=\mathrm{1}-\left(\mathrm{1},-,e^{-\frac{q_{\text{t}}^{*}}{U^{*}}\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}}\right)\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}\\ & {\displaystyle }& {\displaystyle }=\mathrm{1}-\left(\mathrm{1},-,e^{-M_{\mathrm{0}}^{*}\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}}\right)\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}.\end{array}$$

Similarly, equations can be found for the other analytical solutions of the cover function. For the linear case (Eq. 6), we obtain $${\displaystyle }\stackrel{\mathrm{‾}}{A^{*}}=\mathrm{1}+\ln \left\{\mathrm{1},-,\left(\mathrm{1},-,e^{-M_{\mathrm{0}}^{*}\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}}\right),\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}\right\}.$$

For the power law case (Eq. 9), we obtain $${\displaystyle }\stackrel{\mathrm{‾}}{A^{*}}={\left[\mathrm{1},+,\left(\mathrm{1},-,\mathit{\alpha }\right),\ln ,\left\{\mathrm{1},-,\left(\mathrm{1},-,e^{-M_{\mathrm{0}}^{*}\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}}\right),\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}\right\}\right]}^{\frac{\mathrm{1}}{\mathrm{1}-\mathit{\alpha }}}.$$

The exponential cover function essentially leads to a combined linear and exponential relation between $\stackrel{\mathrm{‾}}{A^{*}}$ and $\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}$. Instead of a linear decline as in the original linear cover model (Eq. 1) or a concave-up relationship as in the original exponential model (Eq. 2), the function is convex-up for all solutions (Fig. 6). Adjusting $M_{\mathrm{0}}^{*}$ shifts the lines: decreasing $M_{\mathrm{0}}^{*}$ leads to a delayed onset of cover and vice versa. The former result arises because a lower $M_{\mathrm{0}}^{\ast }$ means that the sediment flux is conveyed through a smaller mass moving at a higher velocity. The original linear cover function (Eq. 1) can be recovered from the exponential model with a high value of $M_{\mathrm{0}}^{*}$, since the exponential term quickly becomes negligible with increasing $\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}$ and the linear term dominates (Fig. 6c). Note that for the linear (Eq. 5) and the power law cases (Eq. 9), high values of $M_{\mathrm{0}}^{*}$ may give $\stackrel{\mathrm{‾}}{A^{*}}=\mathrm{0}$ for $\stackrel{\mathrm{‾}}{Q_{\text{s}}^{*}}<\mathrm{1}$ (Fig. 6b, d), which is consistent with the concept of runaway alluviation. Using the Beta distribution to describe $P$, a numerical solution is necessary, but a wide range of steady-state cover functions can be obtained (Fig. 7). By varying the value of $M_{\mathrm{0}}^{*}$, an even wider range of behaviours can be obtained.

Analytical solutions at steady state for the exposed fraction of the bed ($A^{*}$) as a function of relative sediment supply ($Q^{*}$, cf. Fig. 2). (a) Comparison of the different solutions, keeping $M_{\mathrm{0}}^{*}$ constant at 1. (b) Varying $M_{\mathrm{0}}^{*}$ for the linear case (Eq. 31). (c) Varying $M_{\mathrm{0}}^{*}$ for the exponential case (Eq. 30). (d) Varying $M_{\mathrm{0}}^{*}$ for the power law case with $\mathit{\alpha }=\mathrm{0.1}$ (Eq. 32).

[Figure omitted. See PDF]

Steady-state solutions using the Beta distribution to parameterize $P$ (Eq. 10) for a range of parameters $a$ and $b$ and using $M_{\mathrm{0}}^{*}=\mathrm{1}$ (cf. Fig. 3). The solutions were obtained by iterating the equations to a steady state, using initial conditions of $A^{*}=\mathrm{1}$ and $M_{\text{m}}^{*}=M_{\text{s}}^{*}=\mathrm{0}$.

[Figure omitted. See PDF]

The previous analysis shows that steady-state cover is controlled by the characteristic dimensionless mass $M_{\mathrm{0}}^{*}$, which is equal to the ratio of dimensionless transport capacity and particle speed (Eq. 28). In the following, we relate $M_{\mathrm{0}}^{*}$ to hydraulic variables and argue that it is, in general, not a constant. Converting $M_{\mathrm{0}}^{*}$ to dimensional variables, we can write $${\displaystyle }{M}_{\mathrm{0}}^{*}={\displaystyle \frac{q_{\text{t}}^{*}}{U^{*}}}={\displaystyle \frac{q_{\text{t}}}{M_{\mathrm{0}}U}}.$$

The minimum mass necessary to completely cover the bed per unit area, $M_{\mathrm{0}}$, can be estimated assuming a single layer of closely packed spherical grains residing on the bed (cf. Turowski, 2009), giving $${\displaystyle }{M}_{\mathrm{0}}={\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\text{s}}{D}_{\mathrm{50}}}{\mathrm{3}\sqrt{\mathrm{3}}}}.$$

Here, ${\mathit{\rho }}_{\text{s}}$ is the sediment density and $D_{\mathrm{50}}$ is the median grain size. We use equations derived by Fernandez-Luque and van Beek (1976) from flume experiments that describe transport capacity and particle speed as a function of bed shear stress (see also Lajeunesse et al., 2010, and Meyer-Peter and Mueller, 1948, for similar equations): $${\displaystyle }{q}_{\text{t}}=\mathrm{5.7}{\displaystyle \frac{{\mathit{\rho }}_{\text{s}}\mathit{\rho }}{\left({\mathit{\rho }}_{\text{s}},-,\mathit{\rho }\right)g}}{\left({\displaystyle \frac{\mathit{\tau }}{\mathit{\rho }}},-,{\displaystyle \frac{{\mathit{\tau }}_{\text{c}}}{\mathit{\rho }}}\right)}^{\mathrm{3}/\mathrm{2}},$$ $${\displaystyle }U=\mathrm{11.5}\left({\left({\displaystyle \frac{\mathit{\tau }}{\mathit{\rho }}}\right)}^{\mathrm{1}/\mathrm{2}},-,\mathrm{0.7},{\left({\displaystyle \frac{{\mathit{\tau }}_{\text{c}}}{\mathit{\rho }}}\right)}^{\mathrm{1}/\mathrm{2}}\right).$$

Here, ${\mathit{\tau }}_{\text{c}}$ is the critical bed shear stress for the onset of bedload motion, $g$ is the acceleration due to gravity, and $\mathit{\rho }$ is the water density. Combining Eqs. (34), (35), and (36) to get an equation for $M_{\mathrm{0}}^{*}$ gives $${\displaystyle }{M}_{\mathrm{0}}^{*}={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{{\left(\mathit{\theta },-,{\mathit{\theta }}_{\text{c}}\right)}^{\mathrm{3}/\mathrm{2}}}{{\mathit{\theta }}^{\mathrm{1}/\mathrm{2}}-\mathrm{0.7}{\mathit{\theta }}_{\text{c}}^{\mathrm{1}/\mathrm{2}}}}={\displaystyle \frac{\mathrm{3}\sqrt{\mathrm{3}}{\mathit{\theta }}_{\text{c}}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{{\left(\mathit{\theta },/,{\mathit{\theta }}_{\text{c}},-,\mathrm{1}\right)}^{\mathrm{3}/\mathrm{2}}}{{\left(\mathit{\theta },/,{\mathit{\theta }}_{\text{c}}\right)}^{\mathrm{1}/\mathrm{2}}-\mathrm{0.7}}}.$$

Here, the Shields stress $\mathit{\theta }=\mathit{\tau }/({\mathit{\rho }}_{\text{s}}-\mathit{\rho })g{D}_{\mathrm{50}}$, and ${\mathit{\theta }}_{\text{c}}$ is the corresponding critical Shields stress, and we approximated $\mathrm{5.7}/\mathrm{11.5}=\mathrm{0.496}$ with $\mathrm{1}/\mathrm{2}$ (compare to Eqs. 35, 36). At high $\mathit{\theta }$, when the threshold can be neglected, Eq. (37) reduces to a linear relationship between $M_{\mathrm{0}}^{*}$ and $\mathit{\theta }$. Near the threshold, $M_{\mathrm{0}}^{*}$ is shifted to lower values as ${\mathit{\theta }}_{\text{c}}$ increases (Fig. 8). The systematic variation of $U^{*}$ with the hydraulic driving conditions (Eq. 36) implies that the cover function evolves differently in response to changes in sediment supply and transport capacity. For a first impression, by comparing Eqs. (35) and (36), we assume that particle speed scales with transport capacity raised to the power of one-third (Fig. 9).

The characteristic dimensionless mass $M_{\mathrm{0}}^{*}$ depicted as a function of (a) the Shields stress and (b) the ratio of Shields stress to critical Shields stress (Eq. 37).

[Figure omitted. See PDF]

Variation of the exposed bed fraction as a function of transport capacity, assuming that particle speed scales with transport capacity to the power of one-third.

[Figure omitted. See PDF]

Temporal evolution of cover within a reach

To calculate the temporal evolution of cover on the bed within a single reach, we solved Eqs. (3), (22), (23), and (24) numerically for a section of the bed with homogenous conditions using a simple linear finite difference scheme. In this case sediment input is a boundary condition, while sediment output, mobile and stationary sediment mass, and the fraction of the exposed bed are output variables. In general, a change in sediment supply leads to a gradual adjustment of the output variables towards a new steady state (Fig. 10). It is desirable to obtain expressions for the response time of the system to external perturbation, such as a change in sediment supply or hydraulic conditions. Such a response time could then be compared to the timescales of changes in boundary conditions. For example, during a flood event, both transport capacity and sediment supply change over time. If these changes are slow in comparison to the response time of cover, the bed cover state can essentially keep up with the imposed changes at all times, and therefore steady-state equations (Sect. 3.2) can be used to calculate its evolution. In contrast, if the imposed change is rapid in comparison to the response time, cover may lag behind, and an approach that resolves cover as a dynamic variable is necessary. This may, for example, be important when studying the erosional behaviour of channels in response to floods (see Lague, 2010; Turowski et al., 2013). Unfortunately, a general analytical solution is not possible, but results can be obtained for special cases. We first derive analytical solutions for the response time for a reach without upstream sediment supply and for a system responding to small perturbations in sediment supply or transport capacity (Sect. 3.3.1) and discuss the system behaviour (Sect. 3.3.2). Finally, we apply the concepts to data from a flood in a natural river and demonstrate that, for this specific case, because of the response times, the steady-state relations do not capture cover behaviour.

System timescales

First, consider a reach without upstream sediment supply; i.e. $q_{\text{s}}^{*}=\mathrm{0}$. Such a situation is rare in nature but could be easily created in flume experiments as a model test. Then, the time derivative of stationary mass is given by $${\displaystyle \frac{\partial M_{\text{s}}^{*}}{\partial t^{*}}}=-\left(\mathrm{1},-,e^{-M_{\text{s}}^{*}}\right)q_{\text{t}}^{*}.$$

Using the exponential cover model (Eq. 8), we obtain $${\displaystyle \frac{\mathrm{1}}{A^{*}\left(\mathrm{1},-,A^{*}\right)}}{\displaystyle \frac{\partial A^{*}}{\partial t^{*}}}=q_{\text{t}}^{*}.$$

Equation (39) is separable and can be integrated to obtain $${\displaystyle }\ln \left(A^{*}\right)-\ln \left(\mathrm{1},-,A^{*}\right)=t^{*}{q}_{\text{t}}^{*}+C.$$

Letting $A^{*}(t^{*}=$0) $=A_{\mathrm{0}}^{*}$, where $A_{\mathrm{0}}^{*}$ is the initial cover, the final equation is in the form of a sigmoidal-type function: $${\displaystyle }{A}^{*}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\left(\frac{\mathrm{1}-A_{\mathrm{0}}^{*}}{A_{\mathrm{0}}^{*}}\right)e^{-t^{*}{q}_{\text{t}}^{*}}}}.$$

By making the parameters in the exponent on the right-hand side of Eq. (42) dimensional, we get $${\displaystyle }{t}^{*}{q}_{\text{t}}^{*}={\displaystyle \frac{t}{T}}{\displaystyle \frac{T}{L{M}_{\mathrm{0}}}}{q}_{\text{t}}={\displaystyle \frac{t{q}_{\text{t}}}{L{M}_{\mathrm{0}}}},$$

which allows a characteristic system timescale $T_{\text{E}}$ to be defined as $${\displaystyle }{T}_{\text{E}}={\displaystyle \frac{L{M}_{\mathrm{0}}}{q_{\text{t}}}}.$$

Since this timescale is dependent on the transport capacity $q_{\text{t}}$, we can view it as a timescale associated with the entrainment of sediment from the bed (cf. Eq. 20) – hence the subscript “E” on $T_{\text{E}}$. From Eq. (41), the exposed bed fraction evolves in an asymptotic fashion towards equilibrium (Fig. 11). We can expect that there are other characteristic timescales for the system, for example associated with sediment deposition or downstream sediment evacuation.

We can make some further progress and define a more general system timescale by performing a perturbation analysis (Appendix A1). For small perturbations in either $q_{\text{s}}^{*}$ or $q_{\text{t}}^{*}$, we obtain an exponential term describing the transient evolution, which allows the definition of a system timescale $T_{\text{S}}$ $${\displaystyle }\exp \left\{-,\left(\stackrel{\mathrm{‾}}{q_{\text{t}}^{*}},-,\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}\stackrel{\mathrm{‾}}{U^{*}}}\right),\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}\right),t^{*}\right\}=e^{-\frac{t}{T_{\text{S}}}}.$$

Here, exp denotes the natural exponential function. The characteristic system timescale can then be written as $${\displaystyle }{T}_{\text{S}}={\displaystyle \frac{L{M}_{\mathrm{0}}}{\stackrel{\mathrm{‾}}{q_{\text{t}}}\left(\mathrm{1},-,\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}\stackrel{\mathrm{‾}}{U^{*}}}\right),\frac{\stackrel{\mathrm{‾}}{q_{\text{s}}}}{\stackrel{\mathrm{‾}}{q_{\text{t}}}}\right)}}={\displaystyle \frac{L{M}_{\mathrm{0}}}{\stackrel{\mathrm{‾}}{q_{\text{t}}}}}{e}^{\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}}.$$ Note that for $q_{\text{s}}^{*}=\mathrm{0}$, Eq. (45) reduces to Eq. (43), as would be expected. Since $\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}$ is directly related to steady-state bed exposure $\stackrel{\mathrm{‾}}{A^{*}}$, we can rewrite the equation, for example by assuming the exponential cover function (Eq. 8), as $${\displaystyle }{T}_{\text{S}}={\displaystyle \frac{L{M}_{\mathrm{0}}}{\stackrel{\mathrm{‾}}{q_{\text{t}}}\stackrel{\mathrm{‾}}{A^{*}}}}.$$

Since bed cover is more easily measurable than the mass on the bed, Eq. (46) can help to estimate system timescales in the field. Further, $\stackrel{\mathrm{‾}}{A^{*}}$ varies between 0 and 1, which allows the estimation of a minimum system time using Eq. (43). As $\stackrel{\mathrm{‾}}{A^{*}}$ approaches 0, the system timescale diverges.

To illustrate these additional dependencies, we have used numerical solutions of Eqs. (3), (22), (23), and (24) to calculate the time needed to reach 99.9 % of total adjustment after a step change in transport stage (chosen due to the asymptotic behaviour of the system), analysed across a plausible range of particle speeds $U$ (Fig. 12). Response time decreases as particle speed increases. This reflects elevated downstream evacuation for higher particles speeds, resulting in a smaller mobile particle mass and thus higher entrainment and lower deposition rates. Response time also increases with increasing relative sediment supply $Q_{\text{s}}^{*}$. As the runs start with zero sediment cover and the extent of cover increases with $Q_{\text{s}}^{*}$, at higher $Q_{\text{s}}^{*}$the adjusted cover takes longer to develop.

Temporal evolution of cover for the simple case of a control box with sediment through-flux, based on Eqs. (3), (22), (23), and (24). Relative sediment supply (supply normalized by transport capacity) was specified to 0.25 and increased to 1 at $t^{*}=\mathrm{5}$. The response of sediment output, mobile and stationary sediment mass, and the exposed bed fraction was calculated. Here, we used the exponential function for $P$ (Eq. 8) and $M_{\mathrm{0}}^{*}=U^{*}=\mathrm{1}$. The initial values were $A^{*}=\mathrm{1}$ and $M_{\text{m}}^{*}=M_{\text{s}}^{*}=\mathrm{0}$.

[Figure omitted. See PDF]

Evolution of the exposed bed fraction (removal of sediment cover) over time starting with different initial values of bed exposure, for the special case of no sediment supply; i.e. $q_{\text{s}}^{*}=\mathrm{0}$ (Eq. 41) and $q_{\text{t}}^{*}=\mathrm{1}$.

[Figure omitted. See PDF]

Dimensionless time to reach 99.9 % of the total adjustment in the exposed area as a function of (a) transport stage and (b) particle speed. All simulation were started with $A^{*}=\mathrm{1}$ and $M_{\text{m}}^{*}=M_{\text{s}}^{*}=\mathrm{0}$.

[Figure omitted. See PDF]

Phase shift and gain in response to a cyclic perturbation

The perturbation analysis (Appendix A) gives some insight into the response of cover to cyclic sinusoidal perturbations. Let sediment supply be perturbed in a cyclic way described by an equation of the form $${\displaystyle }{q}_{\text{s}}^{*}=\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}+\mathit{\delta }{q}_{\text{s}}^{*}=\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}+\text{d}\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{p}}\right).$$

Here, the overbar denotes the temporal average, $\mathit{\delta }{q}_{\text{s}}^{*}$ is the time-dependent perturbation, $d$ is the amplitude of the perturbation, and $p$ its period. A similar perturbation can be applied to the transport capacity (see Appendix A). The reaction of the stationary mass and therefore cover can then also be described by sinusoidal function of the form (Appendix A) $${\displaystyle }\mathit{\delta }{M}_{\text{s}}^{*}=G\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{p}},+,\mathit{\phi }\right).$$

Here, $\mathit{\delta }{M}_{\text{s}}^{*}$ is the perturbation of the stationary sediment mass around the temporal average, $G$ is known as the gain, describing the amplitude response, and $\mathit{\phi }$ is the phase shift. If the gain is large, stationary mass reacts strongly to the perturbation; if it is small, the forcing does not leave a signal. The phase shift is negative if the response lags behind the forcing and positive if it leads. The phase shift can be written as $${\displaystyle }\mathit{\phi }={\tan }^{-\mathrm{1}}\left(-,\mathrm{2},\mathit{\pi },{\displaystyle \frac{T_{\text{S}}}{p}}\right).$$

The gain can be written as $${\displaystyle }G={\displaystyle \frac{p}{T_{\text{S}}}}{\displaystyle \frac{Kd}{\sqrt{{\left(\frac{p}{T_{\text{S}}}\right)}^{\mathrm{2}}+\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}}}}.$$

Here, $d$ is the amplitude of the perturbation, and $K$ is a function of the time-averaged values of $q_{\text{s}}$, $q_{\text{t}}$, and $U$ and differs for perturbations in transport capacity and sediment supply (see Appendix A). Thus, the system behaviour can be interpreted as a function of the ratio of the period of perturbation $p$ and the system timescale $T_{\text{s}}$. The period $p$ is large if the forcing parameter, i.e. discharge or sediment supply, varies slowly and small when it varies quickly. According to Eq. (49), the phase shift is equal to $-\mathit{\pi }/\mathrm{2}$ for low values of $p/T_{\text{s}}$ (quickly varying forcing parameter), implying a substantial lag in the adjustment of cover. The phase shift tends to 0 as $p/T_{\text{s}}$ tends to infinity (Fig. 13). The gain varies approximately linearly with $p/T_{\text{s}}$ for small $p/T_{\text{s}}$ (quickly varying forcing parameter), while it is approximately constant at a value of $Kd$ for large $p/T_{\text{s}}$ (slowly varying forcing parameter) (Eq. 50). Thus, if the forcing parameter varies slowly, cover adjustment keeps up at all times.

Phase shift (Eq. 49) and gain (Eq. 50) as a function of the ratio of the period of perturbation $p$ and the system timescale $T_{\text{s}}$. For the calculation, the constant factor in the gain ($Kd$) was set equal to 1.

[Figure omitted. See PDF]

A flood at the Erlenbach

To illustrate the magnitude of the timescales using real data, we use a flood data set from the Erlenbach, a sediment transport observatory in the Swiss Prealps (e.g. Beer et al., 2015). There, near a discharge gauge, bedload transport rates are measured at 1 min resolution using the Swiss Plate Geophone System, a highly developed and fully calibrated surrogate bedload measuring system (e.g. Rickenmann et al., 2012; Wyss et al., 2016). We use data from a flood on 20 June 2007 (Turowski et al., 2009) with the highest peak discharge that has so far been observed at the Erlenbach. The meteorological conditions that triggered this flood and its geomorphic effects have been described in detail elsewhere (Molnar et al., 2010; Turowski et al., 2009, 2013). The Erlenbach does not have a bedrock bed in the sense that bedrock is exposed in the channel bed; however, the data provide a realistic natural time series of discharge and bedload transport over the course of a single event. Rather than predicting bed cover evolution for a natural system, for which we do not currently have data for validation, we use the Erlenbach data to illustrate possible cover behaviour during a fictitious event with different initial sediment cover extents, using natural data to provide realistic boundary conditions.

Using a median grain size of 80 mm, a sediment density of 2650 kg m${}^{-\mathrm{3}}$, and a reach length of 50 m, we obtained $M_{\mathrm{0}}=\mathrm{128}$ kg m${}^{-\mathrm{2}}$. We calculated transport capacity using the equation of Fernandez Luque and van Beek (1976). However, it is known that this and similar equations strongly overestimate measured transport rates in streams such as the Erlenbach (e.g. Nitsche et al., 2011). Consequently, we rescaled by setting the ratio of bedload supply to capacity to 1 at the highest discharge. The exposed fraction was then calculated iteratively assuming $P=A^{*}$ (i.e. the exponential cover formulation, Eq. 8). In a real flood event, water discharge and sediment supply obviously do not follow a small cyclic perturbation (Fig. 13). But we can tentatively relate the observations to the theory by assuming that at each time step, the change in sediment supply can be represented by the commencement of a sinusoidal perturbation with varying period. To estimate the effective period $p$, one needs to take the derivatives of Eq. (47). $${\displaystyle \frac{\text{d}{q}_{\text{s}}^{*}}{\text{d}t}}={\displaystyle \frac{\text{d}\mathit{\delta }{q}_{\text{s}}^{*}}{\text{d}t}}={\displaystyle \frac{\mathrm{2}\mathit{\pi }d}{p}}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{p}}\right)$$

Setting $t=\mathrm{0}$ for the time of interest, we can relate $p$ to the local gradient in bedload supply, which can be measured from the data. $${\displaystyle \frac{\mathrm{2}\mathit{\pi }d}{p}}={\displaystyle \frac{\mathrm{\Delta }{q}_{\text{s}}^{*}}{\mathrm{\Delta }t}}$$

Assuming that all change in the response time is due to changes in the period (i.e. assuming a constant amplitude, $d=\mathrm{1}$), we can obtain a conservative estimate of the range over which $p$ varies over the course of an event. $${\displaystyle }p=\mathrm{2}\mathit{\pi }{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }{q}_{\text{s}}^{*}}}$$

In the exemplary event, the evolution and final value of bed cover depends strongly on its initial value (Fig. 14), indicating that the adjustment is incomplete. The system timescale is generally larger than 1000 s and is inversely related to discharge via the dependence on transport capacity. The $p/T_{\text{s}}$ ratio varies around 1, with low values at the beginning of the flood and large values in the waning hydrograph. Both the high values of the system timescale and the smooth evolution of bed cover over the course of the flood imply that cover development cannot keep up with the variation in the forcing characteristics. This dynamic adjustment of cover, which can lag forcing processes, may thus play an important role in the dynamics of bedrock channels and probably needs to be taken into account in modelling.

Calculated evolution of cover during the largest event observed at the Erlenbach on 20 June 2007 (Turowski et al., 2009). Bedload transport rates were measured with the Swiss Plate geophone sensors calibrated with direct bedload samples (Rickenmann et al., 2012). The final fraction of exposed bedrock is strongly dependent on its initial value.

[Figure omitted. See PDF]

Discussion

Model formulation

In principle, the framework for the cover effect presented here allows the formulation of a general model for bedrock channel morphodynamics without the restrictions of previous models (e.g. Nelson and Seminara, 2011; Zhang et al., 2015). To achieve this, the dependency of $P$ on various control parameters needs to be specified. In general, $P$ should be controlled by local topography, grain size and shape, hydraulic forcing, and the amount of sediment already residing on the bed. Furthermore, the shape of the $P$ function should also be affected by feedbacks between these properties, such as the development of sediment cover altering the local roughness and hence altering hydraulics and local transport capacity (Inoue et al., 2014; Johnson, 2014). In the treatment presented here, we have explicitly accounted only for the impact of the amount of sediment already residing on the bed. However, all of the mentioned effects can be included implicitly by an appropriate choice of $P$. The exact relationships between, say, bed topography and $P$ need to be mapped out experimentally (e.g. Inoue et al., 2014), with theoretical approaches also providing some direction (cf. Johnson, 2014; Zhang et al., 2015). Currently available experimental results (Chatanantavet and Parker, 2008; Finnegan et al., 2007; Hodge and Hoey, 2016; Inoue et al., 2014; Johnson and Whipple, 2007) cover only a small range of the possible parameter space, and, in general, not all necessary parameters to constrain $P$ were reported. Specifically the stationary mass of sediment residing on the bed is usually not reported and can be difficult to determine experimentally but is necessary to determine $P$. Nevertheless, depending on the choice of $P$, our model can yield a wide range of cover functions that encompasses reported functions both from numerical modelling (e.g. Aubert et al., 2016; Hodge and Hoey, 2012; Johnson, 2014) and experiments (Chatanantavet and Parker, 2008; Inoue et al., 2014; Sklar and Dietrich, 2001) (see Figs. 4 and 5).

The dynamic model put forward here is a minimum first-order formulation, and there are some obvious future alterations. We only take account of the static cover effect caused by immobile sediment on the bed. The dynamic cover effect, which arises when moving grains interact at high sediment concentration and thus reduce the number of impacts on the bed (Turowski et al., 2007), could in principle be included in the formulation but would necessitate a second probability function specifically to describe this dynamic cover. It would also be possible to use different $P$ functions for entrainment and deposition, thus introducing hysteresis into cover development. Such hysteresis has been observed in experiments in which the equilibrium sediment cover was a function of the initial extent of sediment cover (Chatanantavet and Parker, 2008; Hodge and Hoey, 2012). Whether such alterations are necessary is best established with targeted laboratory experiments.

Comparison to previous modelling frameworks

We will briefly outline in this section the main differences to previous formulations of cover dynamics in bedrock channels. Thus, the novel aspects of our formulation and the respective advantages and disadvantages will become clear.

Aubert et al. (2015) coupled the movement of spherical particles to the simulation of a turbulent fluid and investigated how cover depends on transport capacity and supply. Similar to what is predicted by our analytical formulation, they found a range of cover function for various model set-ups, including linear and convex-up relationships (compare the results in Fig. 6 to their Fig. 15). Aubert et al. (2015) presented the most detailed physical simulations of bed cover formation so far, and the correspondence between the predictions is encouraging.

Nelson and Seminara (2011, 2012) formulated a morphodynamic model for bedrock channels. They based their formulation on sediment concentration, which is in principle similar to our formulation based on mass. However, Nelson and Seminara (2011, 2012) did not distinguish between mobile and stationary sediment and linked local transport directly to sediment concentration. Further, Nelson and Seminara (2011, 2012) assumed a direct correspondence between sediment concentration and degree of cover, which is equivalent to the linear cover function (Eq. 6). In this case, it is assumed that grains are always deposited on uncovered bed, and the different possible distributions of particles within a grid node are not taken into account. Practically, this implies that the grid size needs to be of the order of the grain size because, strictly, the assumption is only valid if a single grain can cover an entire grid node (cf. Fig. 1). Although different in various details, Inoue et al. (2016) have used essentially the same approach as Nelson and Seminar (2011, 2012) to link bedload concentration, transport, and bed cover. Both of these models allow the 2-D modelling of bedrock channel morphology. Although we have not fully developed such a model in the present paper, our model framework could easily be extended to 2-D problems.

Inoue et al. (2014) formulated a 1-D model for cover dynamics and bedrock erosion. There, they distinguish between stationary and mobile sediment using an Exner equation to capture sediment mass conservation. The degree of bed cover is related to transport rates and sediment mass via a saturation volume, which is related to our characteristic mass $M_{\mathrm{0}}^{*}$ (see Sect. 3.2). A key difference between Inoue et al.'s (2014) model and the ones presented here lies in the sediment mass conservation equations (Eqs. 13 and 14), in which we explicitly take account of both entrainment and deposition. In addition, with the function $P$, describing the relationship between deposited mass and degree of cover, we provide a more flexible framework for complex simulations where the bed needs to be discretized (e.g. 2-D models or reach-scale formulations).

Zhang et al. (2015) formulated a bed cover model specifically for beds with macro-roughness. There, deposited sediment always fills topographic lows from their deepest positions, such that there is a reach-uniform sediment level. While the model provides a fundamentally different approach to what is suggested here, its applicability is limited to very rough beds, and the assumption of a sediment elevation that is independent of the position on the bed seems physically unrealistic. In principle, the probabilistic framework presented here should be able to deal with macro-rough beds, by making the $P$ function (Eq. 3) explicitly dependent on roughness, and thus it allows a more general treatment of the problem of bed cover.

In this paper, we focused on the dynamics of bed cover rather than on the modelling of the dynamics of entire channels. The probabilistic formulation using the parameter $P$ provides a flexible framework to connect the sediment mass residing on the bed with the exposed bedrock fraction. This particular element has not been treated in any of the previous models and could be easily implemented in other approaches dealing with sediment fluxes along and across the stream and the interaction with erosion and, over long timescales, channel morphology. However, it is as yet unclear how flow hydraulics, sediment properties, and other conditions affect $P$, and this should be investigated in targeted laboratory experiments.

Further implications

Based on field data interpretation, Phillips and Jerolmack (2016) argued that bedrock rivers adjust such that, similar to alluvial channels, medium-sized floods are most effective in transporting sediment and that channel geometry therefore can quickly adjust their transport capacity to the applied load and therefore achieve grade (cf. Mackin, 1948). They conclude that bedrock channels can adjust their morphologic parameters (channel width and shape) quickly in response to changing boundary conditions. In contrast, our model suggests that instead bed cover can be adjusted to achieve grade. In steady state, time derivatives need to be equal to 0. Thus, entrainment equals deposition (Eq. 14), implying that the downstream gradient in sediment transport rate is equal to 0 (Eq. 13). When sediment supply or transport capacity change, the exposed bedrock fraction can adjust to achieve a new steady state, and a change in the channel geometry is unnecessary. These changes in sediment cover can occur far more rapidly than changes in width and cross-sectional shape (compare to Eq. 46). Whether a steady state is achieved depends on the relative magnitude of the timescales of perturbation and cover adjustment (see Sect. 3). Our results imply that bedrock channels have two distinct timescales to adjust to changing boundary conditions to achieve grade. Over short times, bed cover is adjusted. This can occur rapidly. Over long timescales, channel width, cross-sectional shape, and slope are adjusted.

Conclusions

The probabilistic view put forward in this paper offers a framework into which diverse data on bed cover, whether obtained from field studies, laboratory experiments, or numerical modelling, can be easily converted to be meaningfully compared. The conversion requires knowledge of the mass of sediment on the bed and the evolution of exposed fraction of the bed. Within the framework, individual data sets can be compared to the exponential benchmark and linear limit cases, enabling physical interpretation. Furthermore, the formulation allows the general dynamic sub-grid modelling of bed cover. Depending on the choice of $P$, the model yields a wide range of possible cover functions. Which of these functions are appropriate for natural rivers and how they vary with factors including topography needs to mapped out experimentally.

It needs to be noted here that the precise formulation of the entrainment and deposition functions also affects steady-state cover relations. When calibrating $P$ on data, it cannot always be decided whether a specific deviation from the benchmark case results from varying entrainment and deposition processes or from changes in the probability function driven for example by variations in roughness. For the prediction of the steady-state cover relations and for the comparison of data sets, this should not matter, but the dynamic evolution of cover could be strongly affected.

The system timescale for cover adjustment is inversely related to transport capacity. This timescale can be long and in many realistic situations, cover cannot instantaneously adjust to changes in the forcing conditions. Thus, dynamic cover adjustment needs to be taken into account when modelling the long-term evolution of bedrock channels.

Our model formulation implies that bedrock channels adjust bed cover to achieve grade. Therefore, bedrock channel evolution is driven by two optimization principles. On short timescales, bed cover adjusts to match the sediment output of a reach to its input. Over long timescales, width and slope of the channel evolve to match long-term incision rate to tectonic uplift or base-level lowering rates.

The model data used in Sect. 2.2 (Fig. 3) can be obtained by contacting the authors. The Erlenbach data used in Sect. 3.3.3 (Fig. 14) belong to the Mountain Hydrology and Mass Movements Group at the Swiss Federal Research Institute for Forest Snow and Landscape Research WSL and have been used with permission.

Perturbation analysis

Here, we derive the effect of a small sinusoidal perturbation of the driving variables, namely sediment supply $q_{\text{s}}^{*}$ and transport capacity $q_{\text{t}}^{*}$, on cover development. The perturbation of the driving variables can be written as $$\begin{array}{}& {\displaystyle }& {\displaystyle }{q}_{\text{s}}^{*}=\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}+\mathit{\delta }{q}_{\text{s}}^{*},& {\displaystyle }& {\displaystyle }{q}_{\text{t}}^{*}=\stackrel{\mathrm{‾}}{q_{\text{t}}^{*}}+\mathit{\delta }{q}_{\text{t}}^{*}.\end{array}$$

Here, the bar denotes the average of the quantity at steady state, while $\mathit{\delta }{q}_{\text{s}}^{*}$ and $\mathit{\delta }{q}_{\text{t}}^{*}$ denote the small perturbation. The exposed area can be similarly written as $${\displaystyle }{A}^{*}=\stackrel{\mathrm{‾}}{A^{*}}+\mathit{\delta }{A}^{*}.$$

Steady-state cover is directly related to the mass on the bed $M_{\text{s}}^{*}$ by Eq. (3), which, as long as $P$ is independent of time, we can rewrite as $${\displaystyle \frac{\text{d}{A}^{*}}{\text{d}t}}=-P{\displaystyle \frac{\text{d}{M}_{\text{s}}^{*}}{\text{d}t}}.$$

Substituting Eq. (A3) and a similar equation for $M_{\text{s}}^{*}$, $${\displaystyle }{M}_{\text{s}}^{*}=\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}+\mathit{\delta }{M}_{\text{s}}^{*},$$ we obtain $${\displaystyle \frac{\text{d}\mathit{\delta }{A}^{*}}{\text{d}t}}=-P{\displaystyle \frac{\text{d}\mathit{\delta }{M}_{\text{s}}^{*}}{\text{d}t}}.$$

Here, the averaged terms drop out as they are independent of time. If $P$ and the steady-state solution for $A^{*}$ are known, a direct relationship between $A^{*}$ and $M_{\text{s}}^{*}$ can be derived. For example, for the exponential cover model (Eq. 8), substituting Eqs. (A3) and (A5), we find $$\begin{array}{cc}& {\displaystyle }\stackrel{\mathrm{‾}}{A^{*}}+\mathit{\delta }{A}^{*}& {\displaystyle }=e^{-\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}-\mathit{\delta }{M}_{\text{s}}^{*}}=e^{-\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}}{e}^{-\mathit{\delta }{M}_{\text{s}}^{*}}=\stackrel{\mathrm{‾}}{A^{*}}{e}^{-\mathit{\delta }{M}_{\text{s}}^{*}}{\displaystyle }& {\displaystyle }\approx \stackrel{\mathrm{‾}}{A^{*}}\left(\mathrm{1},-,\mathit{\delta },M_{\text{s}}^{*}\right).\end{array}$$

Here, since the $\mathit{\delta }$ variables are small, we approximated the exponential term using a Taylor expansion to first order. We obtain perturbation of sediment supply $${\displaystyle }\mathit{\delta }{A}^{*}=-\stackrel{\mathrm{‾}}{A^{*}}\mathit{\delta }{M}_{\text{s}}^{*}.$$

It is therefore sufficient to derive the perturbation solution for $M_{\text{s}}^{*}$, the time evolution of which is given by Eq. (22). Eliminating $M_{\text{m}}^{*}$ using Eq. (24), we obtain $${\displaystyle \frac{\partial M_{\text{s}}^{*}}{\partial t^{*}}}=\left(\mathrm{1},-,e^{-q_{\text{s}}^{*}/U^{*}}\right)q_{\text{s}}^{*}-\left(\mathrm{1},-,e^{-M_{\text{s}}^{*}}\right)q_{\text{t}}^{*}.$$

Perturbation of sediment supply

First, let us look at a perturbation of sediment supply $q_{\text{s}}^{*}$, while other parameters are held constant. Substituting Eqs. (A1) and (A5) into Eq. (A9), we obtain $$\begin{array}{cc}& {\displaystyle \frac{\partial \mathit{\delta }{M}_{\text{s}}^{*}}{\partial t^{*}}}=& {\displaystyle }\left(\mathrm{1},-,e^{-\left(\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}},+,\mathit{\delta },q_{\text{s}}^{*}\right)/U^{*}}\right)\left(\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}},+,\mathit{\delta },q_{\text{s}}^{*}\right){\displaystyle }& {\displaystyle }-\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}-\mathit{\delta }{M}_{\text{s}}^{*}}\right)q_{\text{t}}^{*}.\end{array}$$

Again, since the $\mathit{\delta }$ variables are small, we can replace the relevant exponentials with a Taylor expansion to first order: $${\displaystyle }{e}^{-\mathit{\delta }{q}_{\text{s}}^{*}/U^{*}}\approx \mathrm{1}-{\displaystyle \frac{\mathit{\delta }{q}_{\text{s}}^{*}}{U^{*}}}.$$

A similar approximation applies for the exponential in $M_{\text{s}}^{*}$. Substituting Eq. (A11) into Eq. (A10), expanding the multiplicative terms, dropping terms of second order in the $\mathit{\delta }$ variables and rearranging, we get $$\begin{array}{cc}& {\displaystyle \frac{\partial \mathit{\delta }{M}_{\text{s}}^{*}}{\partial t^{*}}}=& {\displaystyle }\mathit{\delta }{q}_{\text{s}}^{*}\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}{U}^{*}},+,{\displaystyle \frac{\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}}{U^{*}}},e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}}\right){\displaystyle }& {\displaystyle }-\mathit{\delta }{M}_{\text{s}}^{*}\left(q_{\text{t}}^{*},-,\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}}\right),\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}\right).\end{array}$$

The perturbation is assumed to be sinusoidal $${\displaystyle }\mathit{\delta }{q}_{\text{s}}^{*}=d\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{p}}\right).$$

Here, $p$ is the period of the perturbation and $d$ is its amplitude. Note that, to be consistent with the assumptions previously made, $d$ needs to be small in comparison with the average sediment supply. Substituting, Eq. (A12) can be integrated to obtain the solution $$\begin{array}{cc}& {\displaystyle }\mathit{\delta }{M}_{\text{s}}^{*}=& {\displaystyle }{G}_{q_{\text{s}}^{*}}\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{P}},+,{\mathit{\phi }}_{q_{\text{s}}^{*}}\right)+C{\displaystyle }& {\displaystyle }\exp \left\{-,\left(q_{\text{t}}^{*},-,\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}}\right),\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}\right),{\displaystyle \frac{t}{T}}\right\},\end{array}$$ where $C$ is a constant of integration. The gain is given by $${\displaystyle }{G}_{q_{\text{s}}^{*}}={\displaystyle \frac{p}{T}}{\displaystyle \frac{\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}},+,\frac{\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}}{U^{*}},e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}}\right)d}{\sqrt{{\left(q_{\text{t}}^{*},-,\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}}\right),\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}\right)}^{\mathrm{2}}{\left(\frac{p}{T}\right)}^{\mathrm{2}}+\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}}}}$$

and the phase shift by $${\displaystyle }{\mathit{\phi }}_{q_{\text{s}}^{*}}={\tan }^{-\mathrm{1}}\left[-,{\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\frac{p}{T}\left(q_{\text{t}}^{*},-,\left(\mathrm{1},-,e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}}\right),\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}\right)}}\right].$$

Perturbation of transport capacity

The perturbation of the transport capacity $q_{\text{t}}^{*}$ is a little more complicated, since both $q_{\text{t}}^{*}$ and $U^{*}$ are explicitly dependent on hydraulics (e.g. shear stress; see Eqs. 43 and 44), and thus $U^{*}$ is implicitly dependent on $q_{\text{t}}^{*}$ and $\mathit{\delta }{q}_{\text{t}}^{*}$. To circumvent this problem, we expand the exponential term featuring $U^{*}\left(\mathit{\delta }{q}_{\text{t}}^{*}\right)$ in Eq. (A9) using a Taylor series expansion around $\mathit{\delta }{q}_{\text{t}}^{*}=\mathrm{0}$. $$\begin{array}{cc}& {\displaystyle }& {\displaystyle }\exp \left\{-,{\displaystyle \frac{q_{\text{s}}^{*}}{U^{*}\left(\mathit{\delta },q_{\text{t}}^{*}\right)}}\right\}{\displaystyle }\approx & {\displaystyle }\exp \left\{-,{\displaystyle \frac{q_{\text{s}}^{*}}{U^{*}\left(\mathit{\delta },q_{\text{t}}^{*},=,\mathrm{0}\right)}}\right\}\\ {\displaystyle }& {\displaystyle }\left[\mathrm{1},-,{\displaystyle \frac{q_{\text{s}}^{*}}{U^{{*}^{\mathrm{2}}}\left(\mathit{\delta },q_{\text{t}}^{*},=,\mathrm{0}\right)}},{\displaystyle \frac{\partial U^{*}}{\partial \mathit{\delta }{q}_{\text{t}}^{*}}},\left(\mathit{\delta },q_{\text{t}}^{*},=,\mathrm{0}\right),\mathit{\delta },q_{\text{t}}^{*}\right]\end{array}$$

Both $U^{*}$ and its derivative are constants when evaluated at $\mathit{\delta }{q}_{\text{t}}^{*}=\mathrm{0}$. We can thus write $$\begin{array}{cc}{\displaystyle }\exp \left\{-,{\displaystyle \frac{q_{\text{s}}^{*}}{U^{*}}}\right\}& {\displaystyle }=\exp \left\{-,{\displaystyle \frac{q_{\text{s}}^{*}}{\stackrel{\mathrm{‾}}{U^{*}}}}\right\}\left[\mathrm{1},-,{\displaystyle \frac{q_{\text{s}}^{*}}{{\stackrel{\mathrm{‾}}{U^{*}}}^{\mathrm{2}}}},\stackrel{\mathrm{‾}}{\left({\displaystyle \frac{\partial U^{*}}{\partial \mathit{\delta }{q}_{\text{t}}^{*}}}\right)},\mathit{\delta },q_{\text{t}}^{*}\right]\\ & {\displaystyle }& {\displaystyle }=\left[\mathrm{1},-,C_{\mathrm{0}},\mathit{\delta },q_{\text{t}}^{*}\right]e^{q_{\text{s}}^{*}/\stackrel{\mathrm{‾}}{U^{*}}}.\end{array}$$

Here, $C_{\mathrm{0}}$ is a constant. Proceeding as before by substituting Eqs. (A2), (A8), and (A17) into (A9), expanding exponential terms containing $\mathit{\delta }$ variables, dropping terms of second order in the $\mathit{\delta }$ variables, and rearranging, we obtain $$\begin{array}{cc}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\delta }{M}_{\text{s}}^{*}}{\partial t^{*}}}={\displaystyle }& {\displaystyle }\left(B,q_{\text{s}}^{*},e^{-q_{\text{s}}^{*}/\stackrel{\mathrm{‾}}{U^{*}}},+,e^{-\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}},-,\mathrm{1}\right)\mathit{\delta }{q}_{\text{t}}^{*}-\mathit{\delta }{M}_{\text{s}}^{*}\stackrel{\mathrm{‾}}{q_{\text{t}}^{*}}{e}^{-\stackrel{\mathrm{‾}}{M_{\text{s}}^{*}}}.\end{array}$$ A sinusoidal perturbation of the form $${\displaystyle }\mathit{\delta }{q}_{\text{t}}^{*}=d\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{p}}\right)$$ yields the solution $$\begin{array}{cc}& {\displaystyle }\mathit{\delta }{M}_{\text{s}}^{*}& {\displaystyle }=G_{q_{\text{t}}^{*}}\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{P}},+,{\mathit{\phi }}_{q_{\text{t}}^{*}}\right){\displaystyle }& {\displaystyle }+C\exp \left\{-,\left(\stackrel{\mathrm{‾}}{q_{\text{t}}^{*}},-,\left(\mathrm{1},-,e^{-q_{\text{s}}^{*}/\stackrel{\mathrm{‾}}{U^{*}}}\right),q_{\text{s}}^{*}\right),{\displaystyle \frac{t}{p}}\right\}\\ {\displaystyle }& {\displaystyle }\left\{-,\left(\stackrel{\mathrm{‾}}{q_{\text{t}}^{*}},-,\left(\mathrm{1},-,e^{-q_{\text{s}}^{*}/\stackrel{\mathrm{‾}}{U^{*}}}\right),q_{\text{s}}^{*}\right),{\displaystyle \frac{t}{T}}\right\},\end{array}$$ with $${\displaystyle }{G}_{q_{\text{t}}^{*}}={\displaystyle \frac{p}{T}}{\displaystyle \frac{\left(\frac{q_{\text{s}}^{{*}^{\mathrm{2}}}}{{\stackrel{\mathrm{‾}}{U^{*}}}^{\mathrm{2}}},\stackrel{\mathrm{‾}}{\left(\frac{\partial U^{*}}{\partial \mathit{\delta }{q}_{\text{t}}^{*}}\right)},e^{-q_{\text{s}}^{*}/\stackrel{\mathrm{‾}}{U^{*}}},-,\left(\mathrm{1},-,e^{-q_{\text{s}}^{*}/\stackrel{\mathrm{‾}}{U^{*}}}\right),\frac{q_{\text{s}}^{*}}{\stackrel{\mathrm{‾}}{q_{\text{t}}^{*}}}\right)d}{\sqrt{{\stackrel{\mathrm{‾}}{q_{\text{t}}^{*}}}^{\mathrm{2}}{\left(\frac{p}{T}\right)}^{\mathrm{2}}{\left(\mathrm{1},-,\left(\mathrm{1},-,e^{-q_{\text{s}}^{*}/U^{*}}\right),\frac{q_{\text{s}}^{*}}{\stackrel{\mathrm{‾}}{q_{\text{t}}^{*}}}\right)}^{\mathrm{2}}+\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}}}}$$ and $${\displaystyle }\mathit{\phi }={\tan }^{-\mathrm{1}}\left(-,{\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\frac{p}{T}\left(\stackrel{\mathrm{‾}}{q_{\text{t}}^{*}},-,\left(\mathrm{1},-,e^{-q_{\text{s}}^{*}/\stackrel{\mathrm{‾}}{U^{*}}}\right),q_{\text{s}}^{*}\right)}}\right).$$

Summary

Using the system timescale $T_{\text{S}}$, the phase shift and gain can be generally rewritten as $${\displaystyle }\mathit{\phi }={\tan }^{-\mathrm{1}}\left(-,\mathrm{2},\mathit{\pi },{\displaystyle \frac{T_{\text{S}}}{p}}\right),$$ $${\displaystyle }G={\displaystyle \frac{p}{T_{\text{S}}}}{\displaystyle \frac{Kd}{\sqrt{{\left(\frac{p}{T_{\text{S}}}\right)}^{\mathrm{2}}+\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}}}}.$$

Here, $K$ differs for perturbations in sediment supply and transport capacity, given by the equations $${\displaystyle }{K}_{q_{\text{s}}^{*}}=\mathrm{1}-e^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}}+{\displaystyle \frac{\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}}{U^{*}}}{e}^{-\stackrel{\mathrm{‾}}{q_{\text{s}}^{*}}/U^{*}},$$ $${\displaystyle }{K}_{q_{\text{t}}^{*}}={\displaystyle \frac{q_{\text{s}}^{{*}^{\mathrm{2}}}}{{\stackrel{\mathrm{‾}}{U^{*}}}^{\mathrm{2}}}}\stackrel{\mathrm{‾}}{\left({\displaystyle \frac{\partial U^{*}}{\partial \mathit{\delta }{q}_{\text{t}}^{*}}}\right)}{e}^{-q_{\text{s}}^{*}/\stackrel{\mathrm{‾}}{U^{*}}}-\left(\mathrm{1},-,e^{-q_{\text{s}}^{*}/\stackrel{\mathrm{‾}}{U^{*}}}\right){\displaystyle \frac{q_{\text{s}}^{*}}{\stackrel{\mathrm{‾}}{q_{\text{t}}^{*}}}}.$$

Notation

Overbars denote time-averaged quantities.
$a$	Shape parameter in the regularized
	incomplete Beta function
$A^{*}$	Fraction of exposed (uncovered) bed area
$A_c^{*}$	Fraction of covered bed area
$b$	Shape parameter in the regularized
	incomplete Beta function
$B$	Regularized incomplete Beta function
$C$	Constant of integration
$C_{\mathrm{0}}$	Constant (m${}^{\mathrm{2}}$s kg${}^{-\mathrm{1}}$)
$d$	Amplitude of perturbation (kg m${}^{-\mathrm{2}}$s)
$D$	Sediment deposition rate per bed area (kg m${}^{-\mathrm{2}}$s)
$D_{\text{tot}}$	Sediment deposition rate (kg s${}^{-\mathrm{1}}$)
$D^{*}$	Dimensionless sediment deposition rate
$D_{\mathrm{50}}$	Median grain size (m)
$e$	Base of the natural logarithm
$E$	Sediment entrainment rate per bed area (kg m${}^{-\mathrm{2}}$s)
$E_{\text{tot}}$	Sediment entrainment rate (kg s${}^{-\mathrm{1}}$)
$E^{*}$	Dimensionless sediment entrainment rate
$E_{\text{max}}$	Maximal possible dimensionless
	sediment entrainment rate
$g$	Acceleration due to gravity (m s${}^{-\mathrm{2}}$)
$G$	Gain (kg m${}^{-\mathrm{2}}$s)
$I$	Non-dimensional incision rate
$k$	Probability of sediment deposition on uncovered
	parts of the bed, linear implementation
$K$	Parameter in the gain equation
$L$	Characteristic length scale (m)
$M_{\mathrm{0}}$	Minimum mass per area necessary to cover
	the bed (kg m${}^{-\mathrm{2}}$)
$M_{\mathrm{0}}^{*}$	Dimensionless characteristic sediment mass
$M_{\text{m}}$	Mobile sediment mass (kg m${}^{-\mathrm{2}}$)
$M_{\text{m}}^{*}$	Dimensionless mobile sediment mass
$M_{\text{s}}$	Stationary sediment mass (kg m${}^{-\mathrm{2}}$)
$M_{\text{s}}^{*}$	Dimensionless stationary sediment mass
$p$	Period of perturbation (s)
$p_{\text{c}}$	Probability of entrainment, CA model, blocked grains
$p_{\text{i}}$	Probability of entrainment, CA model, free grains
$P$	Probability of sediment deposition
	on uncovered parts of the bed
$q_{\text{s}}$	Mass sediment transport rate per unit width (kg ms${}^{-\mathrm{1}}$)
$q_{\text{s}}^{*}$	Dimensionless sediment transport rate
$q_{\text{t}}$	Mass sediment transport capacity
	per unit width (kg ms${}^{-\mathrm{1}}$)
$q_{\text{t}}^{*}$	Dimensionless transport capacity
$Q_{\text{s}}^{*}$	Relative sediment supply; sediment transport
	rate over transport capacity
$Q_{\text{t}}$	Mass sediment transport capacity (kg s${}^{-\mathrm{1}}$)
$t$	Time variable (s)
$t^{*}$	Dimensionless time
$T$	Characteristic timescale (s)

$T_{\text{E}}$	Characteristic timescale for
	sediment entrainment (s)
$T_{\text{S}}$	Characteristic system timescale (s)
$U$	Sediment speed (m s${}^{-\mathrm{1}}$)
$U^{*}$	Dimensionless sediment speed
$x$	Dimensional streamwise spatial coordinate (m)
$x^{*}$	Dimensionless streamwise spatial coordinate
$y$	Dummy variable
$\mathit{\alpha }$	Exponent
$\mathit{\gamma }$	Fraction of pore space in the sediment
$\mathit{\delta }$	Denotes time-varying component
$\mathrm{\Delta }$in	Sediment supply rate from
	upstream direction (kg s${}^{-\mathrm{1}}$)
$\mathrm{\Delta }$out	Transport rate of sediment leaving in
	the downstream direction (kg s${}^{-\mathrm{1}}$)
$\mathrm{\Delta }{M}_{\text{m}}$	Change in mobile sediment mass (kg)
$\mathrm{\Delta }t$	Change in time (s)
$\mathit{\theta }$	Shields stress
${\mathit{\theta }}_{\text{c}}$	Critical Shields stress
$\mathit{\rho }$	Density of water (kg m${}^{-\mathrm{3}}$)
${\mathit{\rho }}_{\text{s}}$	Density of sediment (kg m${}^{-\mathrm{3}}$)
$\mathit{\tau }$	Bed shear stress (N m${}^{-\mathrm{2}}$)
${\mathit{\tau }}_{\text{c}}$	Critical bed shear stress at the onset of
	bedload motion (N m${}^{-\mathrm{2}}$)

The authors declare that they have no conflict of interest.

Acknowledgements

We thank J. Scheingross and J. Braun for insightful discussions and two anonymous reviewers and associate editor D. Egholm for their comments on the paper. The data from the Erlenbach is owned by and was used with permission of the Mountain Hydrology and Mass Movements Group at the Swiss Federal Research Institute for Forest Snow and Landscape Research WSL.The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association. Edited by: David Lundbek Egholm Reviewed by: two anonymous referees