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.

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' has been later elaborated by Thomas (2001) and Brunsden (2001) but afterwards little attempt has been made in fluvial geomorphology to address this issue in practice. 5 Bull (1991) applies the theory of Brunsden and Thornes (1979) to the impact of a hypothetical temporal succession of disturbance events on sediment storage in a stream-bed. Episodic disturbance events of sediment contribution and sediment export result in an increase or decrease of the stream-bed elevation, respectively. This reaction is not instantaneous, but delayed by a reaction time in which no adjustment takes place. The system processes the disturbance over the relaxation time, which together with the reaction time constitutes the total response time. If there is no further disturbance, the system status remains 10 unchanged over a time of persistence. Sequenced disturbance events might lead to overlaps of reaction times and relaxation times and thus create a complex pattern of system changes. Due to the complexity of geomorphic systems, it might be impossible to discern these temporal phases of system reaction in field studies. Bull's concept is based on a system's trend towards a dynamic equilibrium between the forcing by, and the reworking of, disturbances. If the forcing on the system is constant over time, we can expect the system to eventually reach this dynamic equilibrium. 15 Howard (1982) concluded that if episodic inputs occur with a frequency that matches the inverse of relaxation time, the output of the system will remain in a constant equilibrium with the average value of the forcing. Numerical modelling can be used to incorporate complex properties of the system (e.g. armouring, dynamic changes in grain size composition), and thus to allow for non-linear relations between input and output parameters. This enables us to explore equilibrium conditions in a quasi-grey box system, where the exact relationship that produces equilibrium between input and output must not be 20 understood in detail, but still allows for larger scale understanding of the system's response.
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 on 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 25 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 30 using a one-dimensional morphodynamical 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 conducted to examine the impacts of episodic sediment supply on bed surface evolution and channel adjustment of a gravel bed stream (von Flotow, 2013;Elgueta, 2014;Elgueta-Astaburuaga and Hassan, 2017;Elgueta-Astaburuaga, 2018). Although these experiments provide detailed information on channel adjustment to changes 35 in the sediment supply regime, they are limited in terms of number of experiments and range of scenarios that can 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 are: (1) Does the sequencing of supply events play a role in the reaction of a gravel-bed stream, when several events of specified 5 magnitudes occur in a different ordering?
(2) How will different combinations of episodic sediment supply, obtained by varying their magnitude and frequency, impact channel evolution of a gravel-bed stream? And can the numerical model recreate the channel response that was observed in flume experiments of similar scope?

10
We developed the one-dimensional morphodynamical model BESMo (Bedload Scenario Model) to calculate capacity based sediment transport under different sediment supply regimes. We parameterized the model to match the flume experiments as closely as possible and used measurements of sediment transport rate, surface grain size distribution, and slope to validate the model's performance. Matching our research questions, we then conducted two types of simulations: (1) in 'event sequencing experiments', we simulated different permutations of events to understand the role event succession plays in long term channel 15 response; (2) in 'equilibrium experiments', we used the same model setup and imposed different, but within each run regular supply event frequencies for durations until we achieved a recurrent pattern in slope and grain size adjustment, allowing to identify how the channel adjusts to the supply regime in the long-term.

Model setup
The structure of the model is similar to other models designed to reproduce and interpret data from flume experiments (e.g. 20 Cui and Parker, 2005;Wong and Parker, 2006;Ferrer-Boix and Hassan, 2014). 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 accordance to modelling scenarios.
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. 25 (2006), implementing a threshold Froude number Fr to switch conditions between supercritical and subcritical flow: The Froude number F r is calculated as a function of water depth Q w , gravity g, channel width w r and water depth h (eq. 1). slope S f and bed slope S 0 (eq. 2a): Water depth under supercritical flow conditions is calculated locally assuming steady uniform flow using the Manning-Strickler formulation (eq. 2b & 3), where α r is a coefficient of 8.1 (Parker, 1991) and roughness height k s is calculated in eq. 4 choosing 5 n k = 2 and D s90 as the grain size for which 90% of the surface is finer: This approach simplifies friction slope S f to be equal to bed slope S 0 with bed elevation η b at downstream position x (eq. 5): If the solution of water depth with eq. 2 is 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 5 of at least one channel width. Boundary shear stress τ b is then calculated with the depth-slope product: with ρ being the water density. τ b is then converted to the shear velocity u * , which is used in the sediment routing component: The volumetric unit bedload transport rate per size class q bxi is calculated using the sediment transport function provided by 10 Wilcock and Crowe (2003). The change in bed elevation ∂η b in for each node x per time step t follows from the Exner equation of mass conservation: with λ as bed porosity. The volumetric bedload transport rate per unit width is given as q bx and is calculated per grain size class i in the sediment mixture The model incorporates subsurface stratigraphy using the active layer concept (Parker, 2008), which gives a modified Exner equation: with L a being the active layer thickness, F i the surface frequency of the ith grain size class, f Ii the ith grain size class 20 proportion exchanged between the surface and the subsurface, and q bxi = p bi q bx the volumetric unit bedload transport rate of the ith grain size class with p bi being the ith fraction of the bedload transport rate. The active layer thickness is calculated as L a = n a D s90 , with the parameter n a = 2, representing the scale of bed fluctuations.
The exchange of material between the layers 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): with f i as the fraction of the ith grain size class in the subsurface and α being a coefficient of 0.45. The vertical stratigraphy is stored in layers of a specified height 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 active layer thickness factor n a and active layer exchange ratio α have on the model results, we executed sensitivity runs (Results shown in supplement Figures S1-S4).

5
The Exner equation (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 and allows the calculation of bed elevation as a function of space and time . The model needs an initial bed profile and an initial value of the surface grain size to be solvable. 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 bedload transport function is used to calculated transport rates for each 10 channel cross section. During the model run, changes of the sediment transport rate are depending on changes in the sediment supply, channel slope, and surface grain size distribution.
Using the described model, we simulated a 12 m long and 1 m wide channel in 13 downstream nodes each spaced 1 m apart.
This channel is chosen to match the flume experiments by Elgueta-Astaburuaga and Hassan (2017). The model was set to calculate sediment transport for all nodes in time steps of 10 s. All simulations use a constant water discharge of 0.065 m 3 /s, a 15 geometric mean grain size of 5.64 mm for both initial bed and sediment feed (full distribution shown in Fig. 2a), and an initial slope of 0.022 m/m, also matching the parameters from the flume experiments.

Model validation
We Flow and sediment data are summarized in Table 1. The bed was fixed in the first 1 m downstream of the flume headbox 25 with stones equivalent to about D84 of the experimental bed material. In the remainder of the flume the bed consisted initially of 0.1 m of loose material with particle sizes ranging from 0.5 mm to 64 mm with a D 50 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 bed load texture. Measurements of the water surface elevation were conducted 30 throughout the experiment using a mechanical point gauge with precision 0.001 m. Photos were used to determine bed surface texture 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 of exported material and the number of stones exiting the flume were measured with a camera and a light table at the outlet of flume.
The experiment consisted of seven 40 hour 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 can be studied 5 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 1 pulse, in 2 pulses, or in 4 pulses. 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 1 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 4 pulse runs were similar in their low sediment output and showed different 10 grain mobility. This implies that the system reacts differently at a threshold frequency somewhere between 2 pulses / 40 h and 4 pulses / 40 h, which is interpreted to represent the relaxation time of the bed to a pulse event (Elgueta-Astaburuaga and Hassan, 2017).
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 g and slope S were measured over a 15 central 2 m long section in intervals varying between every 1 and 20 hours (depending on the pulse interval), while the values for the numerical simulation are averaged over the whole reach and recorded every 10 minutes in simulation time. We achieved a best match in slope S, D g , D 90 , and the transport rate by increasing the critical reference Shields stress τ * rm in the Wilcock & Crowe formula by a factor of 2 (see Fig. 3). In some cases the long interval between measurements in the flume experiments might have hidden some short-term slope response to individual sediment pulses (e.g. after hours 80 and 180).  Table 2. Sequencing of events in runs that were simulated as permutations of the original flume experiment. Each of the seven periods was 40 hours long and each run lasted 280 hours in total. There was no sediment input during the no feed runs (0). Within all other period-types, 300kg of sediment was fed over 40 hours, either constantly (C) or in pulses (1, 2, or 4 events). Besides recreating the original flume sequence (OF), we simulated two runs where the pulsed events either occur in order from many pulses to few (MtF) or from few to many (FtM). To explore if the system can rebound from the impact of a certain pulse phase during a constant feed phase, we created two additional runs where this is the case (cMtF and cFtM).

Run name
Event sequence periods Mass fed total [kg] OF Original Flume 0 C 1 4 2 C 0 1500 MtF Many to Few 0 C 4 2 1 C 0 1500 FtM Few to Many 0 C 1 2 4 C 0 1500 cMtF C-Buffered: Many to Few C 4 C 2 C 1 C 2100 cFtM C-Buffered: Few to Many C 1 C 2 C 4 C 2100

Event sequencing experiments
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 (4, 2 or 1 pulses over 40 hours), but the pulse ordering is either from 'Few to Many (FtM)' (i.e. 1, then 2, then 4 pulses) or from 'Many to Few (MtF)' (i.e. 4, then 2, then 1 pulses) per 40 h phase. To allow the system 5 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).

Equilibrium experiments
In our final set of experiments, we kept the frequency and magnitude of pulse events constant to achieve equilibrium slope and grain size conditions. We define equilibrium following Ahnert (1994) meaning that the system reaches a long-term balance 10 between erosional and depositional forces, even though there are periodic changes in the short term (dynamic equilibrium).
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 at first increase the channel slope.
After the supply of sediment ends and material is removed from the channel, the slope will decrease, and the surface grain size 15 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.
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 9 different grain size distributions (Fig. 2a). All distributions 5 have the same mean grain size of 5.64 mm, but differ in the width of the distribution by the standard deviation σ, that 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. 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 10 over different pulse frequencies, with four of the combinations matching the flume experiments. Each pulse was 10 minutes in length. The lowest frequency was chosen to be one pulse every 400 hours (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 a whole number of cycles summed up to 20,000 h. In total we executed 360 simulations, 9 different σ with 40 frequencies each.

15
A basic relation can be drawn between the external forcing on the stream, given by sediment supply and water discharge, and the internal reaction of the stream by adjusting the throughput, storage, or entrainment of sediment. As described earlier, we expect that both the forcing and the reacting processes, if sustained long enough, to achieve a dynamic equilibrium, so that the long-term sediment input equals the long-term sediment output of the channel. In a pulsed supply regime, we expect this relation to be cyclic in the short term, exhibiting a recurrent channel adjustment with the input and the reworking of each individual sediment pulse, as described by Howard (1982). In field studies it might be impossible to quantify these temporal phases, as for example flood hydrographs introduce an additional forcing that renders the attribution of a currently dominating 5 channel reaction to any given input complex. Additionally, in the field the observation time is often shorter than recurrence intervals, requiring historical data that often is unavailable. The flume studies discussed above tried to isolate the causes of any given reaction, but changes in sediment transport and surface structuring happen hand in hand. Furthermore, the system reacts to the external forcing in a complex way that is highly dependent on the initial bed state. Still, our simplified view of process interaction includes enough elements to help understand a complex system by comparing how a model performs in 10 comparison to real world data. 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 15 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 reach temporary or local equilibria, and produce regular patterns of transient adjustment to the supply events. 20 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 1D 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 f luv (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 : A characteristic pulse period time T pp for a sediment supply regime is the inverse of the pulse frequency: We define a reach averaged fluvial export velocity U f luv = q b /l 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 30 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: It is challenging to define the exact point in time the system reaches equilibrium, as multiple aspects of the system response might reach an equilibrium at different times. For example a transport rate equilibrium might be reached before an equilibrium in slope. More importantly, 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 f luv , S, and GSD f luv ). This way 5 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 10 supply regimes. Due to the fixed channel geometry in the 1D model, only the bed slope S and the grain size distribution GSD f luv can adjust to the change in the supply regime.

Results
In the flume D g and D 90 were measured in a 2 m wide middle section, while slope, D g , and D 90 in the simulation are averaged over the full 12 m length. The simulation transport rate was recorded as the sediment output of the lowest node, while the values 15 of the flume experiments were measured using a light table at the flume outlet. We achieved the best match (shown in Fig. 3) with a grain size distribution of width σ = 1.6, a factor 2 increase in τ * rm in the Wilcock & Crowe formula, and a grain exchange ratio α of 0.45 (eq. 11). Supplement Figures S1-S4 show the sensitivity of the model to changes in values of α and n a .

Event sequencing experiments
After obtaining a good match between the model and the flume data, we simulated alternative event sequences as described 20 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 on the slope adjustment over the short term, here about 80 hours after the events. On this shorter time scale, 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 The effect of event sequencing on the channel response in D g is shown in Fig. 4c for all runs in Table 2. The ordering of events has only a weak impact on patterns of adjustment in D g , as maximum grain size conditions are reached within 20-30 hours. Afterwards no further adjustment occurs until the introduction of fine material with the next pulse lowers the surface 5 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 g,surf /D g,subsurf of about 2.2 was reached within 20 hours and then increased towards 2.7 over the following 260 hours. 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.

Equilibrium experiments
Our second set of simulations explored the equilibrium conditions that are reached under different supply regimes. As mentioned earlier, the effect of the supply regime can be constrained to changes in slope S and the surface grain size distribution GSD f luv , which in the following will be characterized by the armouring ratio D g,surf /D g,subsurf . After different times in the simulations, these parameters reach a time-independent periodic adjustment that is illustrated in Fig. 5.
15 Figure 6a shows box plots of the distribution of slope values during the last pulse of all 40 runs with σ = 1.6. We normalized the values of slope with the slope of the constant-feed run of the same grain size distribution width σ, allowing to compare the equilibrium slopes between runs with different σ. The red lines of the box plots represent the normalized mean slopes (S mlp /S const ), which we choose 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 20 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. Figure 6b shows the change in normalized mean slope. Each line represents a different wideness σ of the GSD over 40 runs with increasing T pp .
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 f e (eq. 15). The latter is a representation of how long it would take to 25 remove a layer of sediment as long and wide as the flume (l r and w r ) with the thickness of D g , under the average transport rate U f luv , multiplied by an estimate of the ratio of feed grain size to armoured grain size (D g,supply /D g,armoured ) 2 , which was matched visually and can be interpreted as an inverse of the degree of potential bed armoring similar to D 90 /D 50 (Recking, 2012).
T f e = l r w r D g,supply U f luv D g,supply D g,armoured 2 (15) 30 Figure 6c shows the same data as Fig. 6b, but in non-dimensionalized time T pp /T f e . A ratio of T pp /T f e < 1 can be interpreted as a condition in which the pulsed input of material occurs faster than the fluvial removal of a D g thick theoretical  layer of material under average transport conditions modified by armouring. At a ratio of T pp /T f e > 1, the sediment is fed in time steps longer than the time that the fluvial system needs to remove said theoretical layer. The condition of T pp /T f e = 1 constitutes a threshold in the slope adjustment to pulsed sediment supply. Pulse periods shorter than the fluvial evacuation time (T pp /T f e < 1) lead to equilibrium slopes similar to the constant feed equilibrium slopes (S mlp ≈ S const ). In simulations with a narrow GSD (σ < 0.4), we observe a decrease in S mlp when the pulse period is long (T pp /T f e > 1). In these cases the time 5 between pulses allows more erosion so that less material is stored in the channel and armoring is nearly non-existent. The mean equilibrium slope increases (S mlp > S const ) for runs with wider GSD (σ > 0.7) and long pulse distances. In these simulations,

Time (hours) (b)
Pulse period = 20 hours mean Pulse period = 40 hours mean Figure 5. (a) Slope and (b) armouring ratio for the last 400 h of two out of 40 experiments using σ = 1.6. A run with low event frequency is shown in blue, a run with high event frequency in red. In these example runs the mean slope in equilibrium is higher for the run with longer pulse periods. the longer time between pulses allows an armouring layer to develop which hinders further erosion into material left by earlier pulses. Figure 6d shows the normalized armouring ratio (AR elp /AR const ), which was obtained by dividing the armouring ratio AR elp = D g,surf /D g,subsurf 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, 5 we used an armouring time T ar to non-dimensionalize the grain size adjustment (eq. 16). This visually identified timescale represents how long it would take to remove a layer of sediment as long and wide as the flume (l r and w r ) with the thickness of the supplied D 90 , under the reach averaged transport rate U f luv = q b /l r , multiplied by an estimate of the ratio of the sediment supply D 90 to the D 90 of an armoured bed.
The normalized armouring ratio increases for all runs where the pulse period is longer than the armouring time (T pp /T ar > 1). The increase in armouring ratio is highest for σ = 1 which could mean that the armouring layer does not fully develop in the constant feed runs for these GSDs. Note that even the runs with high pulse frequencies develop an armouring layer, but not necessarily in a way that is restricting erosion to the point where the adjustment of the channel slope is affected. Low frequency of supply events leads to an increase in armouring ratio compared to constant feed runs (AR elp > ARconst), especially for wide GSD's (σ ≥ 0.4).

Comparison to flume results and event sequencing
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 5 adjustment of mean surface grain size was not sensitive to the ordering of the pulsed phases, the evolution of slope differed considerably. In 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 10 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.

Equilibrium experiments
In our second set of experiments we focused on the equilibrium slope and grain size conditions in simulations under different (2) 'Pulse-dominated' behaviour occurred in runs where sediment was fed in low-frequency and high-magnitude events. In case of a narrow GSD (σ < 0.4) or 20 a pulse period that is not much longer than the fluvial evacuation time (1 < T pp /T f e < 20), we observe equilibrium slopes that are up to 20% lower 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 up to 30% higher than the respective constant feed runs (S mlp > S const ). In this case, the supply regime caused an increase in both equilibrium slope and armouring ratio in comparison to constant feed runs that had the same GSD and total sediment feed volume. The concurrence of more intense armouring with higher slopes 25 implies that this 'pulse-dominated' type of channel response occurs due to the long time between pulses, in which the channel bed is starved of sediment supply and 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 f luv < 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, which increases the shear velocity and in return 30 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 f luv ≈ 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.
We defined a threshold between constant-feed-like and pulse-dominated episodic sediment supply regimes in terms of a pulse period T pp that equals the fluvial evacuation time T f e . Applying this concept to the flume experiments is, however, inconclusive. The unity of T pp = T f e ≈ 6 h lies between the constant feed runs and the highest frequency runs (4 pulses with The timescale T f e corresponds to the relaxation time as discussed by Brunsden and Thornes (1979), Howard (1982), and 15 Bull (1991), and is a way to quantify that a system is more sensitive to events that occur with lower return intervals than the relaxation time of the system. The ratio between the 'fluvial evacuation time' T f e and the 'pulse period' T pp constitutes a threshold that can help to quantify a system response time to episodic disturbances.
An interesting finding is that when T pp > T f e , an increase in slope only happens in simulations that had grain size distributions wide enough to allow armouring to occur. In these cases, we can then use the timescale of T ar to determine if the 20 armouring would be significant enough to prevent a lowering 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 even 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 f e does not change abruptly, but the system response slowly tilts to either 25 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.

Implications for river morphology
As we developed the fluvial evacuation time T f e and the armouring timescale T ar purely from observations in numerical 30 simulations, their usefullness 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 celerities. If we get information about the long-term sediment supply volume, grain size supply, and average dimensions for a channel, we can use eq. 15 to calculate T f e and thus infer the matching threshold pulse period where T pp = T f e . 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.
As an example, we used eq. 15 and eq. 16 to calculate the corresponding timescales for East Creek in British Columbia, 5 Canada, using data from the same 'rapids' section after which the flume simulations were designed. Table 4 shows the resulting timescales in seven different scenarios. 'East Creek' and 'Threshold T f e ' 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 all supplied sediment is contributed in annual events with a magnitude that matches the annual fluvial transport.
As the calculated fluvial evacuation time T f e is 2.23 years, this scenario implies that the system would behave constant-feed-10 like, which would only change if the supply events are more than 2.23 years apart on average (i.e. T pp = 2.23 years), as shown in the 'Threshold T f e ' calculation. Due to the high uncertainty in our assumption that the measured values represent equilibrium conditions, we calculated the two more scenarios with double transport rates (2x U f luv ) and double armoured grain size (2x ). 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 used measurements do 15 not reflect long-term conditions, these calculations still can give a rough idea whether a system is constant-feed-like in our classification of channel response to episodic supply regimes.

Conclusions
We characterized an episodic sediment supply regime in terms of event frequency, magnitude, and grain size distribution. To test the effect that different episodic sediment supply regimes can have on the morphology of a mountain stream, we developed 20 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 to what extent 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-25 infrequent), while keeping the total sediment volume feed the same. These simulations showed that different pulse frequency sequences have no strong effect on the long-term slope and GSD of the system. In the short-term large pulse events can dominate the channel response causing an abrupt increase in slope, while the effect of following 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 30 supply volume on initially equal 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 Table 4. Application of the fluvial evacuation time to the rapids reach in East Creek. The system is assumed to be in equilibrium with matching fluvial transport rate and long term sediment supply rate. The time of active fluvial transport is estimated to be 100 hours/year.
We approximated the long-term fluvial transport rate U f luv with three years of data from a sediment trap below the reach. We assumed the subsurface grain size measurements to reflect the average supply GSD and the surface grain size measurements to represent the long-term average state of armouring in the reach. Besides the original 'East Creek' data and 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 9 GSD that differed in the wideness of the distribution σ around the same geometric mean grain size. The channels adjusted to the episodic sediment supply until they reached a state in which each successive pulse led to the same slope and grain size adjustment. We compared this equilibrium state between runs and found a distinctive regime change when the time between pulses T pp got lower than a fluvial evacuation time T f e , which we developed as a measure of the time it takes to remove a D 50 thick layer from the channel surface under average transport conditions, modified by a measure of potential armouring 5 (see eq. 15). We term the condition T pp < T f e as 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. We describe the conditions where T pp > T f e as a pulse-dominated sediment supply regime. In these, 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 10 can develop during the periods of seletive transport of smaller grain sizes and bed load 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 pulsedominated regime lies at the fluvial evacuation time of roughly 2.2 years. This creek probably receives sediment at a lower pulse frequency, which indicates a pulse-dominated regime. The armouring timescale lies around 3. 5  the long-term sediment supply would be introduced over event frequencies between 2.2 and 3.5 years, it would be removed most efficiently and result in a smaller slope.
Channels higher up in the mountains than East Creek can show both a lower fluvial evacuation time (due to higher slope, smaller channel area) and a lower pulse frequency (landslide dominated), which should increase the likelihood of pulsedominated channels. Further study of field cases is needed to strengthen the case for our classification of channel response 5 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 behaviour we observe in low frequency sediment supply regimes, 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 f e 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 10 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.