A Well‐Posed Alternative to the Hirano Active Layer Model for Rivers With Mixed‐Size Sediment

The active layer model (Hirano, 1971) is frequently used for modeling mixed‐size sediment river morphodynamic processes. It assumes that all the dynamics of the bed surface are captured by a homogeneous top layer that interacts with the flow. Although successful in reproducing a wide range of phenomena, it has two problems: (1) It may become mathematically ill‐posed, which causes the model to lose its predictive capabilities, and (2) it does not capture dispersion of tracer sediment. We extend the active layer model by accounting for conservation of the sediment in transport and obtain a new model that overcomes the two problems. We analytically assess the model properties and discover an instability mechanism associated with the formation of waves under conditions in which the active layer model is ill‐posed. Numerical simulations using the new model show that it is capable of reproducing two laboratory experiments conducted under conditions in which the active layer model is ill‐posed. The new model captures the formation of waves and mixing due to an increase in active layer thickness. Simulations of tracer dispersion show that the model reproduces reasonably well a laboratory experiment under conditions in which the effect of temporary burial of sediment due to bedforms is negligible. Simulations of a field experiment illustrate that the model does not capture the effect of temporary burial of sediment by bedforms.


Introduction
A common approach in modeling river morphodynamic changes in space and time consists of solving a set of differential equations that account for the flow and bed changes. As the hydrostatic pressure assumption is typically valid in fluvial problems, the flow is often modeled in two dimensions using the Shallow Water Equations (e.g., Tan, 1992;Vreugdenhil, 1994), which reduce to the Saint-Venant (1871) equations under one-dimensional conditions. Bed elevation changes are accounted for using the Exner (1920) equation and changes in the bed surface texture using the active layer model (Hirano, 1971). Changes in bed elevation and surface texture depend on the gradient of the sediment transport rate, which is a function of the local bed shear stress. This modeling framework has been applied for decades (e.g., Bennett & Nordin, 1977;Ferguson et al., 2015); it is implemented in major software packages (e.g., Mosselman & Sloff, 2007;Sloff & Mosselman, 2012;Vetsch et al., 2006;Villaret et al., 2013) and has proven to be a powerful tool able to model processes from the lab to the field scale (e.g., Cui et al., 2003;Orrú et al., 2016;Williams et al., 2016). However, the above modeling approach has two drawbacks. First, the solution may be invalid as the system of equations may become mathematically ill-posed (Chavarrías et al., 2018;Ribberink, 1987;Stecca et al., 2014). In ill-posed problems, the growth rate of perturbations tends to be infinite as the wavelength decreases (Hadamard, 1923;Joseph & Saut, 1990;Kabanikhin, 2008). This is physically unrealistic, as natural phenomena are subject to short-wavelength perturbations (e.g., noise) that do not grow unbounded. Numerical solutions of ill-posed problems continue to change as the grid is refined (i.e., do not converge with the grid) (Joseph & Saut, 1990). Chavarrías et al. (2018) show an example of the consequences of ill-posedness in river morphodynamic simulations. As the model becomes ill-posed, perturbations due to numerical truncation error grow up to the scale of the flow depth and alter the stratigraphy in an unrealistic manner. Ill-posedness is a symptom of a model not capturing key physical elements (Fowler, 1997;Joseph & Saut, 1990). In our context, the active layer model is ill-posed when key mixing processes that are not included in the model become relevant and the model is incapable of reproducing the actual mixing occurring in nature (Chavarrías, Stecca, et al., 2019). This occurs mainly under degradational conditions and if the active layer is coarser than the substrate (i.e., when modeling degradation of an armored river) (Chavarrías et al., 2018;Ribberink, 1987;Stecca et al., 2014). The active layer model may also be ill-posed under degradational conditions when the active layer is finer than the substrate and under aggradational conditions (Chavarrías et al., 2018). However, the likelihood that the model is ill-posed when modeling degradation of an armored river is larger than for all other conditions. In order to solve for this first limitation and recover well-posedness of the active layer model, Chavarrías, Stecca, et al. (2019) proposed a regularization strategy. The strategy is based on identifying the locations at which the active layer model is ill-posed and locally modifying the celerity at which mixing processes occur. In essence, the strategy alters the celerity of the mixing processes predicted by the active layer model, such that the regularized model is well-posed. Chavarrías, Stecca, et al. (2019) conducted a set of laboratory experiments under conditions in which the active layer model is ill-posed to test the regularization strategy. In particular, the experiments were characterized by degradation of a coarse bed surface into a substrate consisting of fine sediment. The situation appeared to be unstable: Bedforms composed of coarse sediment grew due to entrainment of fine sediment in the troughs and decayed when coarse sediment coming from upstream covered the bed. The oscillations due to the instability mechanism were superimposed on the overall degradational trend. The regularized model captured the overall degradational trend but did not reproduce the oscillations due to the instability mechanism (Chavarrías, Stecca, et al., 2019). Another drawback of the regularization strategy is the fact that it can only be applied if the active layer thickness remains constant with time. This is limiting, as the active layer thickness is usually related to time-varying properties such as a characteristic grain size of the bed surface (Petts et al., 1989;Rahuel et al., 1989) or bedform dimensions (Armanini & di Silvio, 1988;Deigaard & Fredsøe, 1978).
The second limitation of the active layer model is related to sediment dispersion. Under dynamic equilibrium conditions, the active layer model predicts tracer sediment to be advected downstream at a mean speed inversely proportional to the active layer thickness without dispersing (Iwasaki et al., 2017). However, traced sediment particles are observed to disperse as they move downstream both in the field (e.g., Bradley & Tucker, 2012;Bradley, 2017;Drake et al., 1988;Hassan et al., 1991;Nikora et al., 2002;Rathbun et al., 1971;Sayre & Hubbell, 1965) and in the laboratory (e.g., Hill et al., 2010;Martin et al., 2012;Roseberry et al., 2012). The active layer model is presented in two forms (Parker et al., 2000): (1) the flux form, in which the sediment transport rate is computed and changes in bed elevation and bed surface texture depend on the gradient of the sediment transport rate, and (2) in entrainment-deposition form, in which the rates at which sediment is entrained and deposited are computed and changes occur due their difference. The inability of the active layer model to capture dispersion does not depend on the form of the equation, as both forms are equivalent (Parker et al., 2000). For capturing sediment dispersion, it is necessary to account for a mechanism that causes sediment transport not to be dependent on the local conditions only.
Dispersion of tracer sediment is captured when sediment transport is treated in a stochastic manner, which was first shown by Einstein (1936). He described sediment transport as particles traveling in a series of jumps of varying length and frequency. This allowed for deriving the probability that a tracer particle is at a certain location as a function of time, although the model was not linked to flow and bed elevation changes. Einstein's (1936) theory was applied by Sayre and Hubbell (1965) and Habersack (2001), among others, to explain tracer dispersion in field cases in the United States and New Zealand, respectively. Recent refinements of Einstein's (1936) theory include, for instance, the work by Fan et al. (2016Fan et al. ( , 2017, who considered the effect of the travel time in simulating tracer motion. However, just as in Einstein (1936), their approach to sediment transport is not combined with flow and morphodynamic change.
A different approach for modeling tracer dispersion consists of accounting for the number of static and moving particles in the bed. In this way, Lajeunesse et al. (2013Lajeunesse et al. ( , 2017Lajeunesse et al. ( , 2018) model dispersion of a plume of unisize tracer sediment under normal flow conditions. Ancey et al. (2006Ancey et al. ( , 2008 considered the statistics of unisize sediment transport and showed that this can be well captured using a discrete statistic theory (i.e., birth-death discrete Markov processes). A continuum version of the theory was derived by Ancey (2010) and Ancey and Heyman (2014), who accounted for the variability in particle velocity. This led to an advection-diffusion equation modeling the ensemble average of the volume of sediment in transport per unit of bed area. The diffusive behavior of the equation explains tracer dispersion. Bohorquez and Ancey (2015) coupled the model developed by Ancey and Heyman (2014) to the Saint-Venant (1871) flow equations to model the formation of antidunes. Interestingly, an advection-diffusion equation for modeling the volume Our objective is to develop a model that accounts for mixed-size sediment river morphodynamic changes and overcomes the limitations mentioned above. In particular, we aim at a well-posed model that captures tracer dispersion and the instability mechanism due to the entrainment of fine substrate sediment. To this end, we combine the active layer model with the model by Bohorquez and Ancey (2015), after we extend the latter to mixed-size sediment conditions.
In section 2, we present the mathematical model and the numerical technique to solve the set of equations. In section 3, we linearize the system of equations and study the growth rate of perturbations to explain the instability associated with the entrainment of fine substrate sediment with respect to the one present at the bed surface. We follow this approach, as it has explained the formation of other patterns such as ripples (Colombini & Stocchino, 2011;Sumer & Bakioglu, 1984), dunes (Colombini & Stocchino, 2008;Kennedy, 1963), bars (Callander, 1969;Schielen et al., 1993), and meanders (Ikeda et al., 1981;Seminara, 2006). We also use the linear model to corroborate the well-posedness of the system of equations and to study the manner in which tracer sediment disperses. In section 4, we run numerical simulations to test the performance of the model against measured data and to compare the results to those of the active layer model. We model two cases of tracer propagation and two cases under conditions in which the active layer model is ill-posed. In one of the tracer cases, the model performs reasonably well, while the second shows the model deficiencies. The first ill-posed case shows the capabilities of the model in reproducing the instability mechanism, and the second is a case that only our model can reproduce. In section 5, we discuss the results.

Conservation Equations
The SedIment Layers with source-sinK Exchange (SILKE) model for mixed-size sediment river morphodynamics accounts for one-dimensional flow over a bed composed of an arbitrary number of noncohesive sediment fractions. The sediment fractions are characterized by a grain size d k (m), where the subscript k identifies each fraction in increasing size (i.e., d 1 < d 2 < ... < d N ) and N is the total number of size fractions.
We assume that the concentration of solid to liquid discharge is below 6×10 −3 , so that the effect of the concentration of sediment on the flow is negligible (Garegnani et al., 2011(Garegnani et al., , 2013. Under these conditions, unsteady depth-averaged flow is modeled by the Saint-Venant (1871) equations. Here, for simplicity, we assume steady flow and apply the backwater equation (e.g., Chow, 1959): where x (m) is the streamwise coordinate, h (m) the mean flow depth, (m) the bed elevation, Fr = u∕ √ gh (−) is the Froude number, u (m/s) denotes the mean flow velocity, g (m/s 2 ) is the acceleration due to gravity, and S f (−) denotes the friction slope. When using the backwater equation, we implicitly assume that the relevant processes occur on a length scale longer than the flow depth (e.g., Battjes & Labeur, 2017), as for waves shorter than the flow depth hydrostatic flow cannot be assumed. We also assume that changes in the flow field occur on a faster time scale than changes in the bed elevation and the bed surface texture (i.e., quasi-steady flow; Cao & Carling, 2002;Colombini & Stocchino, 2005;De Vries, 1965).
Similarly to Hirano (1971), we assume that sediment can be in three states: (1) in motion, predominantly in the streamwise direction; (2) at the bed surface, where it can be entrained into motion and deposited; or (3) in the substrate, where sediment is immobile and cannot be entrained. Subsequently, we will provide equations describing the conservation of mass of the sediment in each of the three states. Hirano (1971) did not consider conservation of the sediment in transport in deriving the active layer model. This was done by Armanini and di Silvio (1988), who extended the active layer model with two more layers to account for suspended load and bed load. In applying their model, they simplified it by assuming that the sediment transport rate is at capacity and that the temporal change of the volume of sediment in the bed load layer is negligible compared to the divergence of the sediment transport rate and the flux of sediment to the bed load layer. These two assumptions essentially remove the dynamics of the bed load layer, which is the key focus of our analysis. Bohorquez and Ancey (2015) describe the motion of unisize sediment based on a stochastic interpretation of sediment transport. We extend their framework to mixed-size sediment conditions by considering all variables to be grain size dependent. Conservation of the mass of sediment in transport is described by where t (s) denotes the time coordinate, E k (m/s) and D k (m/s) are entrainment and deposition rates of each size fraction k, and v pk (m/s) is the ensemble average of the instantaneous velocity of the moving particles of size fraction k. Parameter k (m 2 /s) is the diffusivity of size fraction k due to the variability in particle velocity from an Eulerian perspective (Furbish, Ball, et al., 2012;Furbish et al., 2017;Roseberry et al., 2012). This is essentially different from the Lagrangian diffusivity obtained from tracer dispersion experiments. Parameter Γ k (m) is the ensemble average of the volume of sediment in transport of size fraction k per unit of bed area. Following Furbish, Haff, et al. (2012), we will term Γ k the particle activity.
The particle activity can be interpreted as a fictitious layer with thickness Γ k . However, here no thickness of the bed load layer is specified as opposed to, for instance, Van Rijn (1984), Luu et al. (2004), Wu (2004), Colombini (2004), and Colombini and Stocchino (2005). Contrary to Armanini and di Silvio (1988), we do not make a distinction between suspended and bed load sediment. In our extension to mixed-size sediment conditions, we implicitly neglect the covariance terms that appear in equation (2) due to the correlation between particle size, velocity, and diffusivity (Furbish, Haff, et al., 2012). Worded differently, we consider sediment of different sizes traveling at different speeds, yet each size fraction is treated independently from the other fractions. The possible implications of these assumptions are assessed in section 5.2.
The sediment transport rate of size fraction k (q bk (m 2 /s)) does not depend on the flow properties only and is equal to This expression is a mixed-size sediment version of the expression derived by Furbish, Haff, et al. (2012), where we have assumed that all variables are grain size dependent.
Following Einstein (1950), we assume that the changes in bed elevation are due to the difference between entrainment and deposition rates. Mass conservation of the bed surface sediment (considering all sediment size fractions) yields the entrainment form of the Exner (1920) equation (Borah et al., 1982;Furbish, Haff, et al., 2012;Nakagawa & Tsujimoto, 1980a;Parker et al., 2000): where p (−) denotes the bed porosity, and E = ∑ N k=1 E k (m/s) and D = ∑ N k=1 D k (m/s) are the total entrainment and deposition rates. For simplicity, mechanisms such as subsidence and uplift, compaction and dilation of sediment are neglected in the above equation (Paola & Voller, 2005).
Following Hirano (1971), we assume that only the top part of the bed, characterized by a certain thickness (i.e., the active layer thickness), interacts with the flow (see recent discussions by Haschenburger, 2017, andAshmore et al., 2018). Sediment in the active layer is homogeneously mixed such that the layer has no vertical stratigraphy. Sediment is entrained from and deposited in the active layer, and a vertical flux of sediment between the active layer and the substrate occurs if the elevation of the interface between the active layer and the substrate changes.
The active layer thickness is typically defined as the vertical range covering a significant amount (e.g., 95%; Blom, 2008;Ribberink, 1987) of the probability density function of bed elevation. In particular, the probability density function of bed elevation is defined at such a time (space) scale that the fluctuating component (around the mean) is related to those physical processes relevant for sediment mixing that are unresolved by the model. As our model does not explicitly describe the formation and propagation of bedforms, mean bed level is obtained averaging over the passage of several bedforms, and the active layer thickness is related to bedform height.
Using the active layer model, one implicitly assumes that the probability of entrainment is constant within the active layer and equal to 0 in the substrate (Parker et al., 2000). This is an approximation; as in reality, Figure 1. Sketch representing the main variables of the active layer model (Hirano, 1971) and the SILKE model. The figure represents a longitudinal section. Both approaches model the bed elevation ( ) and grain size distribution at the bed surface (F ak ) averaged throughout the passage of several bedforms. In both cases, the substrate may have an arbitrary grain size distribution varying in the vertical direction (f sk (z)). The arrows indicate the varying amount of volume fraction content. Both models depend on the volume fraction content of sediment at the interface between the active layer and the substrate ( I k ). While in the active layer model, the sediment transport rate (q bk ) is computed using a closure relation, and the mass of sediment in transport is not accounted for, the SILKE model conserves the sediment in transport (Γ k , represented by small discs in the water column). The figure is adapted from Figure 2 in Chavarrías et al. (2018). there is no discrete distinction between the active and inactive parts of the bed. In reality, the probability of a particle of being entrained depends on its probability of being exposed to the flow (Blom & Parker, 2004;Galvin, 1965;Hassan & Church, 1994;Parker et al., 2000;Ribberink, 1987). Bed elevation fluctuations due to ripples, dunes, bars, or any other sort of fluctuations cause the probability of exposure to vary with vertical position. Particles deposited at lower elevations have a smaller probability of being reentrained. For this reason, using a single discrete active layer, it is not possible to capture the effect of temporary burial of sediment due to bedforms reworking sediment at an elevation significantly lower than the lower limit of the active layer (Blom, 2008;Ribberink, 1987).
The entrainment formulation of the active layer model reads (Parker et al., 2000) ( where M ak = F ak L a (m) is the volume of sediment of size fraction k in the active layer per unit of bed area, L a (m) is the active layer thickness, and F ak ∈ [0, 1] (−) and I k ∈ [0, 1] (−) are the volume fraction content of size fraction k in the active layer and at the interface between the active layer and the substrate, respectively. Conservation of sediment in the substrate for each size fraction yields (Stecca et al., 2014) ( where M sk (m) represents the volume of sediment in the substrate per unit of bed area. By definition, the sum of the volume fraction content over all size fractions equals 1: The main variables of the SILKE model are shown in Figure 1, together with the active layer model for comparison.
In tracer sediment dispersion, three time scales exist, in which different physical processes play a role (Nikora et al., 2002). The shortest (local) scale considers the movement of particles between impacts with the bed surface. Processes at this scale are not included in the model. The intermediate scale considers processes related to a single entrainment and deposition cycle. The variability in particle velocity, taken into consideration by the particle diffusivity, is a relevant process at this scale. At the global scale, sediment dispersion is governed by the rest time of sediment particles and the interaction between the bed surface and the sediment in transport (Furbish et al., 2017;Martin et al., 2012Martin et al., , 2014Pelosi et al., 2014;Voepel et al., 2013). Nikora et al. (2002) sets the limit time between the intermediate and global scales in t = 15d k ∕u * , where u * = √ b ∕ w (m/s) is the shear velocity, b (N/m 2 ) is the bed shear stress, and w (kg/m 3 ) is the water density.
When modeling unisize tracer sediment dispersion under normal flow conditions, the SILKE model is similar to the one developed by Lajeunesse et al. (2013Lajeunesse et al. ( , 2017Lajeunesse et al. ( , 2018, which is based on the sediment balance model developed by Charru et al. (2004), Charru and Hinch (2006), Charru (2006), and Lajeunesse et al. (2010). The difference is that the SILKE model includes particle diffusion, and we maintain the active layer concept, whereas Lajeunesse et al. (2013Lajeunesse et al. ( , 2017Lajeunesse et al. ( , 2018 model the number of static and moving particles. In section 3.4, we will analytically compare both models.

Closure Relations
In this section we introduce the closure relations for friction, entrainment and deposition rates, and particle velocity and diffusivity. We assume a Chézy-type friction such that S f = b ∕( w gh) = C f Fr 2 , where C f (−) is a nondimensional friction coefficient that we assume to be independent of the flow and bed parameters.
We assume that, for a given bed shear stress, the rate at which particles of size fraction k are set into motion per unit of bed area, E k , depends on the volume fraction content of sediment of size fraction k in the active layer F ak (e.g., Parker, 2008). The rate at which particles of size fraction k deposit per unit of bed area, D k , depends, for a given bed shear stress, on the volume of sediment of size fraction k in transport per unit of bed area, Γ k (e.g., Seminara et al., 2002). The particle activity, Γ k , is nondimensionalized using the grain size, and the entrainment and deposition rates are nondimensionalized using the parameter √ Rgd k (Einstein, 1950;Fernandez-Luque & Van Beek, 1976;Kalinske, 1947;Seminara et al., 2002), where R = ( s − w )∕ w (−) is the submerged specific gravity and s (kg/m 3 ) is the sediment density. We assume s =2,650 kg/m 3 and w =1,000 kg/m 3 . It is convenient to define a capacity of entrainment and deposition (Ê k [1/s] andD k [1∕s], respectively), which depend on the variables of the model M ak and Γ k , for the later mathematical analysis of the model ParametersÊ k andD k expressed in terms of the nondimensional entrainment and deposition rates E * k and D * k areÊ The entrainment rate does not depend on the amount of sediment being transported, which implies that we neglect the entrainment due to collisions of particles in transport with the bed (i.e., collective entrainment; Ancey & Heyman, 2014).
The nondimensional particle velocity v * pk (−) is defined as (Fernandez-Luque & Van Beek, 1976;Seminara et al., 2002) We need to specify closure relations for the dimensionless particle entrainment E * k , deposition D * k , and velocity v * pk . However, the sediment transport rate under equilibrium conditions, q bk eq (m 2 /s), is the variable traditionally related to the flow properties (e.g., Ashida & Michiue, 1971;Engelund & Hansen, 1967;Meyer-Peter & Müller, 1948;Wilcock & Crowe, 2003). From equations (2) and (3), we find that, under equilibrium conditions, E k = D k and q bk = v pk Γ k , respectively. Substituting equations (8)-(10) in the above equations, we find the following relation between entrainment, deposition, velocity, and transport: 10.1029/2019JF005081 where q * bk eq = q bk eq ∕F ak √ Rgd 3 k (−) is the nondimensional equilibrium sediment transport rate. One can decide which three variables are specified using closure relations. We choose E * k , v * pk , and q * bk eq . In principle, any combination of closure relations can be used. However, attention needs to be paid when using closure relations which account for a critical Shields (1936) stress in conditions close to incipient motion, as this may create discontinuities in, for instance, q * bk eq if D * k is zero but not the product E * k v * pk . Fernandez-Luque and Van Beek (1976) developed closure relations for unisize sediment conditions. We generalize their expressions to mixed-size sediment conditions. To this end, we consider the bed shear stress for each size fraction and account for hiding effects. These hiding effects represent the facts that (a) fine fractions in a sediment mixture hide behind large grains and experience a larger critical bed shear stress than under unisize conditions and (b) coarse particles are more exposed to the flow than under unisize sediment conditions (and so experience a smaller critical bed shear stress) (Einstein, 1950). The relations by Fernandez-Luque and Van Beek (1976) extended to mixed-size sediment conditions read where Shields (1936) stress, and k (−) is the hiding function. We consider the hiding function by Ashida and Michiue (1971): where D m (m) denotes the geometric mean grain size of the bed surface sediment.
We also apply the sediment transport relation developed by Ashida and Michiue (1971): The combination of the sediment transport relation by Ashida and Michiue (1971) and the entrainment rate and particle velocity relations by Fernandez-Luque and Van Beek (1976) yields a discontinuous deposition rate close to incipient motion (see equation (11)). The discontinuity is avoided using an entrainment function that depends on the excess bed shear stress in a similar manner as the sediment transport rate. As a result, the deposition rate is a constant multiplying the relation for the particle velocity function (see equation (11)). The entrainment function reads where the value of the constant is chosen such that the deposition rate is equal to the one found by Fernandez-Luque and Van Beek (1976). This decision is arbitrary, and ideally, one would derive the constant based on experimental measurements. We lack the necessary data for this purpose.
The diffusivity in the sediment transport k is equal to the product of the variance of the particle velocity and the Lagrangian integral time scale of the particle velocity, which is a measure of the time over which the particle velocity is significantly autocorrelated (e.g., Nieuwstadt et al., 2016). Using high-speed photography, Lajeunesse et al. (2010), Roseberry et al. (2012), , Furbish and Schmeeckle (2013), and Fathel et al. (2015) show that, at least under laboratory conditions, the probability distribution of particle velocity is well approximated by an exponential distribution. Assuming an exponential distribution of particle velocities, the variance of the particle velocity is equal to the square of the average particle velocity. As particle travel time is of the same order of magnitude as the integral time scale (Roseberry et al., 2012), 10.1029/2019JF005081 particle diffusivity is proportional and of the same order of magnitude as the product of the average particle velocity and the step length (Furbish, Ball, et al., 2012).
Various relations for the mean step length have been proposed. Nakagawa and Tsujimoto (1980b) found that, under mixed-size sediment conditions, the step length is 10-30 times the grain size, and it slightly increases with the bed shear stress. Sekine and Kikkawa (1992) derived a relation depending on the bed shear stress and the sediment fall velocity, and Niño et al. (1994) found that the step length varies between 4.5 and 8 times the grain size. Similar results were obtained by Hu & Hui (1996a, 1996b, who included the effect of sediment density. A recent study by Hosseini-Sadabadi et al. (2019) reviews literature data and explains the scatter from the different definitions for motion and the experimental bias inherent to the finite area where measurements are conducted. In our analysis, only the order of magnitude of the diffusivity matters. For this reason, we choose the simple relation by Fernandez-Luque and Van Beek (1976) in which the mean step length is 16 times the grain size, and we test the sensitivity of our results to this choice. Then, the particle diffusivity is

Numerical Solution and Boundary Conditions
To numerically solve the system of equations of the active layer model, we use the numerical research code Elv as explained in Chavarrías, Stecca, et al. (2019). In brief, we solve the backwater equation in combination with the flux form of the Exner (1920) and active layer equations in a decoupled manner. This implies that the equations are assumed to weakly interact with each other, which is acceptable if the Froude number is below approximately 0.7 (Lyn, 1987;Lyn & Altinakar, 2002;Sieben, 1999). The backwater equation is solved using the standard fourth-order Runge-Kutta method, and the Exner (1920) and active layer equations are solved using a first-order upwind scheme (i.e., FTBS (Forward in Time, Backward in Space); Long et al., 2008;Sonke et al., 2003;Zima et al., 2015) using a variable time step to guarantee a CFL (Courant-Friedrichs-Lewy) number (Courant et al., 1928; see, e.g., Toro, 2009) equal to 0.9. The substrate is discretized using Eulerian cells such that only the cell below the active layer (apart from the active layer itself) has a variable thickness. Cells are created when a specified thickness of this cell is reached under aggradational conditions. The implementation is mass conservative independent of the number of cells created or consumed in one time step .
We have extended Elv to solve the system of equations of the SILKE model. Equation (4) for the bed elevation and equation (5) for the volume of sediment of each size fraction in the active layer are solved using a first-order forward Euler scheme. When the diffusion coefficient of equation (2) for the particle activity is set to 0, the equation is solved using an FTBS scheme. When considering diffusion, the difference in the order of magnitude between the advective and diffusive components limits the application of a second-order centered scheme. In order to reduce the computational time, we solve the advective term using an upwind scheme and the diffusive term using centered differences. The time integration is performed using the scheme proposed by Crank and Nicolson (1947). The source term complicates the implementation of an automatic time stepping method. For this reason, the time step is fixed, and the largest possible value that guarantees a monotone solution is found by trial and error.
Boundary conditions are necessary for the backwater and particle activity equations only. Regarding the backwater equation, we impose the downstream water surface elevation, which agrees with the fact that the flow is subcritical. Regarding the particle activity, we consider two possibilities: a recirculating and a feed flume-like simulation (Parker & Wilcock, 1993). In the case of recirculating flume-like simulation, we impose cyclic boundary conditions on the particle activity by copying the downstream value of the particle activity at the upstream end (i.e., a Dirichlet boundary conditions).
In the case of a feed flume-like simulation, we make a distinction between the cases accounting for diffusion and the cases neglecting it. When diffusion is neglected ( k = 0∀k), information in the particle activity equation travels in the downstream direction only. In this case, only an upstream boundary condition is required to solve the particle activity equation. The variable that is usually set is the sediment transport rate, which, in this case, is equivalent to a condition on the particle activity (i.e., a Dirichlet boundary condition; see equation (3)). Diffusion requires both upstream and downstream boundary conditions to solve the equation. Imposing the sediment transport rate at the upstream end involves a condition on the particle activity and its derivative (i.e., a Robin boundary condition; see equation (3)). It is difficult to impose a boundary condition at the downstream end, as there is no physical constraint regarding the sediment transport rate or the particle activity. As in our simulations, the downstream end of the domain is initially under equilibrium conditions; we impose equilibrium sediment transport rate at the downstream end (i.e., a Robin boundary condition) and ensure that the domain is long enough such that this downstream boundary condition does not affect the solution within the domain of interest.

Linear Model
In this section we linearize the SILKE model to gain insight on its properties and characteristics analytically. The reader is referred to section S1 in the supporting information for the details. To linearize the system of equations, we consider a reference state which is a solution of the system of equations. The reference state represents steady uniform straight flow over a flat sloping bed composed of an arbitrary but uniform grain size distribution. A small perturbation is added to the reference state. As the perturbations are small, the products between perturbed variables (i.e., the nonlinear terms) are even smaller and can be neglected. This step limits the time scale in which the solution of the linear model is valid in the case that perturbations grow. When perturbations achieve a significant value, nonlinear terms play a role in the solution. By substituting the reference solution, we obtain a system of equations of the perturbed variables. Further, we assume that the solution of the system of equations can be expressed in the form of normal modes. Worded differently, the solution is a set of waves with a certain real wave number k w (rad/m) and complex angular frequency = r + i i (rad/s). Variable r is the angular frequency and i the attenuation coefficient.
The solution of the system of equations is now reduced to an eigenvalue problem. This means that we can easily compute whether perturbations with a certain wave number will grow or decay. In the following sections, we solve this eigenvalue problem to gain insight on the model behavior.

Instability Mechanism
In this section we analyze the stability of the model. To this end and following the framework set by Callander (1969) (section 1), we study the growth rate of perturbations using the linear model (section 3.1) as a function of the wave number.
Considering a two sediment size fractions case, we reduce the variables of the linear model to four ( , Γ 1 , Γ 2 , and M a1 ). Furthermore, we neglect the variability in particle velocity, which implies that particle diffusion is zero (i.e., k = 0∀k). The model is still too complex to be solved analytically. We observe that the particle activity equation has an advective character (second term on the left-hand side in equation (2)) which is not present in the equations for the bed elevation (equation (4)) and bed surface texture (equation (5)). Moreover, in the source term of all equations (i.e., the right-hand side of equations (2), (4), and (5)), the term multiplying the particle activity (i.e.,D k ) is larger than the term multiplying the volume of sediment in the active layer (i.e.,Ê k ). While the deposition rate scales with the inverse of the grain size, the entrainment rate scales with the inverse of the active layer thickness. These differences hint at the possibility of two different time scales, one associated with changes in particle activity and one with changes in bed elevation and the bed surface texture. Section S2 indeed confirms that there exist two different time scales and that the time scale of changes associated with particle activity is shorter (changes are faster) than the time scale associated with changes in bed elevation and bed surface texture.
We use the existence of two different time scales to decouple the model and obtain an approximate solution of the fast variables ( f 1 and f 2 ) and the slow variables ( s 1 and s 2 ). Independently of the wave number, the growth rate of the fast variables is negative, which means that the fast variables are stable. The growth rate of the slow variables may be both positive or negative. We identify three main parameters controlling the growth or decay of perturbations. The first parameter is Ψ (−) (section S2): where =Ê 1D2 ∕Ê 2D1 is a parameter relating the entrainment and deposition rates of the two size fractions and V = v p1 ∕v p2 relates the particle velocities. Perturbations are stable if Ψ > 0, and perturbations may grow if Ψ < 0 (section S2). As can be larger or smaller than one depending on the flow conditions, it is not simple to generalize the behavior of Ψ. However, parameter Ψ can be interpreted as a measure of the grain Figure 2. Separatrix between the growth domain (G) and decay domain (D) as a function of the wave number (k w ) and Ψ (equation (19)). The wave number is nondimensionalized using the fine grain size. Parameter Ψ is relative to its minimum value Ψ min = 1∕ − 2V. Each line represents a separatrix as a function of the ratio between (a) the grain size of the fine fraction and the active layer thickness and (b) the grain sizes. In (a), the active layer thickness is varied to obtain different conditions, while in (b), we vary the grain size of the fine size fraction. The green markers represent the conditions of two numerical simulations under growth (circle) and decay (square) conditions. In the numerical simulations ,d 1 ∕L a = 0.01 and d 1 ∕d 2 = 0.2, and the separatrix is highlighted using a dashed black line.
size distribution of the active layer sediment relative to the grain size distribution at the interface between the active layer and the substrate. Parameter Ψ decreases when the active layer is coarse compared to the interface between the active layer and the substrate. Parameter Ψ is small under the conditions in which the active layer model is prone to be ill-posed (i.e., degradation into a fine substrate; Chavarrías et al., 2018;Ribberink, 1987;Stecca et al., 2014). Hence, under the conditions in which the active layer model is ill-posed, the SILKE model may be unstable.
The second parameter is the ratio between the grain size of the fine sediment and the active layer thickness, d 1 ∕L a . This parameter relates the particle activity to the volume of sediment in the active layer, as the first scales with the grain size and the second with the active layer thickness.
The third parameter is the ratio between grain sizes, d 1 ∕d 2 . A ratio close to one indicates conditions close to unisize sediment.
In Figure 2, we present the separatrix between growth (G) and decay (D) of perturbations as a function of the wave number based on the reference case of Table 1. The results are presented as a function of the three dimensionless parameters. Considering the reference case, we vary the volume fraction content of fine sediment at the interface between the active layer and the substrate between 0.5 and 1. In Figure 2a, the active layer thickness varies between 3 times the coarse grain size (representative of plane bed conditions ; Parker & Sutherland, 1990;Petts et al., 1989;Rahuel et al., 1989) and 0.2 times the water depth (representative of a  dune-dominated case Deigaard & Fredsøe, 1978;Lee & Odgaard, 1986). In Figure 2b, the characteristic size of the fine sediment varies between 0.0005 and 0.003 m.
A fining of the sediment at the interface (i.e., an increase in I 1 , which causes a decrease in Ψ) increases the instability domain (i.e., there is a larger range of unstable wave numbers). This is seen in Figure 2 from the fact that, given a set of parameters (i.e., given a certain line in the plot, such as, for instance, the dashed black one, corresponding with d 1 ∕L a = 0.01 and d 1 ∕d 2 = 0.2), a low value of Ψ∕Ψ min is associated with a larger range of values of k w d 1 falling in the growth (G) domain. When the fine sediment at the interface is below a minimum threshold, the model is stable, and perturbations do not grow. For instance, for d 1 ∕L a < 0.01, all perturbations decay (i.e., there is no growth (G) domain) if Ψ∕Ψ min > −0.93, which is equivalent to I 1 < 0.8. This shows that the mechanism underlying growth of perturbations is associated with the origin of ill-posedness in the active layer model. The active layer thickness significantly affects the domain in which perturbations grow (Figure 2a). This corresponds to the active layer model, where the active layer thickness plays a significant role in defining the region in which the model is ill-posed (Chavarrías et al., 2018). A decrease in active layer thickness (i.e., an increase in d 1 ∕L a ) decreases the growth domain, which is consistent with the fact that the active layer model is well-posed if the active layer thickness is sufficiently thin (Chavarrías et al., 2018). The active layer model is well-posed if the active layer thickness is sufficiently thick too (Chavarrías et al., 2018). However, we do not see an eventual decrease in the growth domain for increasing values of the active layer thickness (i.e., decreasing values of parameter d 1 ∕L a ), as, in our reference case, this effect occurs for unrealistically large values of the active layer thickness. Worded differently, from a mathematical perspective, a sufficient increase in the active layer thickness will yield a decrease in the growth domain. However, for the parameter values under consideration, this effect falls outside the domain of physical interest. The domain of instability decreases when the ratio of the two grain sizes tends to one (Figure 2b). This shows that the model is stable under unisize conditions, which supports the fact that the mixed-sediment character is a condition that allows for instability.
In this reference case, perturbations with a wave length below approximately 5 m do not grow. Thus, the length scale of growing perturbations is at least 5 times the flow depth (h = 1 m; see Table 1). We conclude that the minimum length scale of growing perturbations is consistent with the assumption of hydrostatic flow in the derivation of the backwater equation (section 2.1). The minimum length scale also indicates that instabilities are strongly linked to the process of adaptation of the particle activity to changing conditions, as this process occurs over lengths characterized by v pk ∕D k (Bohorquez & Ancey, 2016), which is of order 1.
We run two idealized numerical simulations with the parameters of the reference case (Table 1) to test the results of the linear analysis and to verify that under the conditions in which we predict growth, the solution of the system of equations including nonlinear terms shows growth of perturbations. The domain is 200 m long and is discretized using 0.25 m long cells. The initial bed elevation is formed by a sinusoidal perturbation with an amplitude equal to 0.001 m and wave length equal 10 m superimposed on a constant slope equal to 8.2e−04, which is the equilibrium bed slope Blom, Arkesteijn, et al., 2017). The small value of the amplitude of the perturbation guarantees that the linear solution is initially valid. The active layer is initially composed of coarse sediment only. We use periodic boundary conditions. The two simulations differ only regarding the initial composition of the substrate sediment. In one case, the substrate is composed of fine sediment only and the linear analysis predicts growth (circle in Figure 2), while in the second case, the substrate is composed of 80% fine sediment and the linear analysis predicts decay of the perturbation (square in Figure 2). Figure 3 shows the bed elevation relative to the mean longitudinal bed profile as a function of time. A fine substrate causes oscillations to grow (Figure 3a), whereas a slightly coarser substrate yields decay of perturbations ( Figure 3b).
The numerical results confirm that the model presents an instability mechanism strongly linked to the mixed-size character of the sediment. An active layer coarser than the interface between the active layer and the substrate triggers the formation of waves, as was observed by Chavarrías, Stecca, et al. (2019) in laboratory experiments. The numerical results also confirm the validity of the analytical predictor for instability conditions.
We note that we are not the first authors reporting a one-dimensional instability mechanism linked to the mixed-size character of sediment. Tsujimoto and Motohashi (1989) and Tsujimoto (1989aTsujimoto ( , 1989bTsujimoto ( , 1990Tsujimoto ( , 1989c also found an instability mechanism. In their case the instability was driven by a spatial lag between the sediment transport rate and the bed surface grain size distribution. The phase difference originates from considering the sediment transport rate as the convolution integral of the product of the entrainment rate and probability of step length. This instability mechanism explains the formation of alternating fine and coarse stripes in the streamwise direction, as well as in the transversal direction (Ikeda & Iseya, 1986;Iseya & Ikeda, 1987;Kuhnle & Southard, 1988). In both the SILKE model and the model by Tsujimoto and coauthors, the instability mechanism is intrinsically related to mixed-size sediment; it is one dimensional and predicts the formation of alternate fine and coarse stripes. However, the models represent different processes: Substrate sediment only plays a role in the SILKE model. In the SILKE model, a homogeneous mixed-size sediment bed is stable, and waves do not grow, contrary to the model by Tsujimoto and coauthors. The instability mechanism due to entrainment of fine substrate sediment is only captured by the SILKE model.
The fact that for an increasing wave number, the domain in which perturbations decay increases (Figure 2) indicates that the SILKE model may be well-posed. In the following section, we prove that this is true.

Well-Posedness
In this section we show that the SILKE model considering two sediment size fractions is well-posed. Several processes in nature are unstable, which means that perturbations to an equilibrium situation will grow. For instance, a wide alluvial reach is unstable, as, starting from a flat bed, a pattern of alternate bars will grow (Colombini et al., 1987;Schielen et al., 1993). However, all physically realistic instability mechanisms have a wavelength cutoff below which all perturbations decay. Worded differently, for a sufficiently small perturbation, the equilibrium situation must be stable. Indeed, in our previous example, the width-to-depth ratio needs to be above a minimum value for bars to start growing. If this property is not satisfied, infinitely short perturbations grow, and as the growth rate increases exponentially with decreasing perturbation length, the model has, in essence, no physically realistic solution. Worded differently, the model is ill-posed. For this reason, we study well-posedness of the model by studying the growth rate of infinitely short perturbations. We refer to Chavarrías, Schielen, et al. (2019) for a thorough study of ill-posedness and its relation to hyperbolicity of the system of equations and growth of perturbations.
We first consider negligible diffusion in the sediment transport rate. The growth rates of the two fast variables are always negative regardless of the wave number (section 3.2). This implies that they are stable. For a wave number tending to infinity, we find that the angular frequencies of the slow variables are also negative (section S2). To account for the effect of particle diffusivity, we assume the same value of diffusion for both grain sizes. We find that all growth rates are negative. We conclude that, independently from diffusion in the sediment transport, the SILKE model is well-posed.
Interestingly, for intermediate wave numbers, diffusion increases the instability domain, although in the limit for a wave number tending to infinity diffusion dampens oscillations. Perturbations that decay when diffusion in the sediment transport is not accounted for appear to grow when diffusion is considered. This result is consistent with the results of Bohorquez and Ancey (2015). Using their unisize sediment model, they found that diffusion increases the instability domain.

Tracer Sediment Dispersion
In this section we analytically study the propagation of tracer sediment predicted by the SILKE model. While the active layer model predicts no dispersion of tracer sediment (section 1), we expect that in the SILKE model, the interaction between the sediment at the bed surface and the sediment in transport captures sediment dispersion.
In general, particles may disperse in a superdiffusive, subdiffusive, or normal (Fickian) manner (Havlin & Ben-Avraham, 1987;Metzler & Klafter, 2000). Normal diffusion means that the variance of the particle position scales linearly with time. This behavior arises from collisions between a large amount of particles (Einstein, 1905) and is ubiquitous in nature. If variance does not scale linearly with time, particle diffusion is called anomalous. Anomalous diffusion is superdiffusive or subdiffusive depending on whether the relation between the variance and the time is larger or smaller than linear, respectively (Havlin & Ben-Avraham, 1987;Metzler & Klafter, 2000). As the variance characterizes the celerity at which particles spread, anomalous diffusion implies that dispersion is scale dependent (i.e., it depends on the time scale; Bouchaud & Georges, 1990).
We consider a simplified case of tracer dispersion under unisize and normal flow conditions assuming the active layer thickness to be constant. As the measurements of tracer sediment in experiments and field cases usually provide the total amount of tracer sediment without distinguishing between the sediment at the bed surface and the sediment in transport (Lajeunesse et al., 2018), we introduce the variable c (−) that accounts for the total concentration of tracer sediment: where Γ T = Γ 1 + Γ 2 (m) is the total volume of active particles per unit of bed area and F Γ1 = Γ 1 ∕Γ T (−) is the fraction of active tracer particles.
We find an analytical solution of c at first order that is valid in the long term (section S3): where the subscript 0 indicates the approximate value of the solution. This is an advection-diffusion equation in which the advective velocity is and the diffusivity is When = 0, the diffusivity is not equal to 0. This result shows that the model predicts dispersion of tracer sediment as a result of the interaction between the active layer and the particle activity equations. In particular, tracer sediment diffuses in a normal (Fickian) manner as it propagates in the downstream direction. Long particle travel distances, characterized either by a large particle velocity or a small deposition flux, Figure 4. Bed surface fraction content of tracer sediment 9.17 m downstream from the feeding location. The black dots are measured data (Chavarrías, Stecca, et al., 2019), and the lines are the results of numerical simulations. The dash-dotted, continuous, dashed, and dotted lines present the results using the active layer model (Hirano, 1971), the SILKE model without diffusion, the standard amount of diffusion (equation (18)), and 10 times the standard amount of diffusion, respectively. The different colors represent a different active layer thickness.
increase particle diffusion. The propagation velocity of tracer sediment is smaller than the particle velocity, and it depends on the inverse of the active layer thickness. The dependence on the active layer thickness is consistent with the fact that, for a given situation, a thicker active layer represents conditions on a longer time scale (section 2.1), which includes more infrequent mixing with deeper substrate sediment, slowing down the streamwise propagation of tracer sediment (Ferguson et al., 2002).
An advection-diffusion equation for the propagation of tracer sediment in the long term was also obtained by Lajeunesse et al. (2018) using their sediment balance model. In the following, we compare their resulting equation to the one of the SILKE model. The advective velocity in both models is mathematically equivalent if we consider that the number of static and moving particles in their model is represented by the active layer thickness and the particle activity, respectively, in our model. The diffusivity of the model by Lajeunesse et al. (2018) depends on the square of the particle step length and inversely on the particle travel time.
Our results are mathematically equivalent to theirs if we relate the deposition rate to the inverse of the particle travel time. Although mathematically equivalent, the variables represent different properties, such that each model is most suited for reproducing certain physical processes. For instance, the active layer thickness is representative of the range of elevations at which the sediment interacts with the flow and that are responsible for significant vertical mixing, and this property is not physically equivalent to the number of static particles. Moreover, while in the model by Lajeunesse et al. (2018), the number of static particles is a variable of the model under unsteady conditions; in our case the active layer thickness is a priori assigned (or related to flow and bed properties).
An advantage of the SILKE model is that it explicitly accounts for the effect of porosity. If bed particles are loosely packed (i.e., the porosity is large), there are fewer particles per unit of bed area. This decreases the interaction of moving particles with particles at the bed surface, which increases tracer sediment propagation. Moreover, the SILKE model shows the effect of the variability in particle velocity: It increases the diffusive behavior caused by the interaction between particles at the bed surface and in transport.

Tracer Propagation Without Temporary Burial
In this section we model a laboratory experiment conducted by Chavarrías, Stecca, et al. (2019) to study whether the model captures tracer propagation under conditions in which the effect of temporary burial of sediment due to bedforms is negligible (section 2.1). We also investigate the role played by diffusion in sediment transport. The experiment consisted of feeding painted (i.e., tracer) unisize sediment with characteristic grain size equal to 0.0055 m under equilibrium conditions. It was conducted in a 14 m long, 0.4 m wide tilting flume in the Delft University of Technology. The slope was equal to 3.50e−3, and the flow depth was equal to 0.187 m. A regular pattern of small bedforms approximately 2 cm high covered the bed surface. Measurements of bed elevation fluctuations confirm that the active part of the bed was restricted to a narrow range of elevations. Worded differently, the probability that a particle is entrained and deposited decreases fast at lower elevations, such that temporary burial due to bedforms is negligible. The content of tracer sediment at the bed surface at the downstream end of the flume was measured using a submerged camera. Section S4 explains further details of the experiment.
In modeling the experiment, we use the same closure relations and parameters as the ones used by Chavarrías, Stecca, et al. (2019). The sediment transport rate at capacity is computed using the relation by Ashida and Michiue (1971). The domain of interest is 10 m long. We add 40 m at the downstream end to guarantee that the downstream boundary condition for the particle activity does not affect the domain of interest (section 2.3). The domain is discretized using 0.01 m long cells.
In Figure 4, we compare the measured content of tracer sediment to the results predicted using the active layer model and the SILKE model. The active layer model predicts that tracer propagates as a front and does not capture the slow increase in tracer content at the bed surface. The deviation of the solution from a step function is due to numerical diffusion. When considering an active layer thickness equal to 0.01 m, which is reasonable based on measurements of the bed elevation fluctuations, the overall rate of adaptation of the bed surface fraction content is captured by the SILKE model. A value of the active layer thickness equal to 0.02 and 0.005 m consistently predicts tracer sediment to be transported slower and faster than the measured data show, respectively. Although the general trend due to the exchange of particles between the bed surface (i.e., in the active layer) and the particles in motion (i.e., the particle activity) is captured by the SILKE model, the short-term fluctuations in tracer content of the order of minutes due to individual bedforms is not captured. This is because the dynamics due to individual bedforms are resolved neither by the active layer model nor the SILKE model. However, deep burial of sediment at elevations significantly lower than the lower limit of the active layer was negligible, such that temporary trapping of sediment can be neglected and the active layer captures the dynamics of the system.
Diffusion in sediment transport does not play a significant role in the dispersion of tracer sediment ( Figure 4). This is explained by the fact that the time scale at which we are studying dispersion (5 hr) is the global scale (Nikora et al., 2002; section 2.1). The limit time between the intermediate and global scales is of the order of 1 s for the specifications of the present experiment. Dispersion of tracer sediment at the global time scale is captured by the interaction between the active layer and the particle activity equations regardless of the diffusion in sediment transport. A value of the diffusion coefficient 10 times larger than according to equation (18) has hardly any effect on the solution. This shows that the uncertainty in the step length (section 2.2) has no implications as regards to the effect of particle diffusion for the time scales we are considering.

Tracer Propagation With Temporary Burial
In this section we study tracer propagation under conditions with mixed sediment and in which the effect of temporary burial of sediment due to bedforms may not be negligible. The active layer model predicts no dispersion of tracer sediment under unisize conditions (sections 1 and 4.1). Under mixed-size sediment conditions, the active layer model does predict sediment dispersion due to mobility differences between grain size fractions. Contrary to the active layer model, the SILKE model predicts sediment dispersion even without considering diffusion in the sediment transport rate (sections 3.4 and 4.1). Here our objective is twofold. First, we study whether the mixed-size sediment character is enough to capture tracer dispersion using the active layer model. Second, we study whether the SILKE model is capable of reproducing the propagation of tracer dispersion under conditions in which the effect of temporary burial of sediment may play a significant role. To this end, we consider the study conducted by Sayre and Hubbell (1965). They tracked the propagation of tracer sediment along a 550 m long stretch of the North Loup River (Nebraska, USA) during approximately 12 days. The minimum and maximum sizes of the bed sediment were equal to 0.088 and 9.4 mm, respectively, with a geometric mean grain size equal to 0.315 mm. The tracer sediment had a narrower grain size distribution characterized by three grain sizes equal to d 1 = 0.209 mm, d 2 = 0.296 mm, and d 3 = 0.418 mm.
The bed was covered with dunes between 0.30 and 0.46 m height, and based on core samples, the average over all samples of the maximum depth at which tracer was found in each sample was 0.44 m below the bed surface. The flow conditions did not significantly vary during the measurement campaign, and a reasonably 10.1029/2019JF005081 constant flow depth was measured. Figure 5 presents the measured concentration of tracer sediment at the end of the experiment. Details about the conditions of the experiment are presented in section S5.
We select appropriate closure relations and a reasonable value for the skin friction coefficient and active layer thickness to model the field case. We find that the sediment transport relation by Ashida and Michiue (1971) in combination with the skin friction predictor by Engelund and Hansen (1967) performs best with an error of only 3% (section S5). For simplicity, we assume that normal flow prevailed in the study area. We assume that the active layer thickness is equal to the average of the maximum depth below the bed surface at which tracer was found. We discretize the domain using 0.5 m long cells.
We run two numerical simulations using the active layer model to test whether the mixed-size sediment character is enough to capture tracer sediment dispersion. In the first simulation, we consider unisize sediment with a characteristic size equal to the geometric mean grain size of the bed surface sediment. In the second one, we approximate the grain size distribution of the bulk mixture using 11 grain size fractions (which is the number of sieves reported by Sayre & Hubbell, 1965). The tracer sediment had a narrower grain size distribution than the parent material and is represented by three grain size fractions (section S5).
In Figure 5a, we present the total (i.e., considering the mass in both the active layer and in transport) tracer concentration at the final time assuming unisize sediment. The small dispersion in Figure 5a is entirely due to numerical diffusion. In Figure 5b, we observe that accounting for mixed-size sediment causes dispersion. We see three peaks corresponding to the three fractions which, in total, cover a longer river stretch than under unisize conditions. We run two simulations under the same conditions as the ones with the active layer model but using the SILKE model (Figures 5c and 5d). Dispersion as predicted by the SILKE model is larger than when using the active layer model. However, for all cases, the amount of dispersion is an order of magnitude smaller than the measured one. Moreover, the model predicts that tracer sediment disperses in the form of normal (i.e., Fickian) diffusion (section 3.4), whereas the field data suggest that the downstream tail spreads faster than normal (i.e., it decreases linearly in logarithmic scale).
One possible explanation for the discrepancy between the measured and the predicted results in Figure 5 is that Sayre and Hubbell (1965) explain that the measurements of the right-hand side of the river after approximately 120 hr could be affected by the late reentrainment of tracer sediment placed on the left-hand side (see section S5 for details about the measurements). It appears to be that sediment placed on the left-hand side was initially trapped at the location where it was dumped and only later reentrained, reaching the right-hand side. This two-dimensional effect cannot be captured using our one-dimensional model. To study whether the late reentrainment of trapped sediment is responsible for the disagreement between the data and the modeling results, we compare the measured concentrations along the last longitudinal traverse that was unaffected by the reentrainment to the results of the model ( Figure 6). The predicted peak concentration is six times larger than the measured concentration, and the amount of predicted dispersion is significantly smaller than the measured dispersion. We conclude that the late reentrainment of sediment is not responsible for the disagreement between the measured and predicted results.
A second explanation is that the imposed diffusivity in sediment transport is too small. This may be due to the assumption of an exponential distribution of particle velocities or a too small step length (section 2.2). Particle velocities may be better approximated by a different probability distribution that causes a larger diffusivity. Simulations using larger diffusivity constants show that the peak concentration would be captured if the diffusion coefficient is 1,000 times larger than our original value ( Figure 6). Apart from the fact that this amount of diffusion is not physically realistic, the downstream tail is not captured independently of the diffusion coefficient. For this reason, we conclude that the disagreement is not due to a low diffusivity.
We conclude the dominant mechanism responsible for sediment mixing is not captured by the combination of the active layer and the particle activity equations. In particular, this mechanism may be temporary burial due to bedforms (section 2.1). The assumption that the part of the bed that interacts with the flow is well represented by an active layer becomes questionable when large rest times occur due to temporary sediment burial by bedforms (Blom, 2008;Ribberink, 1987).
It seems reasonable that temporary burial due to bedforms explains the fact that part of the tracer sediment propagates slower than predicted, but it does not seem to explain the fact that part of the sediment propagates faster. However, if the effect of temporary burial of sediment is accounted for, the active layer thickness should be reduced, as mixing due to burial previously accounted for by an increased active layer thickness is then accounted for by the burial process. The reduced active layer thickness would cause a faster propagation of tracer sediment.

Ill-Posed Conditions Assuming Constant Active Layer Thickness
The fact that the SILKE model is always well-posed is particularly important when reproducing a situation of degradation into a substrate finer than the bed surface, as under these conditions, the active layer model is prone to be ill-posed (Chavarrías et al., 2018;Ribberink, 1987;Stecca et al., 2014). In this section we model Experiment I4 conducted by Chavarrías, Stecca, et al. (2019) and compare the results to predictions using the active layer model, the regularized model (Chavarrías, Stecca, et al., 2019), and the SILKE model. Experiment I4 consisted of a 4 m long patch of fine sediment (d 1 = 2.1 mm) placed below a 3 cm layer of coarse sediment (d 2 = 5.5 mm). Degradational conditions were imposed by lowering the downstream water level 8 cm during 8 hr. All other parameters were the same as in the experiment considered in section 4.1. Figure 7a shows the bed elevation data. Before fine sediment from the patch was entrained, the bed was covered by bedforms. This is more clearly illustrated by the data upstream of the patch for all time, where bed elevation variation in space for any fixed time is more visible. As degradation proceeded, fine sediment from the substrate started to be entrained at the troughs of bedforms. The entrainment of fine sediment caused the bedforms to grow. Note, for instance, how a deep trough starts to form right at the upstream end of the patch after approximately 3.5 hr. Coarse sediment was deposited in the troughs, which prevented subsequent entrainment of fine sediment. Fine sediment was again available for entrainment as degradation continued, which caused the periodic formation of large bedforms superimposed on the original ones. We refer to Chavarrías, Stecca, et al. (2019) for an in-depth analysis of this experiment. (a) measured, (b) predicted using the active layer model (Hirano, 1971), (c) predicted using the regularized active layer model (Chavarrías, Stecca, et al., 2019), and (d) predicted using the SILKE model. Dashed lines indicate the position of the patch of fine sediment.
The active layer model is ill-posed when reproducing this experiment, and the consequences are clearly seen in the predicted bed elevation (Figure 7b). An unrealistically large oscillation develops as soon as the interface between the active layer and the substrate is composed of fine sediment. Moreover, the solution is grid dependent, and simulations not shown here using a smaller grid size generate ever larger oscillations (see Chavarrías et al., 2018, andChavarrías, Schielen, et al., 2019). The regularized model provides more satisfactory results (Figure 7c). The solution is well-posed and captures the changes in mean bed elevation, but it does not capture the oscillatory behavior when large bedforms entrain fine substrate sediment (Chavarrías, Stecca, et al., 2019). In contrast, the SILKE model shows oscillatory behavior (Figure 7d). The solution is stable upstream of the patch, and waves start to form at the upstream edge of the patch. Fine sediment from the substrate is transferred to the active layer and entrained. The particle activity adapts in the streamwise direction, which causes an imbalance between erosion and deposition and induces oscillations superimposed on the overall degradational trend.
Oscillatory behavior is also observed in the volume fraction content of coarse sediment at the bed surface ( Figure 8). The bed is composed of coarse sediment only, and the growth of large bedforms that lead to the entrainment of fine sediment is associated with a fining of the bed surface ( Figure 8a). The entrainment cycles are less clear downstream from the patch, as fine sediment entrained at the patch mixes with coarse sediment while traveling downstream.
The result of the active layer model is meaningless. It predicts that the bed surface is composed of coarse sediment only. This is because, after the large nonphysical oscillation, coarse sediment from upstream fills the space and forms a new coarse substrate. Hence, further degradation does not entrain fine sediment. The regularized model predicts a constant bed surface grain size distribution as a function of time at each location after a short period of adjustment at the beginning of the run (not visible in Figure 8). The bed surface becomes finer with increasing streamwise position, as fine sediment is entrained from the patch. The instability mechanism of the SILKE model induces oscillatory behavior in the bed surface grain size distribution. However, the variation in bed surface volume fraction content is smaller than the measured one. Similarly to the field tracer case (section 4.2), this indicates that our model lacks some mechanisms responsible for the fast fining and coarsening that we will discuss in section 5.3.
The SILKE model shows a clear advantage with respect to the regularized model, as it contains an instability mechanism responsible for mixing which is not captured by the regularized model. Even though the SILKE model explains more physical processes than the regularized model, both models are well-posed and reproduce the overall degradational trend. For this reason, the regularized model may yield a satisfactory solution if only the general trend is to be captured. However, as the regularized model requires that the active layer thickness is constant with time (Chavarrías, Stecca, et al., 2019), there are cases that cannot be modeled using the regularized model.

Ill-Posed Conditions Assuming Variable Active Layer Thickness
In this section we reproduce a laboratory experiment that cannot be reproduced using either the active layer model (Hirano, 1971) or the regularized active layer model (Chavarrías, Stecca, et al., 2019). We consider Experiment B2 conducted by Blom et al. (2003). Blom et al. (2003) studied vertical mixing of sediment due to dune growth using a trimodal sediment mixture. Barchan dunes migrating over a coarse sediment layer characterized the starting conditions of the experiment. Below the coarse layer, a fine substrate was present, which was not exposed to the flow. The water discharge was increased, which mobilized the coarse layer, led to the entrainment of fine substrate sediment, and caused dunes to grow adapting to a new equilibrium situation. Sediment was recirculated, and normal flow conditions were maintained during the entire experiment. The details of the experiment are described in section S6.
As normal flow conditions were maintained, the mean bed elevation (averaged over the passage of several bedforms) remained constant with time. For this reason, in modeling the experiment using the active layer model, the only source of mixing in the active layer is the flux of sediment from the substrate due to a lowering of the interface between the active layer and the substrate. The lowering of the interface occurs due to an increase of the active layer thickness that represents the increase in dune height.
We have tested the mathematical character of the active layer model using the approach by Chavarrías et al. (2018), and we find that when fine substrate sediment is entrained, the active layer model is ill-posed. We run a numerical simulation of the 50 m long flume discretized using 0.02 m long cells to study the consequences of modeling this experiment using the active layer model. The substrate is discretized using 0.005 m thick cells, and the measured initial grain size distribution is used as initial condition of the substrate. We impose a linear increase of the active layer thickness with time between half the initial and final values of the mean dune height during the first 4 hr of the 25 hr experiment, as this is the approximate period of time over which dunes adapted to the new equilibrium value (Blom et al., 2003). Based on the initial and final values of the sediment transport rate, flow conditions, and volume fraction content in the active layer, we select a sediment transport relation which reasonably approximates the initial and final conditions. We find that the relation by Meyer-Peter and Müller (1948) using the hiding correction by Parker et al. (1982) with a power parameter equal to 0.9 yields a reasonable approximation.
In Figure 9, we present the volume fraction content in the sediment transport rate as a function of time at the downstream end of the flume. The coarsening of the load over the first hour (note the increase in volume fraction content of coarse sediment and decrease of the volume fraction content of fine sediment) is due to the entrainment of the coarse layer underneath the migrating barchan dunes. This is captured by the active layer model. Subsequently, fine sediment in the substrate starts being transferred to the active layer, as the active layer thickness continues to grow. The fining of the active layer causes a fining of the load. As the sediment at the interface between the active layer and the substrate becomes fine, the active layer model becomes ill-posed. This causes oscillations in the solution that are not visible in the figure, as the oscillations due to ill-posedness are not instantly felt at the downstream end. A large wave develops in the center of the domain and causes entrainment of a large amount of fine sediment causing the bed surface to consist of fine sediment only. The wave is felt at the downstream end after approximately 5 hr (note the sudden fining of the load). The solution breaks down after 10 hr, as the flow solver is incapable of dealing with such abrupt changes in bed elevation. Test simulations not shown here exhibit the unrealistic behavior expected from ill-posed simulations: A decrease of the grid size causes the large wave to form faster and arrive downstream earlier, which causes the simulation to break down earlier too.
We simulate the same conditions using the SILKE model without accounting for the effect of diffusion in sediment transport, and we find a stable solution (Figure 9). Initially, the solution is very similar to the one predicted using the active layer model. This is reasonable, as variations in particle activity occur at a shorter time scale than changes in the grain size distribution (section 3.2). Due to the relatively slow increase in active layer thickness, the particle activity is almost at capacity conditions. However, although the effect of accounting for the particle activity seems negligible initially, the consequences become striking after 4 hr. The simulation does not lose its well-posed character, and after the active layer thickness becomes constant with time, there no longer is a flux of sediment between the substrate and the active layer, and the grain size distribution of the active layer becomes constant with time too. For this reason, the grain size distribution of the load becomes constant with time. Figure 9. Volume fraction content of fine, medium, and coarse grain size fractions in the transported sediment as a function of time in Experiment B2 conducted by Blom et al. (2003). The dots, lines, and dot-dashed lines correspond to the measured data, SILKE model, and active layer model, respectively. The color of lines and dots is related to the grain size. The fine, medium, and coarse size fractions are shown in blue, red, and green, respectively. The active layer model and SILKE model predict an almost equal result during the first 4 hr, and the lines are superimposed. The active layer model crashes after approximately 10 hr.
In conclusion, it is crucial to account for the interaction between the bed surface and the sediment in transport occurring at "short" time scales to obtain a well-posed model that captures the dynamics of the "long" time scale. Although the effect of accounting for particle activity may seem negligible at first instance, as changes in the volume fraction content of sediment in the active layer are slow, we have to consider the small short time scale to obtain a physically realistic model at long time scales.

Discussion
In this section we discuss the reasons why the SILKE model is well-posed (section 5.1). Subsequently, we discuss the limitations of the model (section 5.2) and of the model's instability mechanism (section 5.3). We provide remarks on the modeling of the experiment by Blom et al. (2003) in section 5.4.

Physical Reasoning of the Well-Posedness of the Model
In section 3.3, we have shown that the SILKE model is well-posed, a property shared with the regularized model (Chavarrías, Stecca, et al., 2019). Here we clarify why, although being totally different approaches, both the SILKE model and the regularized model yield well-posed problems.
We consider a situation in which the active layer model is ill-posed such as the laboratory experiment in section 4.3. The active layer model becomes ill-posed because it is incapable of reproducing the "short" time scale phenomena related to the entrainment of fine substrate sediment. In the active layer model, the sediment transport is assumed to be at capacity conditions, and this results in a too crude approximation of the cycles of fine sediment entrainment.
The regularization strategy slows down the mixing processes (Chavarrías, Stecca, et al., 2019), which is mathematically equivalent to considering a thicker active layer. A thick active layer represents changes over a "long" time scale that encompasses the effect of large and less frequent bed elevation fluctuations (section 2.1). A slowdown of the mixing processes by the regularization strategy yields a well-posed model, as "short" term processes of entrainment and deposition are filtered out and not resolved.
The SILKE model resolves the "short" time scale by accounting for conservation of the sediment in transport. Sediment transport is not assumed to be at capacity, and this allows for the formation of oscillations that enhance vertical mixing of sediment. Solving for this physical process yields a well-posed model.
In summary, in both cases, the regularized model and the SILKE model yield a well-posed model by considering the processes occurring at a "short" time scale. The regularized model filters the processes, while the SILKE model resolves them.

Limitations of the SILKE Model
In developing the SILKE model, we combine and extend previous modeling strategies. We maintain the active layer concept of Hirano (1971) to capture mixed-size sediment processes at the bed surface. In dealing with the sediment in transport, we consider the variability in particle velocity (Furbish, Haff, et al., 2012;Ancey & Heyman, 2014;Bohorquez & Ancey, 2015). Nevertheless, in our generalization to mixed-size sediment conditions, we neglect the correlation terms that appear in the particle activity equation (Furbish, Haff, et al., 2012) (section 2.1). There is a need to derive a particle activity equation which considers the effect of the correlation terms. This, presumably, would affect both the advective and diffusive components of the particle activity. As we have seen, if one is interested in the "long" time scale, the diffusive component does not significantly affect the solution. For this reason, we do not expect that sediment dispersion at a "long" time scale differs considerably when taking into consideration the correlation terms if these only affect the diffusive component. However, the advective component may also be affected by the correlation terms, and this could have a significant influence.
There is no specific layer to which the active particles are constrained. This implies that the model may account for both bed load and suspended load. However, our closure relations account for bed load only, and the cases we have considered are bed load dominated. In accounting for both bed load and suspended load, it suffices that the closure relations for entrainment, deposition, particle velocity, and diffusion reflect the different transport mode. In contrast, for modeling suspended load using the active layer model, not only the closure relations need to be adapted, but also the active layer model needs to be combined with a model that accounts for transport of suspended matter.
In the field, large rest times are measured due to sediment being trapped in bars, pools, and floodplains (Ferguson et al., 2002;Hassan et al., 1991;Malmon et al., 2003;Pyrce & Ashmore, 2003;Schmidt & Ergenzinger, 1992). Nakagawa and Tsujimoto (1979) argue that under plane bed conditions, the rest period and step length follow an exponential distribution, and ripples and dunes cause them to be nonexponential, as the rest period depends on the elevation at which particles are deposited. The lowest elevation at which tracer sediment is found increases with time, as the probability of a large bedform that deposits sediment lower in the substrate increases (Crickmore & Lean, 1962a, 1962b. This process also decreases the lowest elevation at which tracer sediment is found in the streamwise direction (Galvin, 1965;Hassan et al., 1999). These effects are not included in the model, as all sediment in the active layer has the same probability of being entrained and there is no sediment flux to the substrate under normal flow conditions.
The model does not have a mechanism to account for large resting times due to the temporary burial of sediment by bedforms. This fact causes the model to not properly reproduce tracer dispersion under conditions in which the probability of bed elevation fluctuations is not well represented by a step function (i.e., by the active layer; section 4.2). The existence of large rest times is related to the manner in which sediment disperses. Laboratory and field data suggest that tracer sediment diffuses anomalously (Bradley, 2017;Bradley et al., 2010;Fan et al., 2017;Ferguson et al., 2002;Haschenburger, 2011Haschenburger, , 2013Nikora et al., 2002), whereas the model predicts normal (i.e., Fickian) diffusion. The propagation of tracer sediment is best represented by means of a fractional advection-diffusion equation (Ancey, 2010;Martin et al., 2012;Schumer & Jerolmack, 2009), which originates, for instance, from heavy-tail resting times under unisize sediment conditions (Fan et al., 2016), as well as from a relation between the travel distance and grain size under mixed-size sediment conditions (Ganti et al., 2010;Hill et al., 2010).
To model the increased complexity and the possibility of large rest times, it may be necessary to abandon the active layer model embedded in the SILKE model in favor of a continuous formulation of the bed sediment such as that by Pelosi et al. (2016). The latter predicts anomalous diffusion in the propagation of tracer sediment (i.e., scale dependence) even considering thin-tailed sediment statistics, as it models the trapping of sediment deep in the substrate. The challenge is to combine this formulation with mixed-size sediment, bed elevation change, and mass conservation of the sediment in transport. Recently, Wu et al. (2019) captured anomalous diffusion using an adapted version of the active layer model that accounts for permanent burial of tracer sediment by means of an extra sink term in the active layer equation. This line of research may allow for the possibility of preserving the essence of the active layer model in improving our understanding of morphodynamic processes. A following step may be to consider the effect of a variable active layer thickness due to, for instance, changes in bedform geometry, on tracer dispersion. An analysis similar to the one conducted in section 3.4 including this effect may yield anomalous diffusion as observed in the field.
Finally, we have studied the dynamics under one-dimensional conditions only. The model can be generalized to two-dimensional conditions by considering particle velocity as a vector quantity and particle diffusion to be different in the streamwise and transverse components. We do not expect well-posedness to be affected by the extension to two-dimensional conditions. On the other hand, we expect the dispersion of tracer sediment in the streamwise direction to be affected by the fact that sediment is dispersed in the transverse direction. It would be interesting to repeat the exercise of modeling the experiment by Sayre and Hubbell (1965) using a two-dimensional version of the SILKE model.

Limitations of the Interpretation of the Model Instability
In section 3.2, we found that the SILKE model presents an instability mechanism, and in section 4.3, we related the mechanisms to the formation of waves in the experiments conducted by Chavarrías, Stecca, et al. (2019). In this section we discuss limitations in the interpretation of the mathematical instability.
One is the linearization of the model. The solution is only valid for small perturbations and short times such that nonlinearities do not play a significant role. In this sense, as in all other linear studies, the main outcome of the analysis is the fact that it predicts the initial formation of waves but not their evolution. A second limitation specific to our mathematical analysis is the fact that we assume that the grain size distribution of the sediment at the interface between the active layer and the substrate is constant and equal to the one of the substrate. This is valid under degradational conditions or when the active layer grows. However, as soon as a wave is formed, there exist aggradational (i.e., the lee face) and degradational (i.e., the stoss face) zones. In the aggradational zone, the grain size distribution at the interface between the active layer and the substrate is a combination of the grain size distribution of the sediment in the active layer and in transport (Hoey & Ferguson, 1994;Toro-Escobar et al., 1996). For this reason, the formation of a wave coarsens the substrate sediment. This substrate coarsening mechanism reduces the time over which our linear solution is valid. Nevertheless, numerical simulations confirm the persistence of the instability mechanism at the nonlinear level.
We have assumed steady flow in the analysis. This assumption is reasonable for long waves (i.e., wavelength of several times the flow depth), but it may be interesting to study whether the result is significantly different when considering unsteady flow. Assuming unsteady flow would not only yield a more accurate result for the cases we have considered but also increase the applicability range to transcritical conditions (0.8 ⪅ Fr ⪅ 1.2), in which it is not valid to assume steady flow due to the strong interaction between the bed and the flow. This would allow for studying, for instance, the effect of mixed-size sediment in antidune formation.
We have modeled an experiment conducted under conditions in which the active layer model is ill-posed and the active layer thickness can be assumed to be constant (section 4.3). The model predicts the formation of waves that we relate to the oscillatory behavior measured in the laboratory experiment. However, there are several differences. The model oscillations induce less fining and at a different time scale than it was observed. This is explained by the fact that, in the laboratory experiment, the formation of waves was strongly linked to bedforms. Bedforms developed upstream of the patch of fine sediment (i.e., under unisize conditions) grew when fine sediment was entrained at the troughs. The SILKE model does not predict the formation of bedforms, which implies that the physical mechanism of entrainment of fine sediment at the troughs may be similar but not completely equivalent to the instability causing wave formation in the model.
Moreover, there is also a difference between the waves predicted in the numerical simulation of the experiment and the ones predicted by the linear analysis. In the numerical simulation of the laboratory experiment, waves form at a discontinuity in the grain size distribution of the substrate sediment in the streamwise direction because of the presence of a patch of fine sediment. This mechanism of wave formation is related but not equivalent to the one predicted in the linear analysis, in which the substrate is uniform. As the predicted wave length is of the order of several meters, we would need laboratory experiments in a longer flume to capture these waves. Such laboratory experiments would shed light on the evolution of oscillations for long time scales and the validity of the instability domain.

Modeling of the Experiment by Blom et al. (2003)
One of the reasons to develop the SILKE model is that there are conditions, such as the ones in the experiment conducted by Blom et al. (2003), that cannot be reproduced using either the active layer model (Hirano, 1971) or the regularized model (Chavarrías, Stecca, et al., 2019). However, Blom (2008) modeled this experiment using the active layer model and compared it to the results of the model developed by Blom and Parker (2004) and Blom et al. (2006Blom et al. ( , 2008, in which the bed is treated in a continuous fashion and specifically accounts for mixing due to dunes. This was possible because, by considering normal flow conditions, Blom (2008) neglected all spatial derivatives. As there is no dependence on the spatial coordinate, the active layer equation is no longer a partial differential equation, and it does not suffer from ill-posedness.
The continuous model by Blom and Parker (2004) and Blom et al. (2006Blom et al. ( , 2008 satisfactorily reproduced the experimental data, yet spatial changes are not accounted for, and an instability mechanism cannot be captured. It would be interesting to model spatial changes using the continuous formulation under conditions in which we predict growth of perturbations. The model by Blom and Parker (2004) and Blom et al. (2006Blom et al. ( , 2008 may also predict an instability mechanism and capture the formation of waves that contribute to the vertical mixing caused by dunes. If this mechanism is not present, the continuous formulation could be extended relating the dune size to the grain size distribution of the bed surface. This would increase the vertical mixing when dunes migrate over a fine substrate.

Conclusions
The active layer model (Hirano, 1971) accounting for mixed-size sediment morphodynamic processes can be ill-posed. For this reason, there are problems which cannot be modeled using it. Chavarrías, Stecca, et al. (2019) derived a regularized active layer model that can deal with situations in which the active layer model is ill-posed. However, the regularized active layer model cannot deal with situations in which the active layer thickness changes with time, for instance, when representing the vertical mixing due to changes in dune height. We derive the SILKE model by combining the active layer concept for mass conservation of sediment at the bed surface and the particle activity concept (Ancey & Heyman, 2014;Furbish, Haff, et al., 2012) for mass conservation of transported sediment. The SILKE model is unconditionally well-posed. For this reason, using the SILKE model, we can reproduce situation which we could not before.
A laboratory experiment conducted under conditions in which the active layer model is ill-posed shows that the situation is unstable (Chavarrías, Stecca, et al., 2019). This instability mechanism is captured by the SILKE model only.
Contrary to the active layer model, the SILKE model predicts tracer sediment to advect and diffuse. Moreover, the SILKE model accounts for different grain sizes, contrary to other models predicting sediment dispersion. These properties allow us to satisfactorily model tracer propagation under conditions in which the effect of large rest times caused by temporary burial due to bedforms is negligible. However, when the effect of temporary burial of sediment becomes the dominant mechanism, the model yields an underprediction of sediment dispersion. Moreover, dispersion may be anomalous (non-Fickian), while our model predicts normal diffusion. Thus, although we increase the range of processes that can be modeled using the SILKE model, there is still need for including certain physical mechanisms to achieve a model applicable over a wider range of situations. Netherlands Organization for Scientific Research (NWO) and which is partly funded by the Ministry of Economic Affairs under Grant P12-14 (Perspective Programme). The numerical model Elv can be downloaded online (https://doi.org/10. 5281/zenodo.2595222).