A reduced-complexity model for sediment transport and step-pool morphology

A new particle-based reduced-complexity model to simulate sediment transport and channel morphology in steep streams in presented. The model CAST (Cellular Automaton Sediment Transport) contains phenomenological parameterizations, deterministic or stochastic, of sediment supply, bed load transport, and particle entrainment and deposition in a cellular-automaton space with uniform grain size. The model reproduces a realistic bed morphology and typical fluctuations in transport rates observed in steep channels. Particle hop distances, from entrainment to deposition, are well fitted by exponential distributions, in agreement with field data. The effect of stochasticity in both the entrainment and the input rate is shown. A stochastic parameterization of the entrainment is essential to create and maintain a realistic channel morphology, while the intermittent transport of grains in CAST shreds the input signal and its stochastic variability. A jamming routine has been added to CAST to simulate the grain–grain and grain–bed interactions that lead to particle jamming and step formation in a step-pool stream. The results show that jamming is effective in generating steps in unsteady conditions. Steps are created during high-flow periods and they survive during low flows only in sediment-starved conditions, in agreement with the jammed-state hypothesis of Church and Zimmermann (2007). Reduced-complexity models like CAST give new insights into the dynamics of complex phenomena such as sediment transport and bedform stability and are a useful complement to fully physically based models to test research hypotheses.

In this paper we present a new reduced-complexity stochastic model for step-pool streams based on grain-grain interactions: CAST (Cellular Automaton Sediment Transport). CAST simulates a generic fluvial channel on a cellular-automaton domain, 5 where the bed is formed by an arrangement of particles like in a sandpile model (e.g. Bak et al., 1988;Kadanoff et al., 1989).
The basic processes of bed load transport, particle entrainment and deposition are treated at the grain scale taking advantage of analogies between bed load transport and granular phenomena (e.g. Church, 2009, 2011;Houssais et al., 2015).
The stochastic framework of CAST has two main reasons. First, the goal of the model is not to predict deterministically the morphology of a specific river reach but rather to capture feedbacks related to its dynamics and to test research hypotheses 10 on step formation and stability. Second, both bed stability (Zimmermann et al., 2010) and bed load transport (Einstein, 1937(Einstein, , 1950 have been recognized to be stochastic processes, and recent approaches to sediment transport have successfully followed this framework (e.g. Turowski, 2010;Furbish et al., 2012;Heyman et al., 2013;Ancey and Heymann, 2014;Armanini et al., 2015).
The paper objectives are: (a) to present a new reduced-complexity model that simulates bed load transport and channel 15 morphology at the grain scale and to test the effect of different parameters and stochastic forcing on the model outcomes; (b) to explore the effect of jamming on sediment transport and step formation, in comparison with the framework of the the jammed-state hypothesis of Church and Zimmermann (2007). The paper is organized as follows. First, the model rationale and parameters are presented. Second, the effect of different parameter sets on the model outcomes are explored in steady-state simulations, and the effect of the stochasticity on different variables is shown. Then, the process of jamming is parameterized 20 and its role on step formation and stability is explored in unsteady simulations. Finally the results are discussed and compared to the jammed-state hypothesis.

Model Rationale
CAST operates in 2-D cellular-automaton space, which is a rectangular grid of constant length (X) and width (Y ) corresponding to a generic river reach (see Fig. 1). The domain is discretized such that all the dimensional quantities are expressed as 25 multipliers of particle size d. The model developed in this paper works with uniform-size particles, so d is also equal to the dimension of a cell. For example, a simulation domain having X = 300d and Y = 20d represents a river reach with an average length equal to 300 median diameters and an average width equal to 20 median diameters. Particles in the model domain can be either on the bed or in motion, i.e. part of the bed matrix or of the transport matrix ( Fig. 1). The bed matrix Z is composed of particles piled one above the other like in a sandpile model (e.g. Bak et al., 1988). 30 The local bed elevation at a generic location (i, j), where i[1 : X] is the index for the longitudinal coordinate and j[1 : Y ] is that of the transversal coordinate, is given as the total number of particles Z i,j . Particles can leave the bed matrix as a result of entrainment and can enter the bed matrix as a result of deposition. In the case of entrainment and deposition the local value of elevation Z i,j is reduced or increased by one grain unit d, respectively: Particles in motion are allocated to the transport matrix T R, which consists of two layers, and they move as bed load, i.e.
they are in contact with the bed and interact with it. They also interact with each other by collisions. Particles move with a 5 constant velocity, one cross-section downstream for every time step. In this way particle velocity v p , particle size d and time step ∆t are connected: Although in simulation ∆t is a unitless time, Eq. 2 together with the grain and domain size give it a physical meaning connected to particle velocity. 10 Particles enter the system with a specified input rate I R which is the number of particles entering the system at the upstream boundary for every time step, and they leave the system as output rate O R at the downstream boundary. In analyzing the spatial output of the model we consider only a reduced part of the domain which we will call hereafter the control volume, excluding the first 10 cross-sections upstream and the last 10 downstream, to avoid the influence of the upstream and downstream boundary conditions (see next sections). The first parameter of CAST is the particle input rate I R or the specific input rate i R , defined as the total input rate I R divided by the channel width Y ; i R can assume values in the interval [0 : 1] because no more than one particle can enter a single cell at the upstream boundary in ∆t. The supply of particles to the system can be treated as constant input by specifying a value of 20 i R for the entire simulation or as variable random input by specifying a mean value i R with a temporally-variable term i R (t), uniformly distributed around i R . In CAST , i R represents a generic amount of sediment which is delivered to the channel from all the possible sources (alluvial transport, colluvium activity, bank erosion, etc) rather than a specific transport rate in a given cross-section. The actual input to the system can be considered to be the transport rate measured in the first cross-section of the control volume.

Sediment Transport
Particles are transported as bed load along the channel with a constant velocity (see Equation 2). A particle can join the transport matrix T R when it enters the system as input or once it has been entrained from the bed. A particle can also leave the transport matrix T R when it moves beyond the last cross-section, becoming part of the sediment output O R , or when it deposits on the bed surface. The maximum number of particles being transported from one cross-section to the next one is equal to 2Y , i.e. 2 times the channel width, because the T R matrix has two layers.
Particles move preferentially directly to the downstream cell (90% probability), with a small chance for lateral displacements (5% to the left and 5% to the right). This aims to represent a grain having a transport vector aligned with the dominant flow 5 direction with limited diversion. In absence of observations, we constrain the probability for lateral movements in a reasonable interval of 10%. The model is not sensitive to lateral dispersion, at least in the parameter space we have tested. Along its path, a particle can collide with another particle in transport or collide with one of the two boundaries (left and right banks). In both cases the collision leads to the loss of momentum and cessation of motion and the particle deposits on the bed (see section 2.1.4). When a particle deposits on the bed, it changes the local roughness but without directly displacing other particles, i.e. 10 the CAST does not account for collective entrainment as described by Ancey and Heymann (2014).
Sediment flux in the model, q S , is computed at the total number of particles in T R in the control volume divided by the domain size, i.e. Y · (X − 20). This specific rate q S (t) is computed for every time step.

Particle Entrainment
The key process in CAST is the particle entrainment which is considered to be dependent on the local bed topography and on 15 the flow conditions. The degree of exposure of particles on the bed has been shown to strongly influence sediment entrainment and transport especially in steep streams (e.g Kirchner et al., 1990;Malmaeus and Hassan, 2002;Yager et al., 2012;Prancevic and Lamb, 2015). Moreover many feedbacks exist between bed roughness, flow resistance and particle mobility and transport (e.g. Recking et al., 2008;Wilcox et al., 2011) which makes it reasonable to consider particle entrainment as a stochastic process.

20
The effect of local topography on entrainment is accounted for by calculating the local relative exposure R of a particle on the bed. For a generic cell (i, j) in the domain, the relative exposure R i,j is given by the difference between the elevation of the cell Z i,j and the average elevation of the neighboring cells along the flow direction (the 2 in the same cross-section and the 3 downstream): In the case of a cell located close to one of the banks (i.e. when j = 1 or j = Y ), R i,j is evaluated considering the cell in the same cross-section and the 2 cells downstream.
Entrainment is based on R exceeding a threshold R * for entrainment. The probability of entrainment p E is then defined as CAST can model the entrainment process as deterministic or stochastic (see Fig. 2). In the deterministic case, the threshold is a constant R * = E and the probability of entrainment is: In the stochastic case, the threshold R * is modelled as a random variable with a logistic probability density function f (R * ) 5 and a cumulative distribution function F (R * ): The distribution has a mean µ R * = E and a variance σ 2 R * = π 2 S 2 3 . In this way the entrainment probability p E of a cell with relative exposure R depends on two parameters: the mean value of the threshold distribution E and the variability in the threshold R * proportional to S: Figure 2 shows the deterministic and stochastic parameterizations of entrainment.
Conceptually, the value of E is inversely related to the magnitude of the flow. Large E means low probability of entrainment, typical of low flow conditions and low shear stress, while small E means high p E for the same relative exposure values and so high flow conditions and high shear stress. 15

Particle Deposition
Particles in transport (i.e. belonging to the T R matrix) can be deposited in three cases: (a) when they collide with another particle in motion; in this cases both particles involved in the collision events are deposited, each one in its present location; (b) when they collide with one of the two channel banks; (c) because of their interaction with the bed surface.
The relation between particle deposition and bed surface is modeled using the relative exposure matrix R: particles in motion 20 deposit in areas of local depressions, i.e. cells with R < 0. The deposition process is treated as deterministic, with a threshold function having a fixed value of −0.5 below which the probability of deposition p dep = 1. This simplification allows on one hand to avoid redundant parameters poorly connected to physical processes, and on the other hand to shift the variability in sediment transport to the entrainment process which has a more straightforward relation with hydraulics and local channel bed topography. Case (a) and (b) are much less common than (c) since they both require lateral movement and the presence of an 25 obstacle (i.e. another particle or the channel banks).

Boundary and Initial Conditions
CAST needs one boundary condition for the lateral banks and one for the downstream boundary at the channel outlet. The boundary condition for the banks is deposition when a moving particle in the T R layer collides with one of the two lateral boundaries. The boundary condition for the last cross-section at the downstream boundary is given by fixing its elevation Z(X, j) = 0. This is equivalent to a control section with a fixed elevation (e.g. a check dam or a weir) somewhere downstream. 5 To minimize the effect of this boundary condition on the model outcomes, all spatial variables are computed only over the control volume, which is a reduced portion of the entire channel (see Fig. 1).
In order to avoid long simulation times required to fill the channel with particles, we start every simulation from an initial slope, slightly less than the equilibrium slope, with random noise. The model in not sensitive to this initial condition, i.e.
the final equilibrium slope is only a function of the chosen set of parameters (mainly the input rate i R and the entrainment 10 parameter E). Different initial conditions determine how long the system needs to reach this equilibrium state.

Rough-Bed and Jamming in CAST
CAST operates in two modes, with and without dynamic jamming. The rough-bed model without jamming (CAST RBM ) simulates a generic rough-bed channel where processes of transport, entrainment and deposition are considered regardless of any additional granular effect (except for particle collisions and deposition after collisions with the channel banks). The 15 jamming model (CAST JM ) simulates explicitly the process of jamming (blocking) when the density of particles transported in the same cross-section exceeds a threshold. In this case particle interactions lead to deposition of all grains in that cross-section on the bed. This blocking process is considered permanent, i.e. the jammed particles are locked into channel width spanning structures which cannot be entrained anymore. In the same cross-section, entrainment and deposition of other particles is still possible, except for those grains that have been subjected to jamming. This parameterization aims to represent in a simplified 20 way the additional force chains that keep grains together around keystones, as shown in the jammed-state hypothesis of Church and Zimmermann (2007). Intuitively and similarly to other phenomena where jamming is common (e.g. in hoppers) we set up the jamming threshold equal to the channel width Y . In a one grain-size model like CAST this implies that jamming is happening when the transport layer is full of particles (one entire cross-section full of transported particles).
For every time step, the computation sequence is as follows: (1) Sediment input enters the system in the first cross-section. 25 (2) Particles in transport move one cell downstream (straight, left or right), if they collide with other particles or with one of the banks they deposit, otherwise they remain in transport. In the case of CAST JM the jamming condition is checked for every cross-section. When the number of particles traveling in the same cross-section exceeds the jamming threshold all grains are deposited on the bed and frozen.
(3) For every particle in motion the condition for deposition is checked: if it is satisfied the particle leaves the transport matrix and is deposited on the bed. (4) For every particle in the bed matrix (except for those 30 jammed in stage 2 and deposited in stage 3) the condition for entrainment is checked: if it is satisfied the particle leaves the bed matrix and joins the transport matrix. (5) The boundary condition at the channel outlet is applied.
To set-up the model, we have run a set of simulations using CAST RBM with (a) stochastic entrainment with constant E and S parameters, and (b) constant specific input rate i R . The domain is Y = 20d and X = 300d. The domain size has been chosen not to represent or scale any specific river channel, but rather to observe the features and test the hypotheses we are interested in, at a reasonable computational cost. Simulations with a much larger scale (up to 2 orders of magnitude) have also been 5 performed and we observed no significant scale effect on the final outcome (see Supplementary material). The simulations were run until a steady state was reached. This condition is achieved when for a given combination of i R and E the channel slope does not change, i.e. the point at which the particle count (stored sediment volume) in the channel reaches a steady state and the sediment output is on average equal to the input. We explored the effect of different input rates i R and the entrainment parameters E on the final bed structure. Moreover, since CAST is a stochastic model, we perform 20 realizations for every set 10 of parameters to quantify stochastic variability. The set of parameters used in these steady-state simulations is shown in Tab. 1.
We analyzed the results in terms of: mean relative exposure < R >: the spatially-averaged value of R of all the cells in the control volume domain evaluated with Eq. 3; and standard deviation of the relative exposure σ R ; particle hop distance HD: the step length of a single particle from the point it is entrained (or enters the channel) to the point it is deposited or exits the channel. The main outputs of CAST RBM in a steady-state simulation (E = 1.5; i R = 0.5) at equilibrium are shown in Fig. 3. First, the storage volume (Fig. 3a) exhibits a dynamical equilibrium: the volume oscillates, alternating phases of aggradation and degradation, because the stochastic parameterization of particle entrainment leads to a continuous exchange between the bed and the transport layer. Second, the times series of sediment transport (Fig. 3b) show that even with a constant input rate, 25 sediment transport fluctuates as observed in the field and in the lab (e.g. Recking et al., 2009;Saletti et al., 2015). Moreover, CAST RBM produces a realistic rough-bed morphology (Fig. 3c), with the mean R and standard deviation σ R being a function of the input rate i R and the entrainment parameter E. The input rate and the entrainment parameter also determine the final slope of the channel.
One of the advantages of a RC models like CAST is that it is possible to track the movement of every single particle in the 30 system and so to compute all particle step lengths (measured from entrainment to deposition). This is an important quantity which, since Einstein's probabilistic theory on bed load transport (Einstein, 1937(Einstein, , 1950, needs to be reproduced by any reliable particle-based transport model. Since the word 'step' in this paper refers to the bed structures created by grain arrangements, for the sake of clarity we will call hereafter, following Furbish et al. (2012), particle hop distances (HD) what in literature is usually called step length, i.e. the distance travelled by a particle from entrainment to deposition (Fig. 3d). For every simulation we computed values of HD for all the particles and we find they follow an exponential distribution. In Fig. 4 we show for 5 4 different combinations of i R and E the probability density functions of HD and the exponential fit. The good fit given by this distribution is in agreement with previous studies dealing with particle travel distances (e.g. Hill et al., 2010;Hassan et al., 2013;Schneider et al., 2014). Since in the model no HD distribution is specified a priori, the agreement shows that the phenomenological rules of particle entrainment, transport, and deposition of CAST are realistic. The values of HD obtained in our simulations are always much smaller than the channel length X (e.g. the maximum observed HD is less than X/2), 10 therefore the system scale is not influencing our results. The relation between HD and the model parameters is explored in the next section.

Role of Input Rate and Entrainment Probability
With the steady-state simulations we explored the effect of changing input rate and entrainment parameter on the model outcomes. These two parameters are important because they can be linked to the jammed-state diagram parameters of Church 15 and Zimmermann (2007). The input rate I R is related to the sediment concentration Q S Q , which quantifies the effect of sediment supply on step stability. The entrainment parameter E determines the entrainment probability and is directly related to the transport stage τ τcr which quantifies the effect of the hydraulic forces on step stability. Some of the simulations, characterized by low input rate and high entrainment probability (E = 1, i R < 0.4), yield what we call 'washed-out' case, i.e. the bed matrix remains empty. This represent a limiting case where hydraulic forces are too high 20 and sediment supply is too low to be able to sustain a fluvial channel. This constrained our parameter space to 27 simulations in which a channel was formed and maintained around an equilibrium point.
The stochastic simulation of 20 realizations of each of the 27 parameter sets showed that the mean storage volume V and the mean relative exposure < R > converged to the same values. Also the mean hop distances < HD > and the standard deviations of the relative exposure σ R did not change significantly. 25 The values of 4 key variables for the 27 simulations (parameter combinations), averaged over the 20 realizations, are shown in Fig. 5. The mean relative exposure (< R > in Fig. 5a) is obtained by a spatial average of all the values of R for a given time step, and then temporally averaged over the last 20000 time steps in the equilibrium phase. < R > is directly related to the slope of the channel and the storage volume. It increases with increasing input rate and increasing entrainment parameter: channels with large sediment supply and low entrainment probability are those with larger < R > and larger storage volumes. 30 The same trend can be inferred by looking at the mean standard deviation of R (σ R Fig. 5b), also obtained as a spatial average over the equilibrium phase for every time step.
Specific sediment flux q S (t) is on average equal to the input rate i R at equilibrium, but fluctuating around its mean value, as shown in Fig. 3b. The degree of memory of these fluctuations is captured by the Hurst exponent H qs , whose mean value for the steady state simulations is shown in Fig. 5c. The values of H q S obtained in all the realizations of all the simulations are consistent with those obtained from flume experiments by Saletti et al. (2015), being in the interval [0. 5 : 1]. This identifies a long-memory regime which is stronger in the model (H → 1) when the entrainment probability is high (low E) and the input rate is low. H q S shows large variability in different realizations, although its mean value shows a clear trend with both i R and E (Fig. 5c). 5 The mean particle hop distances (< HD > in Fig. 5d) display two contrasting trends. For values of E ≤ 1.5 (high entrainment probability) < HD > is decreasing consistently for increasing input rates i R as a consequence of large particle activity (collisions between particles become very frequent). For larger values of E (low entrainment probability) the maximum of < HD > is for average values of i R (around 0.4). For i R < 0.4 there is a stronger interaction with the bed (which has larger < R > and σ R for larger E, as shown in Fig. 5a,b) and so more likelihood of particle deposition, whereas, for large i R , collisions 10 again dominate hop distances because of large particle activity. In both cases this leads to a reduction of < HD >.

CAST JM : Particle Jamming
The analysis above showed that CAST RBM is able to reproduce fluvial channels with a spatially variable rough bed and to capture basic and important physical phenomena, such as the variability in sediment flux and the exponential distribution of particle hop distances. However, to simulate the formation and test the stability of steps, we need to take into account the effect 15 of particle jamming. This is a well-studied phenomenon in granular physics that has been advocated to be essential in the step stability process and considered through the jamming ratio (the ratio between the channel width and d 84 of the surface) in the diagram proposed by Church and Zimmermann (2007). In our reduced-complexity model CAST JM we account for the jamming effect by blocking particles when local sediment concentration exceeds the jamming threshold, and depositing them in permanent structures on the bed. 20 Jamming simulations were run with the same parameter sets of the steady-state rough-bed model case. Three different situations occur: 1. When particle activity is too low (low sediment transport) the jamming threshold is rarely (often never) exceeded.
2. When particle activity is too high (high sediment transport) jamming is occurring too often in time and space, and the storage volume of the system keeps increasing because of the large amount of particles depositing upstream of the step 25 structures. As a result an equilibrium channel is never reached.
3. When particle activity is in-between the two previous situations, jamming is occurring at a rate which allows the formation of steps and maintains an approximately equilibrium channel.
The first situation represents the rough-bed case discussed previously. The second one represents a case which is very unlikely to happen in river systems where fluvial sediment transport is rarely going to exceed the jamming threshold and 30 certainly not for very long periods of time (e.g. only during large flood events). For the purpose of this study we focus on the last situation where jamming is effectively creating steps. When particles are jammed and instantly deposited, they trigger a deposition process which is propagating upstream, since the values of relative exposure R will be immediately reduced.
This represents what happens in natural step-pool systems, where steps are created, among other factors, by deposition and clustering of sediment around large boulders called keystones, and deposition between steps continuously changes the channel (e.g. Molnar et al., 2010).
We show the effect of adding jamming to the model by comparing simulations having the same parameter sets and same 5 initial conditions (i R = 0.5 and E = 1.25) in CAST RBM and CAST JM runs. The cumulative number of jammed crosssections shows that jamming is a rather intermittent phenomenon with many long periods of no jamming (Fig. 6a). At the end of this simulation 29 cross-sections were jammed (around 10% of the total). The longitudinal profiles of bed elevation (Fig.   6b) show how CAST JM is able to create step structures and this increases the total slope of the channel and its storage, even if the slope between steps is the same as in the case without jamming (CAST RBM ). The boxplots of the instantaneous (i.e. 10 calculated for every time step) values of R, both the mean < R > (Fig. 6c) and the standard deviation σ R (Fig. 6d), show that the model with jamming yields a rougher and more variable bed. Instantaneous values of specific sediment flux q s (Fig. 6e) show that jamming slightly increases the variability in q S and prevents the formation of an equilibrium slope (which would imply q S i R = 0.5) because the system is still aggrading and increasing its storage. Finally the values of the Hurst exponent of sediment flux H q S (Fig. 6f) for the 20 realizations clearly plot separately in the case with and without jamming, with the 15 latter having much greater values. This longer-term memory is likely due to a combination of sediment pulses created by step collapses and the weak but present trend towards aggradation in the CAST JM simulations.

The Effect of Stochasticity
The processes of particle entrainment and sediment supply can be parameterized as deterministic or stochastic. In the simulations presented in the previous sections we used a stochastic parameterization for particle entrainment and a constant sediment 20 supply. To explore the effect of stochasticity on the model results, in the next two sections we quantify the effect of stochasticity in entrainment and sediment supply explicitly.

Stochasticity in the Entrainment
The entrainment probability in CAST can be parameterized as a deterministic or stochastic process (Section 2.1.3). A stochastic parameterization allows a degree of variability in the entrainment threshold and can be controlled by two parameters (E and 25 S), while the deterministic parameterization has a unique entrainment threshold E (Fig 2).
The comparison for a simulation with i R = 0.5, E = 1 and S = E/5 is shown in Figure 7. When the entrainment process is treated as stochastic, the variability of sediment flux is much larger both in case with and without jamming (Fig. 7a). This is due to the fact that when the channel has reached equilibrium in the deterministic case the interaction between the bed and the transport is very low. All particles below the threshold stay on the bed and those above the threshold are entrained. The reduced 30 particle activity can be inferred also by looking at the distribution of particle hop distances (Fig. 7b). In the deterministic case the distributions are shifted towards larger values because particles interact much less with the bed and travel further downstream. The effect of modeling the entrainment as a deterministic process on the bed morphology itself is that the final configuration of the channel in the threshold case is much steeper (cumulative distribution functions of R plot towards larger values in the T-case) since no entrainment is possible below the threshold: the channel can bear steeper slopes and store more sediment (Fig. 7c). However, this does not translate into a rougher surface; σ R shows that channels where the entrainment is modeled with a threshold function have very low variability around < R >: they tend to look more like steep and uniform 5 ramps than like realistic fluvial channels (Fig. 7d).
This analysis support our choice of modeling the entrainment as a stochastic process. This is not only more physically reasonable because the process of particle displacement is random per se, but it is not possible to obtain a realistic rough-bed morphology in a reduced-complexity model like CAST without a stochastic parameterization of particle entrainment.

10
The effect of stochasticity in the input rate i R is shown on simulations with E = 1.5 and constant i R = 0.5, or random input rate uniformly distributed around the mean value < i R >= 0.5.
The effect of stochasticity in the input rate is much smaller than that of entrainment. The distributions of sediment flux ( Fig. 8a) almost overlap both in the case with and without jamming. Like in the entrainment case, jamming causes a constant increase in sediment storage (and so the median transport rate is slightly below the equilibrium value of 0.5). An overlap of 15 the 4 simulations is also observed when looking at the distributions of particle hop distances and final R (Fig. 8b and c). The simulations with and without jamming plot separately only in the case of σ R (Fig. 8d). The larger values of σ R in the case of variable i R and jamming are due to the fact the variability in the input facilitates the jamming process and increases the relative exposure.
These results highlight that the variability and the fluctuations observed in the sediment output variables of the model do not 20 depend on the variability of the sediment input, but are instead function of the internal dynamics of the system, given by the local grain-grain and grain-bed interactions. In other words CAST acts as a shredding filter of the input forcing (Jerolmack and Paola, 2010;Van De Wiel and Coulthard, 2010).

Unsteady Simulations
Although jamming is effective in generating a step-like morphology under certain steady-state sediment input and entrainment 25 conditions, we recognize that step formation is an intermittent process in which flow variability in time is important. Typically step-pool sequences are partially or totally destroyed during large flood events and then reworked and stabilized during the following low flows periods (e.g. Lenzi, 2001;Turowski et al., 2009;Molnar et al., 2010). We show the effects of changing flow conditions by simulating 4 consecutive floods of equal magnitude (Fig. 9a). In CAST the hydraulic conditions are represented by the entrainment parameter E. Therefore, to simulate a change in the flow, we modify the value of E to represent 30 two extreme cases (Fig. 9b): low flow with E = 2 (low entrainment probability) and high flow with E = 1 (high entrainment probability). Moreover, we explore two different situations: (1) we keep the input rate i R constant, incorporating all the effects of the unsteadiness in the entrainment parameter E (Case I in Fig. 9c), and (2) we change also the input rate i R in response to changes in the flow conditions (Case II in Fig. 9d). To facilitate the comparison, the total sediment input over the entire simulation is the same in Case I and II. To mimic the rising and falling limb of an hydrograph, both E and i R were increased and decreased gradually. The relevant parameters of these unsteady simulations are summarized in Tab. 2. Runs were performed both with and without jamming (i.e. with CAST RBM and CAST JM ), to check if and when steps are formed and how many 5 of them remain stable.

Sediment Storage and Sediment Transport
The temporal pattern of storage volume and sediment flux in the unsteady simulations is shown in Fig. 10. The volume for the rough-bed case (in blue) clearly displays phases of degradation during high flow and aggradation during low flow. Without jamming the channel tends to erode during high flows when the entrainment probability is high and to gain sediment again 10 during low flow when the entrainment probability decreases. This turnover is more evident when the input rate is constant (Fig.   10a). With the effect of jamming the picture changes. During high flows the mobile grains are trapped in the channel in steps, while during low flows the channel increases its storage because of grain deposition between steps and channel infilling. With a variable input (Case II), jamming creates more steps and increases the storage volume which then remains constant during the following low flow phases because of the reduced input rate and low entrainment probability (Fig. 10b). 15 The specific sediment flux when the input is constant (Fig. 10c) shows a large variability for the rough-bed case responding to changes in E during low and high flow conditions. In the jamming case, the response to the change in flow conditions is also present but the jamming process modulates the sediment flux towards the equilibrium conditions rapidly. When the input varies with flow conditions (Fig. 10d), the rough-bed model yields the same pattern as the case with constant input with the difference that here the equilibrium condition is changing during low and high flow (0.3 and 0.6 respectively). In the jamming model 20 instead, the sediment flux is almost instantly in equilibrium with the input rate during low flows, while during high flows, the large input rate, together with the high entrainment probability, causes so many jamming events that inhibit the system to reach an equilibrium state and the channel keeps increasing its storage. We show the statistical distributions of specific sediment flux for the 4 cases in Fig. 11. It can be seen that the distributions are centered around the equilibrium point of 0.4 (especially the jamming case in red for Case I), with the rough-bed model having a more spread function due to the more intense phases of 25 aggradation and degradation. The distributions of Case II (dashed lines) are clearly bimodal because of the two equilibrium sediment input rates (0.3 and 0.6).

Step Formation and Stability
The unsteady flow also has impacts on bed roughness in CAST . The time series of the standard deviation σ R , which represent the degree of roughness of the bed, is shown in Fig. 12 for the unsteady simulation with constant and variable input rate. In 30 both cases jamming produces a rougher surface during high flow which is a indication that step structures, causing a larger departure from the mean R, are being formed. When the input is constant (Fig. 12a), σ R goes back to value of low flow for all the 4 floods, because steps that were formed are being buried by sediment. When instead the input rate is reduced during low flow to simulate sediment-starved conditions (Fig. 12b), σ R decreases but not to its pre-flood value because many of the steps created during high flow can survive and do not get buried in-between floods.
The same can be inferred from the longitudinal profiles of bed elevation of the simulations with jamming (Fig. 13). At the end of every high-flow period, the longitudinal profile shows a stepped morphology due to jamming. In the following low-flow periods the steps were buried in Case I (having input rate i R = 0.4), while in Case II (having a lower input rate during low 5 flow: i R = 0.3) some of them survived because of the sediment-starved conditions.
To quantify this effect directly on step formation, we introduce step density d S , defined as the ratio between the number of cross-sections with steps and the total number of cross-section of the channel. The variable d S can vary between 0 when no steps are present in the channel, and 1 when all the channel morphology is made by steps. The definition of a step is not straightforward, even in the field and in the lab where many different identification algorithms have been proposed (e.g. Milzow 10 et al., 2006;Zimmermann et al., 2008). Since our goal here is not to identify and count the number of steps or to test which step identification algorithm works best, we simply define a step in terms of local departure from the equilibrium channel slope, similarly to the method of Milzow et al. (2006). The steady-state simulations give us the value of the final slope at equilibrium for a given set of parameters (E and i R ). We define that a cross-section in CAST has a step if its local slope is greater than the equilibrium slope by a factor β. The time-series of step density evaluated in this way is shown for different values of β 15 in Fig. 14. The temporal pattern of step density variations is largely independent of β. The time evolution of step density as a function of flow and sediment supply conditions allows us to draw two conclusions. First, there is a clear difference between simulations with and without jamming in that jamming is responsible for step formation, and without it there are practically no steps formed in the channel (blue and red lines in Fig. 14). Second, after steps are generated during high-flow periods due to jamming, they only survive during low flow if the sediment supply decreases (yellow lines in Fig. 14), in sediment-starved 20 conditions. This illustrates the temporal dynamics of step counts as observed in the field (e.g. Molnar et al., 2010) and in flume experiments (e.g. Curran and Wilcock, 2005).

Bed load: a Stochastic, Granular and Shredding Phenomenon
The CAST model without jamming, CAST RBM , simulates bed load transport over a rough-bed at the grain scale, considering 25 particle entrainment as a stochastic process driven by a local exposure. Our model describes bed load from a grain prospective because local granular effects in particle mobility and transport are key for developing a bed morphology, especially in steep and well-structured streams. Kirchner et al. (1990) pointed out the role of granular interactions between gravel particles on a river bed and showed how the erodibility of a grain is controlled by its protrusion and friction angle. Given the associated high variability, they also suggested that, instead of using one single value for the shear stress, a probabilistic approach should 30 be applied. The role played by particle interlocking and partial burial in increasing measured friction angles in steep channels has been shown recently also by Prancevic and Lamb (2015). Many laboratory studies have increased our knowledge of bed load transport exactly by looking at the granular scale (e.g. Lajeneusse et al., 2010;Houssais et al., 2015), suggesting that we might be more successful in describing this phenomenon when borrowing concepts from the granular physics community (e.g. Church and Zimmermann, 2007;Frey and Church, 2011).
CAST assumes a stochastic description of sediment transport, following and corroborating recent research Roseberry et al., 2012;Heyman et al., 2014;Ancey and Heymann, 2014). Our model produces fluctuations in transport rates by the interaction with the bed through entrainment and deposition of individual particles (Ancey and Heymann, 2014) and 5 these fluctuations are observed in our simulations even with a constant input forcing. What in our model is defined as specific sediment flux q S is equivalent to the particle activity as defined in Furbish et al. (2012). In a companion paper Roseberry et al. (2012) found that changes in transport rates are dominated by changes in the number of particles in motion rather than velocity: this justifies our choice of assuming constant particle velocity but varying entrainment threshold in our simulations.
The stochastic parameterization of CAST does not assume a priori any probability distribution for particle hop distances, 10 and yet they turn out to be well fitted by an exponential distribution, in agreement with previous theoretical and field studies (e.g. Hill et al., 2010;Hassan et al., 2013;Schneider et al., 2014). The fact that, despite its simplicity, the model can reproduce in a robust way this important feature, proves at least partially that the local grain-grain and grain-bed interaction rules in CAST are appropriate and that the phenomenological descriptions of the simulated phenomena are going in the right direction.
Finally, CAST reproduces also the shredding effect sometimes visible in sediment transport (Jerolmack and Paola, 2010). 15 The measured variability of sediment flux and its fluctuations are dictated by the internal dynamics of the system and the degree of fluctuations in the input forcing does not always affect the sediment flux in a clear way (see Fig. 8). Our results then show the RCM potential of modeling of bed load transport as a stochastic phenomenon at the grain scale. Interactions between individual particles can give rise to, or at least strongly impact, the variability observed in natural fluvial systems.
Reduced-complexity models like CAST can serve to model these interactions and their effects and can be used to gain new 20 insights into the complex dynamics of sediment transport and to test new research hypotheses.

Step Formation and Stability: a Granular Problem
We showed with CAST JM that dynamic jamming of particles in motion is effective in forming steps (see Fig. 13). In fact only by including the jamming process we did generate step-pool like morphologies in our numerical experiments. Moreover, once steps are formed, they remain stable if the flow conditions change (i.e. the entrainment probability decreases) and the supply of 25 sediment is low enough to avoid that these steps are buried by particles (as shown in Fig. 14). These results are consistent with the main ideas of the jammed state hypothesis of Church and Zimmermann (2007) who theorized and showed experimentally (Zimmermann et al., 2010) that step stability needs (a) jamming, expressed as a low width to diameter ratio so as to enhance granular forces, (b) low flow stage, in order to avoid the mobilization of keystones, and (c) sediment-starved conditions, because a too high sediment concentration would bury the steps. Despite its simplifications, especially uniform sediment and no explicit 30 flow parametrization, CAST JM can reproduce these observations and support the jammed state hypothesis for step stability.
We did not observe any specific wavelength of step occurrence, as usually predicted by hydraulic-based theories on step formation (e.g. Whittaker and Jaeggi, 1982). Given the stochastic nature of CAST , steps due to jamming are not formed with a regular spacing. Moreover, as it has been shown by more recent experimental (Curran, 2007;Zimmermann et al., 2010) and field studies (Zimmermann and Church, 2001;Molnar et al., 2010), step occurrence is mainly driven by the random location of boulders (i.e. keystones) around which sediment deposits and clusters.

Outlook
Our modeling approach has by definition some simplifications and limitations which we think can be improved in future research. First, the uniform size of the sediment prevents us to model specifically any grain-size effect that might indeed be very 5 important in steep-channel dynamics. We partially incorporated these effects in the stochastic parameterization of entrainment: the fact for the same value of relative exposure R and entrainment parameter E some particles are displaced and some are not accounts also for differences in their dimension and weight. Also the jamming process may be grain-size dependent, as well as the particle velocity. Second, the parameterization of changing flow conditions is done indirectly, summarized entirely in the entrainment parameter E. This is done mainly because we are not aiming to model discharge, flow, shear stress on the 10 bed, but rather transfer their effects onto the probability of entraining grains. However, future improvements of CAST could include a more direct relation between hydraulic stresses on the bed and the E and S parameters in our model. Third, the granular interactions (i.e. collisions) among particles always lead to deposition, which might not always be realistic, at least in fluvial systems with particles having different sizes and shapes. The same can be said about interactions with the banks of the channel. However, in a uniform-size case this assumption does not seem to be too strong. Furthermore, we did not account for 15 the transfer of momentum that could happen when a particle is deposited and so enhancing the probability of entrainment of the surrounding grains (i.e. "collective entrainment" as in Ancey and Heymann (2014)). With a model like CAST the relative importance of this phenomenon in the entrainment process could be evaluated. Finally, the choice of representing the jamming of grains on the bed as a permanent process is a limitation. In a model having only a single grain size, this choice has been made to account for the additional granular forces that are making step structures more stable around a keystone. In future 20 research, especially in a model which accounts for different grain fractions, the role of step formation and stability may be transferred to the coarsest grains to which jamming threshold will apply.
In conclusion, in our opinion the strongest limitation of the current model is the absence of sediment-sorting and other grainsize effects. All these phenomena will be incorporated in the next version of the model which will have different grain-size fractions. 25

Conclusions
We presented a new particle-based reduced-complexity model CAST (Cellular Automaton Sediment Transport) that simulates bed load transport and changes in channel morphology including the processes of jamming and step formation. The model simulates grain-grain and grain-bed interactions with uniform-size particles and can have stochastic or deterministic parameterizations for sediment input rate and particle entrainment. With only few parameters, it it possible to simulate channels with 30 different sediment supply and flow conditions. At steady state, CAST can reproduce a realistic bed morphology and typical fluctuations in transport rates, whose memory features are consistent with previous experimental data. Moreover, particle hop distances are well-fitted by exponential distributions, in agreement with field observations. One of the main results is the role played by stochasticity both in the entrainment and in the input rate. A stochastic input rate does not change the final outcome of the model compared to a constant input having the same mean. However, if the entrainment is modeled deterministically, the resulting channel does not have the typical variable bed roughness encountered in real fluvial systems. The dynamical effect of particle jamming was added to test under which conditions steps are formed and remain stable in 5 steep channels. The effect of jamming has been tested in unsteady simulations where the entrainment probability and the input rate have been changed to simulate a sequence of high-flow and low-flow periods. CAST generates step structures during highflow periods that survive during low flows in simulations with sediment-starved conditions, in agreement with the jammed-state hypothesis. Our results support the jammed-state hypothesis as a framework to explain step formation and stability and, more in general, they show the potential of reduced complexity-models at a grain scale with stochastic parameterizations. We are of 10 the opinion that models such as CAST can give new insights into the dynamics of complex phenomena like sediment transport and step formation and be useful to test research hypotheses in fluvial geomorphology.  Particles can be either in the bed matrix or in the transport matrix. They can enter the transport matrix as sediment input from the upper boundary or by entrainment from the bed, while they can leave the transport matrix as sediment output or by deposition on the bed. Sediment is input at the upstream boundary and simulated as sediment yield leaving the downstream boundary.  Bed elevation in the final configuration: the model produces a rough bed with particles having different exposures R. (d) Probability density function of particle hop distances with an exponential distribution estimated over more than 2 million simulated particle paths.    Low flow duration T low 70000 70000