Fluvial response to changes in the magnitude and frequency of sediment supply in a 1-D model

In steep headwater reaches, episodic mass movements can deliver large volumes of sediment to fluvial channels. If these inputs of sediment occur with a high frequency and magnitude, the capacity of the stream to rework the supplied material can be exceeded for a significant amount of time. To study the equilibrium conditions in a channel following different episodic sediment supply regimes (defined by grain size distribution, frequency, and magnitude of events), we simulate sediment transport through an idealized reach with our numerical 1-D model “BESMo” (Bedload Scenario Model). The model performs well in replicating flume experiments of a similar scope (where sediment was fed constantly, in one, two, or four pulses) and allowed the exploration of alternative event sequences. We show that in these experiments, the order of events is not important in the long term, as the channel quickly recovers even from high magnitude events. In longer equilibrium simulations, we imposed different supply regimes on a channel, which after some time leads to an adjustment of slope, grain size, and sediment transport that is in equilibrium with the respective forcing conditions. We observe two modes of channel adjustment to episodic sediment supply. (1) High-frequency supply regimes lead to equilibrium slopes and armouring ratios that are like conditions in constant-feed simulations. In these cases, the period between pulses is shorter than a “fluvial evacuation time”, which we approximate as the time it takes to export a pulse of sediment under average transport conditions. (2) In low-frequency regimes the pulse period (i.e., recurrence interval) exceeds the “fluvial evacuation time”, leading to higher armouring ratios due to the longer exposure of the bed surface to flow. If the grain size distribution of the bed is fine and armouring weak, the model predicts a decrease in the average channel slope. The ratio between the “fluvial evacuation time” and the pulse period constitutes a threshold that can help to quantify how a system responds to episodic disturbances.

1 Introduction

Mass movements in mountainous regions often deliver sediment directly to the stream network, resulting in coupled conditions that can trigger immediate channel responses during relatively large delivery events. Notably, delivery events can reset the local channel profile and govern construction and maintenance of channel bed architecture downstream of delivery points. The local response rate and trajectory following a delivery event is a function of the prevailing watershed flow regime, the magnitude of the delivery event, gradients in channel width (Ferrer-Boix et al., 2016), and in some instances the concentrated activity of aquatic species such as salmon (Hassan et al., 2008a). While lowland river systems have been the focus of a substantial body of research, less work has been carried out within steep mountain streams, particularly concerning mountain channel responses to changes in flow or sediment supply regimes.

In mountain streams, large, episodic inputs may temporarily dominate channel processes and morphology, significantly altering sediment transport and storage within the stream channel (Hassan et al., 2005, 2008b). Lisle et al. (1997) conducted flume experiments and numerical modelling that showed that sediment pulses are mainly reworked in situ, in contrast to a downstream translation in the form of a sediment wave. This finding is supported by Lisle et al. (2001), where little evidence of sediment waves was found in the field. Lisle and Church (2002) suggested that a stream channel responds to changes in the sediment supply by altering both storage and sediment transport rates. They describe that after a sediment pulse occurs, a first phase with low armouring rates allows for high transport rates in reworking the introduced material, corresponding to supply limited conditions. This is followed by a second phase in which armouring develops and transport rates decrease, corresponding to transport limited conditions. If the fluvial capacity is too low to evacuate the mass or grain size of material in the current hydrological regime, lag sediment can remain in the channel and dominate the local morphology for a long time (Benda et al., 2005; Brummer and Montgomery, 2006). Patterns of cyclic behaviour, associated with the rapid input of sediment from external sources, have been described in a number of field observations (Hassan et al., 2008b; Hoffman and Gabet, 2007; Madej, 1999, 2001; Madej and Ozaki, 1996; Miller and Benda, 2000; Roberts and Church, 1986) and some experimental studies (An et al., 2017a; Cui et al., 2003; Elgueta, 2014; Elgueta-Astaburuaga and Hassan, 2017; Ferrer-Boix and Hassan, 2014, 2015; Johnson et al., 2015; Sklar et al., 2009; Venditti et al., 2010; von Flotow, 2013). These studies observed fining of the bed surface, higher mobility and thus increased transport rates following such episodic sediment supply events. A reverse trend for coarsening and stabilizing of the bed was also noted as the supply was exhausted and a decrease in sediment transport rates followed (Church et al., 1998; Dietrich et al., 1989; Hassan and Church, 2000; Nelson et al., 2009). Cui and Parker (2005) support these findings using numerical modelling and further point out that abrasion can play an important role in the reworking of sediment pulses. Field observations (Benda, 1990; Pryor et al., 2011), and flume experiments (Luzi, 2014; Pryor et al., 2011) also document cycles of aggradation and degradation due to changes in the sediment supply. The observations discussed above suggest that changes in the sediment supply rate may lead to significant changes in bed elevation, bed surface texture, channel stability, and bed morphology. An analytical model developed by Blom et al. (2017) showed that the local channel geometry and surface grain size composition is mainly governed by long-term mean sediment supply rates and not by short-term changes in supply conditions. In contrast to cyclic sediment supply, cyclic hydrographs were found to mainly affect sediment transport rates and have a lesser impact on bed surface texture and channel morphology (Parker et al., 2007). Wong and Parker (2006) reported that cyclic hydrographs cause a part of the channel bed to undergo cyclic aggradation and degradation forming a hydrograph boundary layer. Cyclic sediment supply causes a similar effect which is termed the sedimentograph boundary layer (An et al., 2017b).

These findings imply that the morphological impact of the sediment pulse is most prevalent at the point of entrance, while the downstream portion mostly conveys the subsequently eroded material. The time needed for channel adjustments to occur after a large sediment input event depends on the amount and the texture of the delivered material. Brummer and Montgomery (2006) reported that 2 years after the supply event, the most mobile fractions (e.g., the fine fractions) were evacuated by a series of moderate floods while the largest grain sizes remained in the channel as lag deposits because flows were below their flow competence. Further, they showed that selective transport of sediment led to the development of an armour layer after only a few flow events. This armour layer protects the supplied material in the bed subsurface, increases bed stability, and causes lateral erosion and channel widening.

We expect the sediment transport rate in a channel to reach a long-term balance between erosional and depositional forces, even though there can be periodic changes in the short term. This state is defined as a “dynamic equilibrium” following Ahnert (1994). If the external forcing on the system changes, the channel will be in a transient state of adjustment towards a new dynamic equilibrium. Little attention has been directed to the question of what effect a change in the frequency of sediment supply events may have on the response of alluvial streams. Brunsden and Thornes (1979) proposed that if the frequency (i.e., recurrence interval) of disturbing events is shorter than the time necessary for a system to adjust to new boundary conditions (“relaxation time”), then transience will dominate the system and it may never achieve equilibrium (Brunsden, 1980). Wolman and Miller (1960) suggest that mountain channels experiencing direct inputs of sediment are good examples of such systems, where form is defined by extreme events rather than events of intermediate magnitude and frequency. The concept of this so-called “temporal sensitivity” was later elaborated on by Thomas (2001) and Brunsden (2001), although since their studies little attempt has been made in fluvial geomorphology to address this issue in practice. Bull (1991) applies the theory of Brunsden and Thornes (1979) to the impact of a hypothetical temporal succession of disturbance events on sediment storage, which can either increase or decrease the stream-bed elevation. The system processes the disturbance over the relaxation time, which together with a potential reaction time constitutes the total response time. If there is no further disturbance, the system status remains unchanged over a time of persistence. The concept from Bull (1991) is based on a system's trend towards a dynamic equilibrium between the forcing by, and the reworking of, disturbances. Howard (1982) concluded that if episodic inputs occur with a frequency that matches the inverse of the relaxation time, the output of the system will remain in a constant equilibrium with the average value of the forcing. Flume based insights about equilibrium conditions and timescales of adjustment to changes in sediment supply rates are discussed in some studies (Elgueta-Astaburuaga and Hassan, 2017; Pryor et al., 2011), but response times are only quantified in few cases (Podolak and Wilcock, 2013).

The paragraphs above illustrate how the frequency at which events occur may be fundamental in defining the response of a fluvial system to a change in boundary conditions. Therefore, it appears that event frequency should be a central aspect in investigations regarding the effect of episodic sediment supply on streams. Consequently, our understanding of involved processes remains incomplete. For example, it is uncertain whether the freshly delivered sediment that buries and is transferred over the bed surface is simply removed by the subsequent floods or whether there is some exchange between armoured and structured bed and the fine and mobile deposits. Furthermore, Hassan and Zimmermann (2012) asserted that it is important to study how quickly internal changes in grain size, channel morphology, and sediment storage occur when the stream shifts between cycles of aggradation and degradation.

Our main research objective is to describe the impact of episodic sediment supply on channel bed evolution in simulations using a 1-D morphodynamic numerical model for a bed of multiple grain sizes. We use the model to recreate conditions from experiments conducted at the Mountain Channel Hydraulic Experimental Laboratory, University of British Columbia, where a set of experiments were carried out to examine the impacts of episodic sediment supply on bed surface evolution and channel adjustment of a gravel bed stream (Elgueta, 2014; Elgueta-Astaburuaga, 2018; Elgueta-Astaburuaga and Hassan, 2017; von Flotow, 2013). Although these experiments provide detailed information on channel adjustment to changes in the sediment supply regime, they are limited in terms of the number of experiments and the range of scenarios that could be conducted. The performance of the model is tested against the experimental results obtained in the laboratory and then used to further explore controls and responses of the fluvial system to changes in flow and sediment supply regimes. The specific research questions addressed in this study are as follows:

1. Can the numerical model recreate the channel response that was observed in flume experiments of similar scope?

2. Does the sequencing of supply events play a role in the reaction of a gravel-bed stream, when several events of specified magnitudes occur in a different order?

3. How will different combinations of episodic sediment supply, obtained by varying their magnitude and frequency, impact channel evolution of a gravel-bed stream?

2 Methods

We applied the 1-D morphodynamic model BESMo (Bedload Scenario Model) to calculate capacity based sediment transport under different sediment supply regimes. We chose values for model parameters to match the flume experiments as closely as possible and used measurements of sediment transport rate, surface grain size distribution, and slope to calibrate the model. Matching our research questions, we then conducted two types of simulations:

1. In “event sequencing simulations”, we simulated different permutations of events to understand the role event succession plays in long-term channel response.

2. In “equilibrium simulations”, we used the same model setup and imposed different, but within each run regular, supply event frequencies. These simulations were run until we achieved a recurrent pattern in slope and grain size adjustment, allowing us to identify how the channel adjusts to the supply regime in the long term.

2.1 Model setup

The structure of the model is similar to other models designed to reproduce and interpret data from flume experiments (An et al., 2017a; Cui and Parker, 2005; Ferrer-Boix and Hassan, 2014; Wong and Parker, 2006). Figure 1 gives an overview of the implemented model components and their basic interaction. The model can be subdivided into a “hydraulic part” and a “sediment part”, both of which are subject to “external forcing” that varies in according to the modelling scenarios.

Figure 1 Flowchart stating the main components of the model and the flow of information between them. The temporal loop is advanced as the new elevation affects the slope in the flow model. The components are coloured as follows: blue – flow related, dark yellow – sediment volume related, peach – particle size related, and green – geometry related.

We set up the modelling environment to run on a Compute Canada research cluster, which allows us to simulate many different input conditions in parallel and compare results quickly. We use a backwater flow model as suggested by Cui et al. (2006), implementing a threshold Froude number (Fr) to switch conditions between supercritical and subcritical flow:

$$\begin{array}{}\text{(1)}& \mathit{Fr}=\sqrt{{\displaystyle \frac{Q_{\mathrm{w}}^{\mathrm{2}}}{g{w}_{\mathrm{r}}^{\mathrm{2}}{h}^{\mathrm{3}}}}}.\end{array}$$

The Froude number (Fr) is calculated as a function of discharge (Q_w), gravity (g), channel width (w_r) and water depth (h) (Eq. 1). The threshold Fr=0.9 simplifies the calculation of flow conditions, allowing us to spatially iterate through the nodes only once from downstream to upstream. In the case of subcritical flow, the water depth is solved locally as a function of downstream friction slope (S_f) and bed slope (S₀) (Eq. 2a):

$$\begin{array}{}\text{(2a)}& {\displaystyle \frac{\mathrm{d}h}{\mathrm{d}x}}={\displaystyle \frac{S_{\mathrm{0}}-S_{\mathrm{f}}}{\mathrm{1}-{\mathit{Fr}}^{\mathrm{2}}}}& {\displaystyle }\text{for subcritical: }\mathit{Fr}<\mathrm{0.9}\text{(2b)}& {\displaystyle }h=({\displaystyle \frac{n^{\mathrm{2}}{Q}_{\mathrm{w}}^{\mathrm{2}}}{S_{\mathrm{f}}}}{)}^{\mathrm{3}/\mathrm{10}}& {\displaystyle }\text{for supercritical: }\mathit{Fr}\ge \mathrm{0.9}.\end{array}$$

Water depth under supercritical flow conditions is calculated locally assuming steady uniform flow using the Manning–Strickler formulation (Eqs. 2b and 3), where α_r is a coefficient of 8.1 (Parker, 1991) and roughness height (k_s):

$$\begin{array}{}\text{(3)}& n={\displaystyle \frac{k_s^{\mathrm{1}/\mathrm{6}}}{{\mathit{\alpha }}_{\mathrm{r}}{g}^{\mathrm{1}/\mathrm{2}}}}.\end{array}$$

k_s is calculated using the constants n_k and D_s90, the surface grain size for which 90 % of the surface is finer:

$$\begin{array}{}\text{(4)}& k_s=n_k{D}_{\mathrm{s}\mathrm{90}}.\end{array}$$

In the case of steady and uniform flow the bed slope S₀ is equal to the friction slope (S_f) with bed elevation (η_b) at downstream position x:

$$\begin{array}{}\text{(5)}& S_{\mathrm{f}}={\displaystyle \frac{-d{\mathit{\eta }}_{\mathrm{b}}}{\mathrm{d}x}}=S_{\mathrm{0}}.\end{array}$$

If the solution of the water depth with Eq. (2a) is numerically unstable on the current node distribution, the model subdivides the channel into more nodes and reiterates the subdivision until a stable backwater curve is found. This approach does not properly represent the location of hydraulic jumps (Cui et al., 2006), which should not be a problem as we average conditions over a node spacing of at least one channel width. Boundary shear stress (τ_b) is then calculated with the depth-slope product:

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{\mathrm{b}}=\mathit{\rho }gh{S}_{\mathrm{f}},\end{array}$$

where ρ is the water density. τ_b is then converted to the shear velocity u^*, which is used in the sediment routing component:

$$\begin{array}{}\text{(7)}& u^{*}=\sqrt{{\mathit{\tau }}_{\mathrm{b}}/\mathit{\rho }}.\end{array}$$

The volumetric unit bedload transport rate per size class q_bi is calculated using the sediment transport function provided by Wilcock and Crowe (2003). The change in bed elevation ∂η_b for each node x per time step t follows from the Exner equation of mass conservation:

$$\begin{array}{}\text{(8)}& (\mathrm{1}-\mathit{\lambda }){\displaystyle \frac{\partial {\mathit{\eta }}_{\mathrm{b}}}{\partial t}}=-{\displaystyle \frac{\partial q_{\mathrm{b}}}{\partial x}},\end{array}$$

where λ is bed porosity. The volumetric bedload transport rate per unit width is given as q_b and is calculated per grain size class i in the sediment mixture of n size classes:

$$\begin{array}{}\text{(9)}& q_{\mathrm{b}}=\sum _{i=\mathrm{1}}^n{q}_{\mathrm{b}i}.\end{array}$$

The model incorporates subsurface stratigraphy using the active layer concept (Parker, 2008), which gives the Hirano equation:

$$\begin{array}{}\text{(10)}& (\mathrm{1}-\mathit{\lambda })\left[{\displaystyle \frac{\partial }{\partial t}}\left(L_{\mathrm{a}}{F}_i\right)-f_{Ii}{\displaystyle \frac{\partial L_{\mathrm{a}}}{\partial t}}\right]=-\left({\displaystyle \frac{\partial q_{\mathrm{b}i}}{\partial x}}-f_{Ii}{\displaystyle \frac{\partial q_{\mathrm{b}}}{\partial x}}\right),\end{array}$$

where L_a is the active layer thickness, F_i is the surface frequency of the ith grain size class, f_Ii is the ith grain size class proportion exchanged between the surface and the subsurface, and q_bi=p_biq_b is the volumetric unit bedload transport rate of the ith grain size class where p_bi represents the ith fraction of the bedload transport rate. The active layer thickness is calculated as L_a=n_aD_s90, with the parameter n_a, representing the scale of bed fluctuations.

The grain size distribution of the sediment flux between the active layer and the substrate is either calculated from the subsurface texture when the bed degrades, or from a linear combination of surface and bedload grain size distributions when the bed aggrades (Hoey and Ferguson, 1994):

$$\begin{array}{}\text{(11)}& f_{Ii}=\left\{\begin{array}{ll}{f}_i& \text{for }\frac{\partial {\mathit{\eta }}_{\mathrm{b}}}{\partial t}<\mathrm{0}\\ \mathit{\alpha }{F}_i+(\mathrm{1}-\mathit{\alpha })p_{\mathrm{b}i}& \text{for }\frac{\partial {\mathit{\eta }}_{\mathrm{b}}}{\partial t}>\mathrm{0},\end{array}\right.\end{array}$$

where f_i is the fraction of the ith grain size class in the subsurface, and α is a constant. The vertical stratigraphy is stored in 10 cm high layers following Viparelli et al. (2010). By keeping track of grain size distributions within the surface and subsurface layers, the model can preserve the history of phases of erosion or aggradation. This allows emergent properties such as armouring layers to occur. To study the combined effect that the active layer thickness factor n_a and the active layer exchange ratio α have on the model results, we executed sensitivity runs shown in the Supplement to this paper (Supplement Figs. S1–S4).

The Exner equation (Eq. 8) in combination with the expression for the friction slope in Eq. (2) and the sediment transport function by Wilcock and Crowe (2003) form an non-linear advection–diffusion system that allows the calculation of bed elevation as a function of space and time (An et al., 2017a). An upwind scheme was used for the numerical discretization. The model needs an initial bed profile and an initial value of the surface grain size to be solvable. Sediment boundary conditions are given by the sediment feed rate and grain size distributions on the inlet, and a fixed bed elevation at the outlet of the simulated reach. The flow boundary condition is a water surface height of 0.1 m over normal flow at the outlet. The bedload transport function is used to calculate transport rates for each channel cross section. During the model run, changes of the sediment transport rate are dependent on changes in the sediment supply, channel slope, and surface grain size distribution.

2.2 Model calibration

We used data from flume experiments to calibrate the model. The objective of the flume experiments was to measure the adjustment of an alluvial steep channel to different frequencies of sediment supply. The experiments were carried out in a water recirculating flume which is 18 m long, 1 m wide, and 1 m deep in the Mountain Channel Hydraulic Experimental Laboratory at the University of British Columbia. Here we will provide a brief summary of the flume setup and experimental design, for more details see von Flotow (2013), Elgueta (2014), Ferrer-Boix and Hassan (2015), and Elgueta-Astaburuaga and Hassan (2017).

The experiment consisted of seven 40 h long runs with different sediment supply frequencies, while keeping the total sediment input the same at 300 kg per run. The experiments were run continuously, i.e., the bed surface at the end of run 1 was the starting condition for run 2 and so on. For all runs, flow was held constant so that the sediment feed regime could be studied with no changes in flow regime (Table 1). The difference between the runs was the spreading of the supply over a changing input frequency, which was either constant, in one pulse, in two pulses, or in four pulses. The bed was fixed in the first 1 m downstream section of the flume head box with stones equivalent to about D_s84 of the experimental bed material. In the remainder of the flume the bed initially consisted of 0.1 m of loose material with particle sizes ranging from 0.5 to 64 mm with a D_s50 of 5.64 mm, matching downscaled (by a factor of 3) conditions of a study reach in East Creek, British Columbia, Canada. The flume slope was set to 0.022 m m⁻¹. Measurements include water depth, water surface slope, water velocity, bed surface slope, bed surface particle size distribution, bed elevation, sediment transport rate, and bedload texture. Measurements of the water surface elevation were conducted throughout the experiment using a mechanical point gauge with 0.001 m precision. Photos were used to manually sample bed surface grain size distributions and the bed elevation was recorded with a green laser scanner at a 2 mm resolution. The bed surface scans were used to measure the bed surface slope along the thalweg, i.e., the line of lowest elevation along the flume. Flow velocity measurements were conducted using an ADV profiler. The grain size and count of particles exiting the flume were recorded with a camera and a light table at the outlet of flume (Zimmermann et al., 2008). The transport rate measurements were done at 30 Hz and validated after the experiments by total exported weight.

Using the model described, we simulated a 12 m long and 1 m wide channel in 13 downstream nodes each spaced 1 m apart. The model was set to calculate sediment transport for all nodes in time steps of 10 s. All simulations used a constant water discharge of 0.065 m³ s⁻¹ and a geometric mean grain size of 5.64 mm for both initial bed and sediment feed (full distribution shown in Fig. 2a). We chose to use a normally distributed approximation of the flume grain size with a width of σ=1.6 to be consistent with distributions used for the equilibrium simulations. This distribution and the original flume GSD are statistically the same. The initial channel slope was 0.022 m m⁻¹, also matching the parameters from the flume experiments.

We calibrated the model by visually reducing the difference between measured and simulated values of bed slope (S), surface grain size parameters (D_sg, D_s90), and transport rate (q_b). In the calibration runs more importance was given to recreating S, D_sg, and q_b than to a good match in D_s90. We first increased the reference Shields stress (${\mathit{\tau }}_{\mathrm{rm}}^{*}$) in the Wilcock and Crowe formula to roughly match simulated and measured q_b. Afterwards, we varied the grain exchange ratio (α) (Eq. 11) and the coefficient of the active layer thickness (n_a). As we achieved a good visual match between simulation and flume measurements, we did not see the need to calibrate more parameters.

Table 1 Overview of runs in the flume experiments. All runs were 40 h long. The texture of the initial bed mixture and the sediment feed were identical. Plots of resulting slope, D_sg, D_s90, and sediment transport are shown in Fig. 3. The symbols in column two signify the feed regime where 0 represents no feed, C represents constant feed, and 1, 2, and 4 represent the number of pulses in 40 h.

Initial slope 2.2 %; mixture D_s50 =5.64 mm; mixture D_s90 =11.2 mm, Q_w=65 L s⁻¹; duration of each run =40 h; total feed =300 kg per run. NA represents not available.

2.3 Event sequencing simulations

We explored the role that the sequencing of the pulse events could have on the flume study by simulating the “original flume” event sequence and comparing the result to alternative sequencing of events (see Table 2). The alternative event sequences are using the same pulse distributions (four pulses, two pulses, or one pulse over 40 h), but the pulse order is either from “few to many” (FtM) (i.e., one pulse, then two pulses, then four pulses) or from “many to few” (MtF) (i.e., four pulses, then two pulses, then one pulse) per 40 h phase. To allow the system more time to recover from pulse events, we simulated two more cases where each pulsed phase is buffered from the next one by a 40 h constant-feed phase. These runs are called “c-buffered: many to few” (cMtF) and “c-buffered: few to many” (cFtM).

Table 2 Sequencing of events in runs that were simulated as permutations of the original flume experiment. Each of the seven periods was 40 h long and each run lasted 280 h in total. There was no sediment input during the no feed runs (0). Within all other period types, 300 kg of sediment was fed over 40 h, either constantly (C) or in pulses (one, two, or four events). Besides recreating the original flume sequence (OF), we simulated two runs where the pulsed events occur either in order from many pulses to few (MtF) or from few to many (FtM). To explore if the system could rebound from the impact of a certain pulse phase during a constant-feed phase, we created two additional runs where this was the case (cMtF and cFtM), which led to a 600 kg higher total sediment feed.

2.4 Equilibrium simulations

In our final set of experiments, we kept the frequency and magnitude of pulse events constant to achieve equilibrium slope and grain size conditions. The use of numerical modelling allows for the comparison of many simulations with differing grain size distributions, pulse frequencies, and pulse magnitudes. We expect a channel under episodic sediment supply to adjust synchronously to the frequency of external forcing events. The added sediment volume from a supply event will increase the channel slope at first. After the supply of sediment ends and material is removed from the channel, the slope will decrease, and the surface grain size will begin to reflect the sediment starved conditions. As the long-term sediment input equals the long-term sediment output, the channel will eventually achieve a condition where the capacity to erode material (through increased slope in conjunction with changes in the surface grain size) equals the depositional forcing (i.e., long-term sediment input) of the supply regime. In this state the adjustment of channel slope and grain size to each sediment input event will return to the same values after every pulse. All runs achieved this equilibrium condition within 20 000 simulation hours.

Figure 2 (a) Variation of the grain size distribution between model runs (lines) and values from the flume experiments (circles). All distributions have a geometric mean grain size (D_sg) of 5.64 mm. The different σ values represent the width of the distribution, calculated for normal distributions in phi-scaled sediment sizes. The data for the flume experiments are roughly matched with σ=1.6. (b) Combinations of pulse magnitude and pulse period (or recurrence interval) used in the equilibrium model runs. As the total mass of supplied material is the same for all simulations, runs with high pulse frequencies (i.e., low pulse periods T_pp) are smaller in magnitude. Four of the combinations match the flume runs (circles).

To find the equilibrium slope resulting from different sediment supply regimes, we simulated different combinations of sediment supply frequency (F_pulse) and magnitude (M_pulse) for nine different grain size distributions (Fig. 2a). All distributions have the same mean grain size of 5.64 mm, but differ in the width of the distribution by the standard deviation (σ), which was chosen for a phi-scaled, normally distributed sample. While σ=0.05 represents a nearly uniform sediment mixture, σ=1.6 roughly matches the grain size distribution of the flume experiments. We used 11 grain size classes in the simulations and set the initial GSD to the feed GSD. Under the hydraulic conditions applied all grain sizes are initially mobile, which might change during the simulations due to the effect of armouring and changes in bed slope. Figure 2b shows the combinations of frequency and magnitude used in this study. Each model run delivered the same input mass over the simulation time (150 000 kg over 20 000 h). We then distributed this total mass over different pulse frequencies, with four of the combinations matching the flume experiments. Each pulse was 10 min in length. The lowest frequency was chosen to be one pulse every 400 h (pulse period time: T_pp=400 h), and the highest frequency was constant feed (one 10 min long pulse every 10 min). We selected a range of 40 pulse frequencies for which the whole number of cycles summed to 20 000 h. In total we executed 360 simulations, 9 different σ with 40 frequencies each.

3 Results

3.1 Model calibration

Elgueta-Astaburuaga and Hassan (2017) describe the flume results we used for model calibration in detail. We will give a short summary of the findings here. An initial run without sediment feed over an unstructured bed showed a high sediment output while the bed armoured. This run was similar in output grain sizes and grain mobility to the run with one sediment pulse over a structured bed. This shows that active restructuring of the bed occurred in both of these runs. On the other end of the spectrum, the constant feed and four pulse runs were similar in their low sediment output and showed different grain mobility. This implies that the system reacts differently at a threshold frequency somewhere between two pulses per 40 h and four pulses per 40 h, which Elgueta-Astaburuaga and Hassan (2017) interpret as the relaxation time of the bed to a pulse event.

We used these experiments to test the ability of the numerical model to recreate the flume results in discharge, pulse frequency, pulse magnitude, slope, and grain sizes. In the experiment, mean grain size (D_sg) and slope (S) were measured over a central 2 m long section to avoid a bias of bed surface measurements due to inflow and outflow conditions. The measurement intervals varied between every 1 and every 20 h (depending on the pulse interval). The values for the numerical simulation are averaged over the whole reach and were recorded every 10 min during the simulation time. We achieved a best match in slope S, D_sg, D_s90, and the transport rate q_b by increasing the reference Shields stress ${\mathit{\tau }}_{\mathrm{rm}}^{*}$ in the Wilcock and Crowe formula by a factor of 2. Wilcock (2001) suggests taking the same approach of increasing the threshold shear stress to match a sediment transport calculation to field data. Other researchers use this method to calibrate models to field and flume data (Chartrand et al., 2015).

Figure 3 shows the comparison between flume measurements and the calibrated model results with the grain exchange ratio α=0.45 (Eq. 11) and the active layer thickness factor n_a=2. The sensitivity of the model to changes in α and n_a is shown in Supplement Figs. S1–S4. Due to the long interval between measurements in the flume experiments, some short-term slope responses to individual sediment pulses might be hidden (e.g., after hours 80 and 180 in Fig. 3a). The model underpredicts both the slope and the mean surface grain size (D_sg) (see Fig. 3b) in the first 60 h, while the coarse grain size fractions (D_s90) and average transport rate are over predicted for the first 30 h (see Fig. 3c and d). This might be due to imperfect initial conditions or boundary effects in the flume experiments. For the rest of the simulation both slope and D_sg show good agreement with the flume results. D₉₀ is overpredicted in the model, but this is seen as a minor issue because the simulated transport rates mainly depend on D_sg and the slope. The transport rate in the simulation lags behind the light table data, which might be due to our numerical implementation of diffusion. As the model matches the average transport rates well, we did not see the need to improve the temporal agreement.

Figure 3 Comparison of (a) slope, (b) mean surface D_sg, (c) surface D_s90, and (d) sediment transport rate between the numerical simulation and the flume experiments.

3.2 Event sequencing simulations

After obtaining a good match between the model and the flume data, we simulated alternative event sequences as described in Table 2. Figure 4a shows the adjustment of slope in runs that preserved the same sediment feed volume and had the same duration as the flume experiments (OF, MtF, and FtM), but the frequency of events is ordered differently. At the end of the simulations, all runs approach the same slope value of 0.022 m m⁻¹, which shows that the main factor determining the long-term slope is the total volume of sediment fed. The sequencing of events seems to play a role in the slope adjustment over the short term, here about 80 h after the events. On this short timescale, large pulses increase the slope quickly, while the smaller, more frequent pulses lead to a more gradual adjustment of slope. Figure 4b shows the runs where pulse phases were buffered by constant-feed phases (cMtF and cFtM), increasing the total sediment feed by 600 kg. This did not change the pattern of adjustment significantly compared to the earlier runs, as the constant-feed phases only prolonged the effect of the previous pulse phase.

The effect of event sequencing on the channel response in D_sg is shown in Fig. 4c for all runs in Table 2. The order of events only has a weak impact on patterns of adjustment in D_sg, as maximum grain size conditions are reached within 20–30 h. Afterwards no further adjustment occurs until the introduction of fine material with the next pulse lowers the surface grain size again. This means that armouring of the channel surface happens quickly in relation to the time between pulse events. The subsurface is made of the same grain size distribution as the sediment feed, so its mean size is 5.56 mm. Therefore, we can infer that an armouring ratio (D_sg∕D_subg) of about 2.2 was reached within 20 h and then increased towards 2.7 over the following 260 h. A finer, less armoured surface at the time of each supply event is followed by a coarsening of the surface as the finer grain sizes are more mobile, and thus more easily evacuated in the time without feed between pulses.

Figure 4 Comparison of slope and mean surface grain size D_sg from runs with different event sequencing (see Table 2) (σ=1.6, factor of 2 increase of ${\mathit{\tau }}_{\mathrm{rm}}^{*}$). Event sequences are (a) rearranged pulse phases to the flume experiments and (b) a setup where the pulsed phases are buffered with constant-feed phases. Panel (c) shows D_sg for the runs in both (a) and (b).

3.3 Equilibrium simulations

Our second set of simulations explored the equilibrium conditions that are reached under different supply regimes. As we will explain in the discussion, the effect of the supply regime can be constrained to changes in slope (S) and the surface grain size distribution (GSD_fluv), which in the following will be characterized by the armouring ratio between surface and subsurface mean grain size (D_sg∕D_subg). After different times in the simulations, these parameters reach a time-independent periodic adjustment that is illustrated in Fig. 5 for the last 400 h of two simulations. Figure 6a shows box plots of the distribution of slope values during the last pulse of all 40 runs with σ=1.6. The presented normalized slopes (S∕S_const) indicate how the slope for each pulse frequency compares to the constant-feed slope of runs with the same grain size distribution. Figure 6b shows the change in normalized mean slope (S_mlp∕S_const) with pulse frequency, which corresponds to the red lines in Fig. 6a and is our main indicator for the equilibrium state of the channel slope. Each line represents a different width σ of the GSD over 40 runs with increasing T_pp.

Figure 5 (a) Slope and (b) armouring ratio for the last 400 h of 2 out of 40 experiments using σ=1.6. A run with a low event frequency is shown in blue, and a run with a high event frequency is shown in red. In these example runs the mean slope in equilibrium is higher for the run with longer pulse periods.

Figure 6 (a) Distribution of the ratios of slope during the last pulse to the constant-feed slope for all 40 runs with σ=1.6. We normalized the values of the slope in the last pulse with the slope of the constant-feed run of the same grain size distribution width (σ), allowing for the comparison of the equilibrium slopes between runs with different σ. The red lines represent the normalized mean slopes (S_mlp∕S_const), which we chose as the main indicator for the equilibrium state of the channel slope. Note that data from longer pulse periods T_pp will contain more data points for the box plots, as there are more slope values recorded during the longer time between pulses (sampling every 10 min). The red crosses are outliers in the distribution of slope values, illustrating the extreme slope values during the time right after the pulse was introduced into the channel. (b) Mean slope ratios for all runs grouped by width of GSD (σ).

4 Discussion

4.1 Extension of flume results with the numerical model

The numerical model shows good agreement with the temporal response of mean surface grain size and slope from the flume experiments. As the initial bed grain size distributions were well mixed, the good match in mean surface grain size also implies a good match in the armouring ratios. A series of runs with an alternative sequencing of events showed that, while the adjustment of mean surface grain size was not sensitive to the order of the pulsed phases, the evolution of slope differed considerably. In the cases where the first pulsed phase consisted of one large magnitude event (FtM and cFtM), the slope increased quickly and the following higher-frequency, lower magnitude pulse phases did not modify the system considerably. In contrast, cases where multiple smaller event phases occurred first (MtF and cMtF), the slope increased more gradually. All runs ended at about the same slope after the 280 h simulation time, which implies that while the low frequency, large magnitude events strongly alter the channel in the short term, the sequencing of events does not play an important role in the long run. The constant-feed-buffered runs (cFtM and cMtF) show similar behaviour to their unbuffered counterparts, as the constant-feed phases preserve the bed state of the previous pulse phase. The main driver of the slope adjustment is the total sediment feed, which is consistent with findings by Blom et al. (2017).

4.2 Development of timescales from the equilibrium simulations

The definition of a sediment supply regime can be based on different aspects of sediment input into a stream, either from outside or within the channel. For simplicity, we restrict the definition of a sediment supply regime to the input of material into the fluvial system from outside the active channel. Sediment supply from storage close to the channel can be viewed as external supply if it only occurs episodically (e.g., less than yearly) in large flooding events. This view allows us to describe the sediment supply regime by a combination of frequency, magnitude, and grain size distribution of sediment supply events over a multi-event time frame (as in Benda and Dunne, 1997a). The time frame must be long enough to contain enough sediment supply events to allow the stream bed to adjust to the external forcing by changing its internal configuration of sediment storage and bed structuring. Even though a natural channel might never reach an equilibrium to a certain sediment supply regime, it might produce regular patterns of transient adjustment to the supply events.

In our case, the forcing on the system is a combination of pulse frequency (F_pulse), pulse magnitude (M_pulse), pulse grain size distribution (GSD_pulse), and water discharge (Q_w). Due to the simple geometry and the lack of bedforms in a 1-D numerical model, the fluvial reaction to the forcing is restricted to the bedload transport rate (q_b), channel slope (S), and channel grain size distribution (GSD_fluv) (see Table 3). Pulse frequency and magnitude can be combined in the virtual pulse velocity (U_pulse), with the magnitude normalized by reach width (w_r) and length (l_r):

$$\begin{array}{}\text{(12)}& U_{\mathrm{pulse}}=F_{\mathrm{pulse}}\cdot M_{\mathrm{pulse}}/(w_{\mathrm{r}}\cdot l_{\mathrm{r}}).\end{array}$$

A characteristic pulse period time (T_pp) for a sediment supply regime is the inverse of the pulse frequency:

$$\begin{array}{}\text{(13)}& T_{\mathrm{pp}}=\mathrm{1}/F_{\mathrm{pulse}}.\end{array}$$

We define a reach averaged fluvial export velocity $U_{\mathrm{fluv}}=q_{\mathrm{b}}/l_{\mathrm{r}}$ in order to compare the sediment export of material to the pulsed sediment input. Considering a multi-event time frame T_sim(≫T_pp) for each simulation, we can assume that the fluvial system will adjust to the external forcing over time, leading to an adjustment of the reach averaged fluvial sediment transport to match the external forcing of the virtual pulse velocity:

$$\begin{array}{}\text{(14)}& T_{\mathrm{sim}}\gg T_{\mathrm{pp}}:U_{\mathrm{fluv}}\approx U_{\mathrm{pulse}}.\end{array}$$

Table 3 List of forcing and reacting parameters, and timescales in our simulations with abbreviations and their dimension (T represents time and L represents length).

The time it takes to achieve an equilibrium state is highly dependent on the initial conditions of the simulations. Instead of using a process rate threshold to find an equilibrium time, we can infer a simulation-time independent relation between the supply regime (i.e., T_pp) and the state of the system in equilibrium (i.e., U_fluv, S, and GSD_fluv). This way we only have to run simulations for long enough to verify equilibrium with the respective supply regime (in our case 20 000 h), and then observe properties of the channel at the very end of the simulation, even though equilibrium might have been achieved earlier.

When comparing the system state in equilibrium between many different supply regimes, while keeping the initial fluvial reworking capacity constant (mainly geometry and Q_w), we can identify how the equilibrium conditions change under different supply regimes. Due to the fixed channel geometry in the 1-D model, only the bed slope (S) and the grain size distribution (GSD_fluv) can adjust to the change in the supply regime.

To better compare the temporal adjustment to sediment pulses of different magnitudes and frequencies, we non-dimensionalized the pulse period T_pp with a fluvial evacuation time T_fe:

$$\begin{array}{}\text{(15)}& T_{\mathrm{fe}}={\displaystyle \frac{l_{\mathrm{r}}{w}_{\mathrm{r}}{D}_{\mathrm{fg}}}{U_{\mathrm{fluv}}}}({\displaystyle \frac{D_{\mathrm{fg}}}{D_{\mathrm{ag}}}}{)}^{\mathrm{2}}.\end{array}$$

This timescale is a representation of how long it would take to remove a layer of sediment as long and as wide as the flume (l_r and w_r) with the thickness of the median feed grain size D_fg, under the average transport rate U_fluv, multiplied by an estimate of the ratio of feed grain size to armoured grain size (D_fg∕D_ag)², which can be interpreted as an inverse of the degree of potential bed armouring similar to D_s90∕D_s50 (Recking, 2012). We developed Eq. (15) visually by matching the inflection of lines in Fig. 7a to $T_{\mathrm{pp}}/T_{\mathrm{fe}}=\mathrm{1}$.

Figure 7 (a) Mean slope ratios in non-dimensional timescale. $T_{\mathrm{pp}}/T_{\mathrm{fe}}<\mathrm{1}$: model runs with a high-frequency sediment supply show similar equilibrium slopes and armouring ratios as conditions of constant sediment feed (S_mlp≈S_const). $T_{\mathrm{pp}}/T_{\mathrm{fe}}>\mathrm{1}$: if the GSD is narrow (σ<0.4) or the pulse period is not much longer than the fluvial evacuation time, we observe lower equilibrium slopes than in constant-feed runs (S_mlp<S_const). Runs with either very low frequency of supply events or wide GSD (σ≥0.4) show equilibrium slopes that are higher than the respective constant-feed runs (S_mlp>S_const). (b) Relative armouring ratio in non-dimensional timescale. $T_{\mathrm{pp}}/T_{\mathrm{ar}}>\mathrm{1}$: low frequency of supply events leads to an increase in armouring ratio compared to constant-feed runs (AR_elp>AR_const), especially for wide GSDs (σ≥0.4).

Figure 7a shows the same data as Fig. 6b, but in non-dimensionalized time T_pp∕T_fe. A ratio of $T_{\mathrm{pp}}/T_{\mathrm{fe}}<\mathrm{1}$ can be interpreted as a condition in which the pulsed input of material occurs faster than the fluvial removal of a D_fg thick theoretical layer of material under average transport conditions modified by armouring. At a ratio of $T_{\mathrm{pp}}/T_{\mathrm{fe}}>\mathrm{1}$, the sediment is fed in time steps longer than the time that the fluvial system needs to remove said theoretical layer.

Figure 7b shows the normalized armouring ratio (AR_elp∕AR_const), which was obtained by dividing the armouring ratio AR${}_{\mathrm{elp}}=D_{\mathrm{sg}}/D_{\mathrm{subg}}$ at the end of the last pulse for each simulation with the armouring ratio at the end of the corresponding constant-feed run. Similarly to the fluvial evacuation time that we used to non-dimensionalize the slope adjustment, we used an armouring time T_ar to non-dimensionalize the grain size adjustment:

$$\begin{array}{}\text{(16)}& T_{\mathrm{ar}}={\displaystyle \frac{l_{\mathrm{r}}{w}_{\mathrm{r}}{D}_{f\mathrm{90}}}{U_{\mathrm{fluv}}}}\mathrm{0.5}({\displaystyle \frac{D_{f\mathrm{90}}}{D_{\mathrm{a}\mathrm{90}}}}{)}^{\mathrm{2}}.\end{array}$$

This timescale represents how long it would take to remove a layer of sediment as long and as wide as the flume (l_r and w_r) with the thickness of the supplied D_f90, under the reach averaged transport rate $U_{\mathrm{fluv}}=q_{\mathrm{b}}/l_{\mathrm{r}}$, multiplied by an estimate of the ratio of the sediment supply D_f90 to the D_a90 of an armoured bed. We developed Eq. (16) visually by matching the inflection of lines in Fig. (7b) to $T_{\mathrm{pp}}/T_{\mathrm{ar}}=\mathrm{1}$.

4.3 Interpretation of the equilibrium simulations

The condition of $T_{\mathrm{pp}}/T_{\mathrm{fe}}=\mathrm{1}$ constitutes a threshold in the slope adjustment to pulsed sediment supply. Pulse periods shorter than the fluvial evacuation time ($T_{\mathrm{pp}}/T_{\mathrm{fe}}<\mathrm{1}$) lead to equilibrium slopes similar to the constant-feed equilibrium slopes (S_mlp≈S_const). In the case of narrow GSDs (σ<0.4), simulations with pulse periods longer than the fluvial evacuation time ($T_{\mathrm{pp}}/T_{\mathrm{fe}}>\mathrm{1}$) show an up to 20 % lower slope than in the constant-feed equivalent run (S_mlp<S_const). In contrast, runs with either very low frequency of supply events or wide GSD (σ≥0.4) show equilibrium slopes that are up to 30 % higher than the respective constant-feed runs (S_mlp > S_const). Simulations with a wider range of material (σ>1) show a drop in S_mlp∕S_const right at $T_{\mathrm{pp}}/T_{\mathrm{fe}}=\mathrm{1}$, but then an increase in S_mlp∕S_const at longer pulse distances. We interpret conditions of lower S_mlp to be less armoured, as they coincide with lower values of AR_elp∕AR_const as shown in Fig. 7b. The simulations showing an intermittent decrease of S_mlp seem to be in a state in which material can be efficiently exported from the system without having intense armouring limiting the slope adjustment. These large pulses increase the slope rapidly, leading to high shear velocities which causes high transport rates.

We interpret the cause for increasing slope and armouring ratios for long pulse periods in the following way. A longer time between pulses (T_pp>T_ar) causes more intense armouring, which shows that the channel bed is starved of sediment between pulses, as finer grain fractions are removed from the surface. This leads to the development of an armouring layer, which restricts the incision into lower deposits, initially limiting the sediment output of the system (U_fluv<U_pulse). This imbalance between input and output of material leads to an increased sediment storage over time, which due to the fixed geometry of the channel leads to an increase in slope. This can increase the shear velocity and in return leads to higher sediment output rates. This response loop between armouring and slope adjustment will continue until the sediment output matches the long-term sediment supply (U_fluv≈U_pulse). Note that due to the restricted geometry of our model, slope and grain size are the main parameters in the channel that can change in response to the sediment supply regime. It is possible that other morphological adjustments (e.g., channel width) could compensate for the transport rate disequilibrium in a similar fashion.

The non-dimensionalized presentation of the simulation results shows two distinct modes of adjustment of the fluvial system to episodic sediment supply regimes: (1) “constant-feed-like” behaviour in runs where supply events are of high frequency and low magnitude. Under this kind of forcing the equilibrium slopes and armouring ratios are similar to equilibrium conditions of runs with constant sediment feed. (2) “Pulse-dominated” behaviour occurred in runs where sediment was fed in low frequency and high magnitude events.

An interesting finding is that when T_pp>T_fe, an increase in slope only occurs in simulations that have grain size distributions wide enough to allow armouring to occur. In these cases, we could use the timescale of T_ar to determine if the armouring would be significant enough to prevent a decrease of the equilibrium slope. Even though armouring develops very quickly in both our simulations and the flume experiments, its long-term persistence in the time between sediment pulses is what governs the channel response. Hence, the grain size distribution in a series of supply events can be more important for the channel response in the long term than the frequencies and magnitudes of the individual events themselves.

It is notable that the channel response at T_pp=T_fe does not change abruptly, but the system response slowly tilts to either pulse-dominated on the one end, or constant-feed-like on the other end of the spectrum. While all constant-feed-like channels (for each σ) have very similar equilibrium properties, all pulse-dominated channels are different in both slope and armouring ratios.

4.4 Implications of the equilibrium simulations

Applying the threshold of $T_{\mathrm{pp}}/T_{\mathrm{fe}}=\mathrm{1}$ between constant-feed-like and pulse-dominated supply regimes to the flume experiments is inconclusive. The unity of $T_{\mathrm{pp}}=T_{\mathrm{fe}}\approx \mathrm{6}$ h lies between the constant-feed runs and the highest frequency runs (four pulses with T_pp=10 h), which means that we have no experimental constant-feed-like pulsed regime where T_pp<T_fe. Elgueta-Astaburuaga and Hassan (2017) found that the four-pulse phase caused a sediment transport response that was similar to the constant-feed runs, implying that T_fe would lie between 10 h and 20 h. The flume experiments were not executed long enough to reach equilibrium, which complicates the attribution of a specific channel response to a specific forcing. Besides recreating the 280 h flume experiment in the model, the 360 equilibrium simulations include four configurations that repeat the constant-feed and three pulse periods for 20 000 h. These four configurations only reached the equilibrium after about 12 000 h of simulated time, which is very long in comparison to the conditions of the flume experiments where each supply regime only lasted 40 h.

Table 4 Application of the fluvial evacuation time to the rapids reach in East Creek. The system is assumed to be in equilibrium with a matching fluvial transport rate and long-term sediment supply rate. The time of active fluvial transport is estimated to be 100 h yr⁻¹. We approximated the long-term fluvial transport rate U_fluv with three years of data from a sediment trap below the reach. We assumed that the subsurface grain size measurements reflect the average supply GSD and the surface grain size measurements represent the long-term average state of armouring in the reach. Besides the original “East Creek” data and the two “threshold” pulse frequency fits, we assumed four more scenarios with doubled armoured grain sizes or doubled fluvial transport rates to give a rough estimate of error bounds.

For all calculations: supply D_f90 =90 mm, supply D_f50 =32 mm, channel width =2.5 m, channel length = 72 m. Some frequencies calculated to match. ^a T_pp=T_fe and ^b T_pp=T_ar

The numerical model was calibrated with only one set of flume experiments. As our study mainly focusses on comparing different simulation results from the same model, the applicability of our results to other flume studies or field cases is uncertain. However, we are confident that the numerical model is an adequate tool to gain insight on the effect of episodic sediment supply on fluvial channels in a general sense. For example if there was an inaccuracy in the calculation of the shear velocity, it would affect all model runs and thus be counterbalanced by relating the changed resulting slope to the constant-feed slope (S_mlp∕S_const), and the changed armouring ratio to the constant-feed armour ratio (AR_elp∕AR_const).

We tested this by executing two additional batches of simulations with a 25 % decrease and a 25 % increase in the total mass fed respectively. If the experimental design would strongly affect the threshold between constant-feed-like and pulse-dominated conditions, we would expect these simulations to plot differently than the data in Fig. 7b. But as shown in Supplement Fig. S5 for the case of σ=1.6, the change in the slope ratio S_mlp∕S_const is relatively insensitive to the total feed volume. We also expect this to be the case when changing the grain size distributions to include large particles that are initially immobile with the applied discharge. In such simulations, the slope would increase to a point where these initially immobile grain sizes become mobile at a very high equilibrium slope. As we compare all results to the corresponding constant-feed equilibrium slope of the same grain size distribution, these conditions might collapse on the existing data as well. However, it is possible that the channel parameters would become extreme in a way that the empirically derived transport function by Wilcock and Crowe would no longer be realistically applicable.

As we developed the fluvial evacuation time T_fe and the armouring timescale T_ar purely from observations in numerical simulations, their usefulness remains to be proven in the field. If such a threshold behaviour between episodic sediment supply event frequency and the fluvial adjustment of a channel exists, it should be possible to find signatures in channel morphology, sediment storage volume, or channel slope when comparing streams subject to different pulse periods and with different fluvial transport capacities. It is possible to use Eq. (15) with information about the long-term sediment supply volume, grain size supply, and average channel dimensions to calculate T_fe and thus infer the matching threshold pulse period where T_pp=T_fe. If the long-term sediment supply occurs in more frequent events than this threshold, the system can be assumed to experience constant-feed-like sediment supply. If the supply frequency is lower, we would expect to find a morphological signature of a pulse-dominated supply regime.

To provide an example application of the developed timescales, we applied Eqs. (15) and (16) to data from East Creek, which is a small creek in the Fraser watershed close to Vancouver in BC, Canada. The “rapids” channel section of this creek was used to design the flume experiments as a 1 to 6 Froude scaled model. Table 4 shows the resulting timescales in seven different scenarios that use the reported values for sediment supply, grain size distribution, and channel dimensions from Cienciala and Hassan (2013) and Papangelakis and Hassan (2016). In East Creek, we assume that of all the supplied sediment is contributed by annual events with a magnitude that matches the annual fluvial transport. As the calculated fluvial evacuation time T_fe is 2.23 years, this scenario implies that the system would behave in a constant-feed-like manner, which would only change if the supply events were more than 2.23 years apart on average (i.e., T_pp=2.23 years), as shown in the “threshold T_fe” calculation. Due to the high uncertainty in our assumption that the measured values represent equilibrium conditions, we calculated two more scenarios with double transport rates (2×U_fluv) and double armoured grain size (2×D_a50 and D_a90). The last three calculations in Table 4 show which pulse frequency is needed to match the armouring timescale (T_ar). While we do not know if East Creek is in equilibrium with the sediment supply regime and the measurements used do not reflect long-term conditions, these calculations can still give a rough idea of whether a system is constant-feed-like in our classification of channel response to episodic supply regimes.

5 Summary and conclusions

We characterized an episodic sediment supply regime in terms of event frequency, magnitude, and supplied grain size distribution. To test the effect that different episodic sediment supply regimes can have on the morphology of a mountain stream, we developed a numerical model to recreate and extend simulations from flume experiments. The model performs well in recreating the flume experiments in both slope and grain size distributions (GSD), which are the two variables that represent morphological adjustment in our model. Channel width is fixed and bedforms are assumed to be absent, even though bedforms did occur in the flume experiments.

To understand the extent to which event succession plays a role in the flume experiments, we simulated alternative pulse successions of large-to-small events (i.e., infrequent-to-frequent) and small-to-large events (i.e., frequent-to-infrequent), while keeping the total sediment volume feed the same. These simulations show that different pulse frequency sequences have no strong effect on the long-term slope and GSD of the bed surface. In the short term large pulse events can dominate the channel response causing an abrupt increase in slope, while the effect of subsequent smaller events is subdued as the channel is still adjusted to the large pulse. If smaller events dominate at first, the channel adjusts more gradually.

In our second set of simulations, we imposed different episodic sediment supply regimes with the same total sediment supply volume on the same initial channel geometries with constant discharge. While being kept constant within a run, the episodic supply regimes differed in event frequencies, magnitudes, and GSD. We simulated 40 different event frequencies for which the sum of event magnitudes matched an overall equal total sediment supply. All 40 pulse configurations were calculated for 9 GSD that differed in the width of the distribution σ around the same geometric mean grain size. The channels adjusted to the episodic sediment supply until they reached an equilibrium state in which each successive pulse led to the same slope and grain size adjustment. We compared this state between runs and found a distinctive regime change when the time between pulses (T_pp) became lower than a fluvial evacuation time (T_fe), which we developed as a measure of the time it takes to remove a D_s50 thick layer from the channel surface under average transport conditions, modified by a measure of potential armouring (see Eq. 15).

The condition of T_pp<T_fe causes a constant-feed-like sediment supply regime, as the model runs show similar slopes and surface grain size distributions as constant-feed runs of the same GSD. When T_pp>T_fe the sediment supply regime becomes pulse-dominated. Under these conditions, we observed a lower relative slope in cases where the GSD is narrow (σ < 0.4), as the long time between pulses in combination with a low armouring potential allows more erosion in the reach, ultimately lowering the equilibrium slope. If the GSD is wide enough to allow armouring (σ≥ 0.4), a stronger armouring layer can develop during the periods of selective transport of smaller grain sizes and bedload starvation. This limits the minimum slope and increases sediment storage (and thus slope) in the long term.

The application of the episodic supply regime classification to data from East Creek shows that the threshold to a pulse-dominated regime lies at the fluvial evacuation time of roughly 2.2 years. This creek probably receives sediment at a lower interval, which indicates a pulse-dominated regime. The armouring timescale lies around 3.5 years, indicating that if the long-term sediment supply was introduced over event frequencies between 2.2 and 3.5 years, it would be removed most efficiently and result in a lower slope.

Steeper channels than East Creek could show both a lower fluvial evacuation time (due to higher slope, smaller channel area) and a lower pulse frequency (more landslide dominated), which could make these channels more likely to be pulse-dominated. Further study of field cases is needed to strengthen the case for our classification of channel response types to episodic supply regimes. In natural rivers, there are further modes of adjustment that the system can undergo after receiving sediment pulses, for example changes in bed forms or the storage of excess material in sediment bodies along the channel. Still, the condition when a channel is receiving more material per pulse than what can be exported in the same time frame (i.e., the ratio of T_pp∕T_fe is above 1), should be observable in natural rivers as irregularities in channel long profiles due to increased sediment storage. In our model, we only supplied grain sizes that were transportable by the imposed flow conditions. In field streams it can be that the biggest clasts (e.g., boulders) are only transportable by extreme flow events, which would further increase the slope of reaches with high sediment supply.