Ill-posedness in modeling mixed sediment river morphodynamics

In this paper we analyze the Hirano active layer model used in mixed sediment river morphodynamics concerning its ill-posedness. Ill-posedness causes the solution to be unstable to short-wave perturbations. This implies that the solution presents spurious oscillations, the amplitude of which depends on the domain discretization. Ill-posedness not only produces physically unrealistic results but may also cause failure of numerical simulations. By considering a two-fraction sediment mixture we obtain analytical expressions for the mathematical characterization of the model. Using these we show that the ill-posed domain is larger than what was found in previous analyses, not only comprising cases of bed degradation into a substrate finer than the active layer but also in aggradational cases. Furthermore, by analyzing a three-fraction model we observe ill-posedness under conditions of bed degradation into a coarse substrate. We observe that oscillations in the numerical solution of ill-posed simulations grow until the model becomes well-posed, as the spurious mixing of the active layer sediment and substrate sediment acts as a regularization mechanism. Finally we conduct an eigenstructure analysis of a simplified vertically continuous model for mixed sediment for which we show that ill-posedness occurs in a wider range of conditions than the active layer model. © 2018 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The mixed character of the sediment is a property necessary to explain physical phenomena such as downstream fining ( Sternberg, 1875;Blom et al., 2016 ), the gravel sand transition zone ( Yatsu, 1955;Blom et al., 2017 ), the formation of bedload sheets ( Seminara et al., 1996 ), or bed surface armoring . Hirano (1971) was the first to develop a mass conservation model for mixed-size sediment. The model assumes that the topmost part of the bed, i.e. the active layer, interacts with the flow and is instantaneously mixed. Below the active layer lies the substrate which can have vertical stratification. In this schematic representation of the morphodynamic processes only the active layer sediment is affected by entrainment and depositional processes. A vertical flux of sediment originates from changes in elevation of the interface between the active layer and the substrate.
One of the critical aspects of the active layer model is the fact that the vertical extent of the active layer, or active layer thickness, shall be a priori assigned. However, it cannot be physically measured, as it stems from the above schematic representation ( Siviglia * Corresponding author. E-mail address: v.chavarriasborras@tudelft.nl (V. Chavarrías). et al., 2017;Church and Haschenburger, 2017 ). The active layer thickness is related to the time scale of the process under consideration ( Bennett and Nordin, 1977;Rahuel et al., 1989;Sieben, 1997;Wu, 2007 ). In plane bed conditions and short time scales the active layer thickness is assumed to be proportional to the size of a characteristic coarse fraction in the bed, for instance, D 84 or D 90 (e.g., Petts et al., 1989;Rahuel et al., 1989;Parker and Sutherland, 1990 ). If bed forms are predominant and the time scale under consideration involves the mixing induced by the passage of several bed forms, the active layer thickness is typically related to a characteristic bed form height (e.g., Deigaard and Fredsøe, 1978;Lee and Odgaard, 1986;Armanini and Silvio, 1988 ). The active layer thickness may vary over space and time, although often it is assumed to be a uniform constant.
The active layer modeling framework has proven to be able to represent a wide variety of physical phenomena such as bed surface armoring (e.g., Park and Jain, 1987 ) and the morphodynamics of gravel-bed rivers (e.g., Vogel et al., 1992 ) and tidal basins (e.g., Carniello et al., 2012 ). Moreover, it is implemented in a large amount of software packages such as Telemac ( Villaret et al., 2013 ), Delft3D ( Sloff and Mosselman, 2012 ), and BASEMENT ( Vetsch et al., 2006 ). In the latter case, a perturbation in bed elevation introduces another wave, which is mainly related to the bed surface grain size distribution. Yet, each wave perturbs the flow, bed elevation, and bed surface grain size distribution. The arrows indicate the direction of propagation of the perturbations under subcritical flow conditions. The words "water", "bed", and "sorting" refer to a perturbation in water flow, bed level, and surface grain size distribution, respectively.
The mathematical representation of river morphodynamics should be well-posed. This means that the mathematical problem must have a unique solution which depends continuously on the data ( Hadamard, 1923 ). If the solution does not depend continuously on the data, the model is unfit to represent the corresponding physics.
Despite its widespread use, the active layer model has one major mathematical shortcoming: the model can change its mathematical character under some parameter settings. Therefore the mathematical problem that represents the physics of river morphodynamics can become ill-posed. This fact was first recognized by Ribberink (1987) . To this end he simplified the active layer model by considering an equation for the mean grain size of the active layer sediment rather than one active layer equation for each grain size fraction. He found that under aggradational conditions the problem is unconditionally well-posed and the system may become ill-posed under degradational conditions if the substrate is finer than the active layer (i.e. degradation in an armored river). Ribberink (1987) included a third layer between the active layer and the substrate to model the effects of dunes exceptionally larger than the average dune height. Although this model includes more physical mechanisms and improves the prediction of mixed sediment processes in dune-dominated cases, it may still become ill-posed ( Sieben, 1994 ).
To understand the conditions in which the active layer model becomes ill-posed we focus on how information propagates along a river. We first consider a certain reach characterized by normal flow and immobile sediment. A perturbation of the flow propagates along the river in the form of two waves traveling at speeds equal to u ± gh where u [m/s] denotes the mean flow velocity, h [m] the flow depth, and g [m/s 2 ] is the acceleration due to gravity. If sediment is mobile, yet uniform, a perturbation in bed elevation (e.g., a sediment hump) will propagate with a speed that is termed the "bed celerity" ( de Vries, 1965;Lyn and Altinakar, 2002;Stecca et al., 2014 ). As the bed elevation affects the flow, the bed elevation perturbation also induces a perturbation of the flow. Thus, under unisize sediment conditions, a perturbation of the bed elevation leads to three waves ( Fig. 1 a). Although each of the waves perturb both bed elevation and flow, two of the waves perturb mainly the flow without much change in bed level if the Froude number ( F r = u/ gh ) is sufficiently small ( de Vries, 1973;Needham, 1990;Zanré and Needham, 1994 ).
The consideration of mixed sediment (of two size fractions to simplify the example) introduces another celerity which is termed the "sorting celerity" ( Suzuki, 1976;Ribberink, 1987;Stecca et al., 2014 ). Thus, under mixed sediment conditions (with two grain sizes), a perturbation of bed elevation causes four waves. Although each wave perturbs the flow, bed elevation, and surface grain size distribution, two of these perturb mainly the flow, one mainly the bed level, and one mainly the surface grain size distribution ( Ribberink, 1987;Stecca et al., 2014 ) ( Fig. 1 b). Sieben (1994) identified a region of parameters where, for a sediment mixture consisting of two grain size classes under bed degradation into a substrate finer than the active layer, the model is unconditionally ill-posed. This occurs when the "sorting celerity" equals the "bed celerity". This was confirmed by Stecca et al. (2014) , who observed, through numerical computation of the system eigenvalues, such model behavior also in case of more than two sediment fractions.
Furthermore, Stecca et al. (2014) analytically confirmed the outcomes of Ribberink's analysis using a more realistic unsteady model for two sediment size fractions. They considered grain size selectivity of the bedload but hiding in a limited manner. Hiding accounts for the fact that grain size fractions finer than a characteristic mean grain size of the mixture hide behind larger grains and so they experience a larger critical bed shear stress compared to the unisize case ( Einstein, 1950;Komar, 1987a;1987b ). The opposite happens for coarse sediment fractions, which experience a larger exposure to the flow than in a unisize case. In their analysis Stecca et al. (2014) showed that the model can become illposed under degradational conditions if and only if the substrate is finer than a reference grain size distribution which is related to the grain size distribution of the bedload, instead of the active layer (as in Ribberink's (1987) analysis).
To overcome the problem of setting the active layer thickness, Parker et al. (20 0 0) developed a stochastic framework without the need for a distinction between the active and inactive parts of the bed. Blom and Parker (2004) , Blom et al. (2006) , and Blom et al. (2008) developed a model that accounts for dune sorting and the variability of bed elevation based on the stochastic framework developed by Parker et al. (20 0 0) . The model associates a probability of grain size selective entrainment to all elevations within the bed, and hence allows for sediment at any elevation to be entrained and contribute to the bedload discharge. Viparelli et al. (2017) developed a simplified vertically continuous model assuming slow changes in bed elevation and a steady probability distribution of entrainment, deposition, and bed elevation, which make their model suitable for large space and time domains. So far the well-posedness of the continuous model has never been assessed.
Our main objective is to analyze the problem of ill-posedness of the active layer model used for mixed sediment morphodynamics. The present paper provides four key improvements with respect to presently available knowledge: (i) we obtain analytical expressions to characterize a simplified model (i.e., to find whether it is ill-posed or well-posed) with two sediment fractions only, (ii) we study the effect of model parameter choice on ill-posedness, (iii) we find new (previously neglected) ill-posed domains, and (iv) we study the consequences of ill-posedness in numerical simulations. Our second objective is to mathematically characterize the vertically continuous model developed by Viparelli et al. (2017) . In the next section we present the general set of equations for modeling mixed sediment river morphodynamics using (a) the active layer model and (b) the vertically continuous model developed by Viparelli et al. (2017) . The models are simplified and analyzed in Section 3 . We analyze the effect of model parameters on the illposedness of the active layer model in Section 4 . In Section 5 we study the consequences of ill-posedness using numerical runs. In Section 6 we relax and study the simplifications of our analysis.

Model equations
In this section we present the equations used to model river morphodynamics. These equations represent one-dimensional hydrostatic flow over a mobile bed composed of an arbitrary number N of non-cohesive sediment fractions 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 ).
In the following section we describe the flow equations. As previous research has clarified ( Ribberink, 1987;Stecca et al., 2014 ), a key parameter in determining well-posedness of the active layer model is the active layer thickness. In this paper we both consider a model with constant, and with unsteady (time-varying) active layer thickness. While the well-posedness of the model with constant active layer thickness has been analyzed in previous work ( Ribberink, 1987;Stecca et al., 2014 ), to our knowledge no analysis of the well-posedness of the model with unsteady active layer thickness is available, although use of such a model is documented in the literature (e.g. Karim et al., 1983 ). The equations of the adapted active layer model are presented in Section 2.2 . In Section 2.3 we present the vertically continuous model derived by Viparelli et al. (2017) . The closure relations for both models are treated in Section 2.4 . In Section 2.5 we present a compact matrix formulation of the model equations.

Flow equations
The flow is described by the 1D Shallow Water Equations (i.e., the Saint-Venant equations, Saint-Venant (1871) ) considering constant width. Assuming steady flow conditions the water discharge is uniform and conservation of momentum reduces to the socalled backwater equation ( Appendix A ). When assuming steady flow over a movable bed composed of sediment of different sizes we implicitly assume that the flow adapts instantaneously to perturbations in bed elevation and grain size distribution. Worded differently, we assume that flow perturbations propagate infinitely fast relative to perturbations in bed elevation and surface grain size distribution. This assumption is referred to in literature as the "quasi-steady flow assumption" ( de Vries, 1965;Zhang and Kahawita, 1987;1990;Cao and Carling, 2002a ). The quasi-steady flow assumption is acceptable provided that the Froude number is sufficiently small, 1 − F r 2 = O(1) ( de Vries, 1973;Sieben, 1999;Lyn and Altinakar, 2002 ). Note that in this context the term "quasisteady" has a meaning different from, for instance, its use in the modeling of flood waves where "quasi-steady" refers to negligible inertia in the momentum balance.

Adapted active layer model equations
The conservation of the total amount of sediment in the bed is formulated by the Exner equation ( Exner, 1920 ). The active layer equation describes mass conservation for each size fraction ( Hirano, 1971 ). Appendix B presents the details of the active layer model.
To analyze the model with unsteady active layer thickness, we first need to set a closure relation expressing the thickness change in time. We consider an empirical empirical power relation between dune height H [m] and flow depth h [m] ( Yalin, 1964;Gill, 1971 ): Assuming that the active layer thickness L a [m] is equal to the mean dune height ( Blom, 2008 ), we relate the active layer thickness to the flow depth as follows: To obtain an equation for the active layer thickness variation we differentiate the constitutive law, Eq. (2) , with respect to time and then substitute the continuity Eq. (A.1) in it: Substitution of Eq. (3) into the active layer Eq. (B.2) yields the following adapted active layer equation:  [-] is the volume fraction content of size fraction k at the interface between the active layer and the substrate. In Fig. 2 we show a schematic representation of the main variables of the active layer model.

Simplified vertically continuous model equations
The conserved quantity in the vertically continuous model (similar to M ak in the active layer model) is the product of the cumulative probability of bed elevation ( P e [-]) and the volume fraction content of a specific grain size class k ( f k [-]) ( Parker et al., 20 0 0;Pelosi et al., 2014 ). The vertical coordinate is z [m]. To simplify the problem, the probability distribution depends on a second vertical coordinate y = z − η which is centered at the mean bed elevation.
Assuming slow changes in mean bed elevation and a constant (in time and space) probability distribution of bed elevation, Viparelli et al. (2017) obtain an equation for the change in time of the volume fraction content, f k : where p e [ m −1 ] is the probability density function of bed elevation ( Fig. 2 ). As in the active layer model, information is only advected in streamwise direction, i.e., the conservation equation does not include divergence terms in the y direction and the only independent variable in space is x . In contrast to the active layer model, there is no inactive substrate and sediment at all elevations plays a role. This is illustrated by the dependence of the probability function on the y coordinate and the gradient in the y direction in Eq. (5) . Thus, although the system of equations is one-dimensional, the mathematical character of the model is a property depending not only on the streamwise coordinate x but also on the vertical coordinate y .

Closure relations
We apply the Chézy law for the friction slope. Thus, the friction slope is proportional to the square of the mean flow velocity divided by the flow depth, S f = C f u 2 / (gh ) , where C f [-] is a nondimensional friction coefficient. For simplicity, we assume a constant nondimensional friction coefficient that is independent of the flow and bed parameters.  ( Hirano, 1971 ) and the vertically continuous model proposed by Viparelli et al. (2017) .

Table 1
Values of parameters in Eq. (6) according to several authors.
Meyer-Peter and Müller (1948) 8 1.5 0.047 Engelund and Hansen (1967) 0.05/ C f 2.5 0 Fernandez-Luque and Van Beek (1976) 5.7 1.5 0.037 -0.0455 Wong and Parker (2006) (1) 4.93 1.6 0.047 Wong and Parker (2006) (2) 3.97 1.5 0.0495 We apply a generalized form of the Meyer-Peter and Müller (1948) transport relation, in which the sediment transport rate is a power function of the excess bed shear stress: where q * bk [ −] is a nondimensional sediment transport rate [-] is the nondimensional bed shear stress of size fraction k , also known as Shields (1936) parameter, θ c [-] is the nondimensional critical bed shear stress, and ξ k [-] is the hiding coefficient. Table 1 summarizes appropriate values of A, B , and θ c according to several authors. A common hiding function is the one due to Egiazaroff (1965) ( Appendix C ). A simpler relation was developed by  : where D m [m] is a characteristic mean grain size and b is a nondimensional parameter. A value of b > 1 ( Dhamotharan et al., 1980;Misri et al., 1984;Kuhnle, 1993 ) implies that hiding is so strong that the coarser fraction(s) in the mixture is (are) more mobile than the finer one(s), i.e., reverse mobility ( Solari and Parker, 20 0 0 ). A final closure relation is required (only for the active layer model) for the volume fraction content of sediment of size fraction k at the interface between the active layer and the substrate, f I k . When the interface lowers the texture at the interface is equal to that at the topmost part of the substrate. When the interface elevation increases various relations can be applied for f I k . Hirano (1971) proposed that during aggradation the grain size distribution at the interface is equal to the one of the active layer. According to Parker (1991) also the bedload sediment plays a role in the aggradational flux to the substrate. Hoey and Ferguson (1994) combined both concepts in one parameter α [-] spanning the range [0,1] that describes the contribution of the active layer relative to the one of the bedload: where p k = q bk /q b [-] is the fraction of sediment transport rate of size fraction k .

Matrix formulation
In this section we introduce a matrix formulation to asses the well-posedness of the system of equations.
A system of partial differential equations (PDEs) can be mathematically classified as being of a hyperbolic, elliptic, or mixed type (e.g., Courant and Hilbert, 1989 ). To this end we write the problem in matrix-vector form (e.g., Toro, 2001 ): This equation is the one-dimensional quasi-linear non-conservative form of the advection equation. Q is the vector of dependent variables, A is the system matrix, and S is the vector of source terms. A system is hyperbolic at a point ( x, t ) if all the eigenvalues of matrix A are real. Physical propagation problems are modeled with hyperbolic systems of equations. If all eigenvalues are complex, the system is termed elliptic. Elliptic systems model equilibrium physical problems. If matrix A has both real and complex eigenvalues it is a mixed-type system.
A space-time dependent problem, in which we prescribe boundary conditions as a function of time and an initial condition (as is the case in modeling river morphodynamics), governed by an elliptic set of equations is ill-posed ( Hadamard, 1923;Joseph and Saut, 1990;Kabanikhin, 2008 ). This is confirmed by a perturbation analysis that shows that, if all eigenvalues of matrix A are real, perturbations of a reference state are bounded (Appendix A in the supplementary material). However, if there is at least one complex eigenvalue (or, precisely, at least two, because of the complex conjugate), perturbations grow exponentially. The exponential growth depends on the product of the imaginary part of the eigenvalues and the wave number of the perturbation, which implies that the solution of an ill-posed problem is unstable to short perturbations. Attempts to numerically integrate an ill-posed problem therefore produce results that continue to change as the grid is refined ( Woodhouse et al., 2012;Barker et al., 2015 ), as in numerical solutions perturbations always exist due to at least truncation errors. In numerical simulations the wave number of the shortest possible perturbation is inversely related to the horizontal discretization ( x ).
By only using the eigenvalues of A to characterize the system of equations we are neglecting the effect of friction (Appendix A in the supplementary material). Yet, this suffices here as friction becomes relevant for small wave numbers only (Appendix A in the supplementary material) and the most critical wave numbers, as regards to oscillation growth, are the large ones.
As a single complex pair of eigenvalues makes the problem illposed, we do not make a distinction between the number of complex eigenvalues. We term a problem with at least a pair of complex eigenvalues as elliptic.
We recast in matrix-vector form the Saint-Venant equations, (A .1) and (A .2) , the Exner equation, (B.1) , the active layer thickness equation, (3) , and the adapted active layer equation, (4) . The vector of dependent variables is: the system matrix is: (11) and the vector of source terms is: The brackets ([ ]) highlight those terms that are vectors or matrices.
We also recast in matrix-vector form the Saint-Venant equations, (A.1) and (A.2) , the Exner equation, (B.1) , and the conservation equation of the vertically continuous model, Eq. (5) . The vector of dependent variables is: the system matrix is: and the vector of source terms is:

Characterization of the mathematical models
In this section we analyze the mathematical character of the models described in Section 2 . Eigenvalues computed numerically can be obtained for an unlimited number of fractions. Here we study a simple case assuming steady flow and two size fractions to obtain analytical expressions of the eigenvalues.
As in our case the temporal change of the active layer thickness depends on the spatial gradient of the water discharge per unit width, Eq. (3) , the steady flow assumption implies a constant active layer thickness. Yet, in a numerical simulation where the steady flow assumption is used but the upstream discharge varies with time (i.e., alternating steady flow), the active layer thickness may vary with time. However, in such a case the perturbations due to a change in active layer thickness propagate infinitely fast relative to the perturbations in bed elevation and surface grain size distribution.
The implications of more than two sediment size fractions and an active layer thickness as a function of the flow depth are studied in Section 6 .

Steady active layer model consisting of two size fractions
the vector of source terms ( S alS2 ) reads: and the system matrix ( A alS2 ) is: We define ψ [-] as: which is a parameter related to the intensity of total bedload in the flow and ranges between 0 (null sediment discharge, i.e., fixed bed) and O(10 −2 ) (high sediment discharge), (e.g., de Vries, 1965;Lyn and Altinakar, 2002;Stecca et al., 2014 ).
The parameter γ 1 [-] is a measure of the fraction content of sediment in transport relative to the fraction content of sediment at the interface between the active layer and the substrate ( Stecca et al., 2014 ): where c 1 ∈ [0,1] [-] is a parameter expressing the increase in the sediment transport intensity of the fine fraction relative to the total sediment transport intensity ( Stecca et al., 2014 ): We now introduce the parameter χ 1 [-] which is a nondimensional measure of the derivative of the total sediment transport rate with respect to the volume of fine sediment in the active layer: The parameter μ 1,1 [-] is defined as: We obtain the eigenvalues of the system matrix finding the roots of its second degree characteristic polynomial. The eigenvalues are nondimensionalized dividing by the flow velocity: where the discriminant is: The eigenvalues of the system carry coupled information on both the bed elevation and the surface grain size distribution which shows that a perturbation in bed elevation causes a perturbation in surface grain size distribution and vice versa ( Section 1 ). Yet, we identify two nondimensional celerities that approximate the changes in bed elevation ( λ b ) and in surface grain size distribution ( λ s 1 ) independently.
The bed celerity, which is independent of the active layer thickness, was first derived by de Vries (1965) for unisize sediment: We define the nondimensional sorting celerity as: This sorting celerity differs from the one of Ribberink (1987) as he considered a perturbation in the mean grain size while here the sorting celerity relates to a perturbation in the volume fraction content of each grain size fraction individually. The proposed expression for the sorting celerity in Eq. (28) is a generalization of the expression proposed by Stecca et al. (2014) , as we have relaxed Stecca's assumption of limited hiding.
The mathematical character of the model depends on the sign of the discriminant Δ alS2 , Eq. (26) . If Δ alS2 > 0 the two eigenvalues are real and the system is hyperbolic. If Δ alS2 < 0 the eigenvalues are complex and the system is elliptic. A large difference between the bed celerity and sorting celerity reduces the likelihood that the model becomes elliptic. Hyperbolicity is guaranteed if γ 1 > 0. If γ 1 < 0 and the bed and sorting celerities are equal, ellipticity is guaranteed ( Sieben, 1994;Stecca et al., 2014 ).
Assuming that reverse mobility does not occur ( Section 2.4 ), c 1 is larger than the volume fraction content of fine sediment in the active layer ( F a 1 ) due to the grain size selectivity of the sediment transport relation ( Stecca et al., 2014 ). If we also assume that the sediment transferred to the substrate in aggradational conditions has the same grain size distribution as the active layer ( Hirano, 1971 ), then the parameter γ 1 is always positive in aggradational conditions. Only a substrate finer than the active layer yields a negative value of the parameter γ 1 . Thus, a two-fraction active layer model can only be ill-posed if the bed degrades into a substrate that is finer than the active layer (a result also found by Stecca et al. (2014) considering unsteady flow).
In Sections 4.1 and 4.2 we assess the relaxation of the assumptions that reverse mobility does not occur and that the aggradational flux to the substrate has the same grain size distribution as the active layer.

Steady vertically continuous model consisting of two size fractions
We apply the same procedure used to analyze the active layer model to the vertically continuous model ( Section 2.3 ). In this manner we obtain the discriminant of the eigenvalues ( Appendix D ): where λ sc 1 , g 1 , and m 1,1 are the equivalents to λ s 1 , γ 1 , and μ 1,1 of the active layer model ( Appendix D ). Similar to the active layer model ( Section 3.1 ), the continuous model is hyperbolic and well-posed if Δ vcS2 > 0 and vice versa. Although the expression of the discriminant of the vertically continuous model, Eq. (29) , is similar to the one of the active layer model, Eq. (26) , there is an essential difference between the two. In the active layer model the discriminant is a function of the streamwise position, Δ alS2 ( x ), yet in the continuous model the discriminant is also a function of the vertical coordinate, Δ vcS2 ( x, y ). Thus, ellipticity or hyperbolicity is a property not only of the streamwise coordinate but also of the elevation in the bed ( Section 2.3 ). Hyperbolicity is guaranteed if g 1 > 0 but, contrary to the active layer model, this parameter can be negative both under aggradational and degradational conditions. Due to grain size selective transport we can assure that, if reverse mobility conditions do not prevail, the concentration c 1 is larger than the volume fraction content representative of the bed surface F b 1 , Eq. (C.4) . However, F b 1 is a weighted average of all sediment and for this reason there is no guarantee that for all bed elevations the average volume fraction content ( F b 1 ) is larger than the local volume fraction content in the bed sediment ( f 1 ). Moreover, as there is no distinction between aggradational and degradational cases, the domain in which the model is likely to be ill-posed is larger than for the active layer model. The presence of fine sediment at the locations having larger probability of entrainment in combination with a "smooth" vertical variation (small derivative) of the volume fraction content of fine sediment reduces the likelihood of the model becoming elliptic.

Active layer model parameter study
In this section we assess the effects of various model parameters on the mathematical character of the active layer model. To this end we study the analytical expressions of the eigenvalues of the steady model considering two sediment size fractions obtained in Section 3.1 .

Hiding
Given the fact that ill-posedness arises when considering different grain sizes in the mixture and that a larger difference between grain sizes increases the ill-posed domain, intuitively hiding should reduce the likelihood of ill-posedness. Its effect, however, is opposite as we will show here.
To explain this counter-intuitive result we analyze the term in the characteristic polynomial intrinsically related to hiding. This term is the derivative of the sediment transport rate of fine and coarse sediment with respect to the volume of fine sediment in the active layer ( ∂ q bk / ∂ M a 1 ). It can be considered as the summation of two terms: We name the first and second terms on the right-hand side the "presence term" and the "hiding term", respectively. The "presence term" explains that an increase in the volume fraction content of the fine sediment in the active layer implies both (1) an increase of the sediment transport rate of the fine fraction as its presence at the bed surface is larger, and (2) a consequent decrease of the sediment transport rate of the coarse fraction because its presence at the bed surface decreases. The "hiding term" indicates the fact that a variation of the volume fraction content of the fine sediment changes the sediment transport rate of both fine and coarse fractions due to a change in the mean grain size of the sediment mixture. The "presence term" is positive for the fine fraction and negative for the coarse fraction. The "hiding term" is always positive.
In a situation where hiding is negligible, an increase of the characteristic size of the coarse fraction or decrease of the fine fraction, which is associated with a larger likelihood that the model is elliptic, causes an increase of the "presence term" and thus of ∂ q bk / ∂ M a 1 . With respect to such a situation, hiding decreases the "presence term" of both fine and coarse sediment (reducing the likelihood of ellipticity) but introduces the positive contribution of the "hiding term". Overall, the "hiding term" may dominate, which increases the value of ∂ q bk / ∂ M a 1 and thus of the likelihood of ellipticity.
Interestingly, in degradational conditions into a fine substrate (a situation prone to be elliptic) the hiding term dominates. Thus, hiding increases the likelihood that the model is elliptic in degradational conditions into a substrate finer than the active layer.
In Fig. 3 a we show the effect of hiding on the discriminant of the steady active layer model considering two size fractions, Eq.  26) . We consider the reference case described in Table 2 . The sediment transport rate is computed using the relation derived by Meyer-Peter and Müller (1948) . To obtain different values of hiding we vary parameter b in the power law hiding function in Eq.
(7) between 0 and 1 (purple line in Fig. 3 a). The yellow line in Fig. 3 a is obtained varying the characteristic grain size of the fine fraction between 0.001 m and 0.004 m using the Egiazaroff hiding relation, Eq. (C.5) . The discriminant decreases for increasing hiding independent from the hiding function. Besides these cases under degradational conditions, we may encounter problems even under aggradation. In fact, if hiding is so strong that reverse mobility is induced, then one of the assumptions of the analysis by Stecca et al. (2014) , may not be fulfilled. In detail, it may happen that the reference content c 1 related to the fine sediment in the bedload, Eq. (21) , is not greater than the content of fines in the active layer F a 1 , which was their assumption under grain size selective transport. When reverse mobility instead determines conditions such that c 1 < F a 1 , then the discriminant, Eq. (26) , may be negative, and the model may become elliptic even under aggradational conditions.

Aggradational flux to the substrate
The sediment transferred to the substrate under aggradational conditions using the model by Hoey and Ferguson (1994) ( Section 2.4 ) is always finer than the sediment in the active layer. This is because the bedload is finer than the bed surface due to grain size selective processes (provided that reverse mobility does not dominate). Thus, application of the model by Hoey and Ferguson (1994) implies that under aggradational conditions the interface between the active layer and the substrate is finer than the active layer. This means that the condition γ 1 > 0 ( Section 3.1 ) may not be fulfilled under aggradational conditions, which implies that the model may become ill-posed. Therefore, a larger contribution of the bedload to the aggradational flux to the substrate (smaller value of the parameter α in Eq. (8) ) implies a larger likelihood of the model becoming elliptic (Appendix B in the supplementary material).
However, in a hypothetical aggrading case in which the grain size distribution transferred to the substrate is fully composed of bedload sediment ( α = 0 ), the relative content of the fine fraction in the vertical sediment flux, γ 1 ( Eq. (20) ), that controls the size of the ill-posed domain ( Section 3.1 ), is still not as small as it can be found under degradational cases (Appendix B in the supplementary material). Thus, ill-posed cases are expected to occur primarily under degradational conditions into a fine substrate.

Prefactor in a sediment transport relation and morphodynamic factor
The discriminant ( Δ alS2 ) of the steady active layer model for two size fractions, Eq. (26) , can be written as Δ alS2 = A 2 Δ alS2 , where alS2 is the discriminant for a unit prefactor (i.e., A = 1 ). The prefactor A increases or decreases the discriminant but does not change its sign, and so it does not change the character of the mathematical system. This is confirmed by Fig. 3 b, which shows the effect of varying the prefactor A in the reference case described in Section 4.1 . Since morphodynamic time scales are usually several orders of magnitude larger than the time scales of the flow ( Section 2.1 ), computations usually cover a significant number of years. The computational time is sometimes reduced using a morphodynamic factor that multiplies the divergence of the sediment transport rate ( Latteux, 1995;Roelvink, 2006;Ranasinghe et al., 2011 ). This factor can also be considered as a multiplication of the sediment transport rate and therefore has the same effect as the prefactor A . Thus, the use of a morphodynamic factor does not change the mathematical character of the model. This result is obtained assuming quasi-steady flow. While the prefactor in the sediment transport relation rarely varies by more than an order of magnitude, simulations may be run with morphodynamic acceleration factors O(10 2 ) . In these latter cases, the quasi-steady flow assumption may not be acceptable, which limits the extension of our analysis.

Exponent and critical Shields stress in a sediment transport relation
The discriminant, Eq. (26) , tends to 0 − with increasing values of B if the effective Shields stress for all sediment fractions is smaller than 1, or to ∞ if the effective Shields stress is larger than 1 for at least one fraction. Thus, it is difficult to generalize the effect of the exponent. Its variation from a reference situation can both make the system hyperbolic if the reference situation is elliptic or vice versa. In Fig. 3 c we show the discriminant as a function of B for the same reference cases as in Section 4.3 . The hyperbolic situation when using Meyer-Peter and Müller (1948) becomes elliptic if the value of the exponent B increases towards the value in Engelund and Hansen (1967) .
The effect of the critical Shields stress, θ c in Eq. (6) , on the discriminant is similar to the effect of the exponent, as its variation can both make a previously hyperbolic case elliptic or vice versa. Fig. 3 d shows how a decrease of the critical shear stress when using the sediment transport relation by Meyer-Peter and Müller (1948) increases the discriminant reducing the likelihood of elliptic behavior.

Active layer thickness
The discriminant of the eigenvalues, Eq. (26) , can be written as a second degree polynomial of the inverse of the active layer thickness, i.e., Δ alS2 = a 1 (1 /L a ) 2 + a 2 (1 /L a ) + a 3 where a 1 > 0, a 2 , and a 3 are coefficients independent of the active layer thickness. This implies that: (1) the model is well-posed for a sufficiently thin active layer, (2) the model is well-posed for a sufficiently thick active layer, and (3) there exists one ill-posed domain only (regarding the active layer thickness). These results of the two-fractions model confirm previous results based on the simplified active layer model ( Ribberink, 1987;Sieben, 1994 ).
The inverse of the roots of the second degree polynomial are the limit values of the active layer thickness that ensure that the model is well-posed: where we have used the notation ( ) for the variables with unit active layer thickness (i.e., L a = 1 m). Given the facts that the active layer thickness is one of the most empirical parameters of the system of equations ( Section 1 ) and that river morphodynamic models often require calibration (e.g., Cao and Carling, 2002b ), Eq. (31) can be applied to select a certain value for the active layer thickness to avoid a situation that is prone to be ill-posed.

Consequences of ill-posedness
In this section we analyze the consequences of ill-posedness using numerical simulations. Our aim is to provide modellers with the tools to detect occurrence of ill-posedness in their results and understand how the observed unrealistic model behavior changes First we make numerical runs to qualitatively observe the consequences of the non-linear effects that are neglected in the perturbation analysis ( Section 5.1 ). In Section 5.2 we conduct a sensitivity analysis to generalize the consequences observed in the previous section.

Numerical examples
The linear perturbation analysis shown in Section 2.5 indicates that perturbations grow unboundedly if the model is elliptic. In this section we run four numerical simulations at flume scale to analyze the effects of the neglected non-linear terms. The simulations are one-dimensional for the flow and are computed using the Delft3D software package ( Lesser et al., 2004 ), which solves the unsteady Shallow Water Equations in combination with the active layer model. For simplicity the active layer thickness is assumed constant. Under aggradational conditions the sediment transferred to the substrate is composed of only the active layer sediment (i.e., α = 1 in Eq. (8) ). Substrate stratigraphy is stored using a bookkeeping system ( Viparelli et al., 2010;Stecca et al., 2016 ). All simulations start from equilibrium conditions under coarse sediment feeding. A lowering of the base level is imposed, which causes degradation into a fine substrate. We consider a well-posed reference case (Simulation 1) that (initially) has the same parameters as the reference case of the parameter study of Section 4 ( Table 2 ). Then, the active layer thickness is changed (Simulation 2) and hiding is considered (Simulation 3). Simulation 4 is equal to Simulation 3 except for its morphodynamic factor. Table 3 summarizes the differences between the four simulations. The boundary conditions that are in equilibrium with the initial condition  as well as other parameters are described in Table 4 .
In Fig. 4 a we plot the discriminant of the eigenvalues of the quasi-steady active layer model, Eq. (26) , at the initial time as a function of the active layer thickness. Note that the conditions of Simulation 1 yield a well-posed model ( Δ alS2 > 0) while the conditions of Simulation 2 yield an ill-posed model ( Δ alS2 < 0). In Fig. 4 d-e we show the evolution of the bed elevation for Simulations 1 and 2, respectively. The ill-posed Simulation 2 shows an oscillatory behavior that is not present in the well-posed Simulation 1. We have not imposed any initial perturbation, which implies that numerical noise is sufficient to trigger the oscillatory behavior. Simulation 3 is the same as Simulation 1, yet the sediment transport rate now accounts for hiding using the  function, Eq. (7) , with exponent b equal to 0.8. Simulation 3 is ill-posed ( Fig. 4 b) and the solution shows oscillations ( Fig. 4 f) just as in the ill-posed Simulation 2. Simulation 4 is the same as Simulation 3 except for its morphodynamic factor equal to two, which decreases the value of the discriminant ( Fig. 4 c), causing oscillations to develop faster ( Fig. 4 g).
For all cases oscillations do not occur at the upstream end of the domain. This is because oscillations require time (and so space) to grow. In all cases oscillations grow until a maximum amplitude is reached and then propagate downstream. This maximum amplitude is such that the conditions are at the brink of ill-posedness and well-posedness. Worded differently, downstream from the location where the amplitude is maximum the model is ill-posed and it is well-posed upstream from it.
The oscillations are associated with degradation and subsequent aggradation. The deposited sediment has the same grain size distribution as the active layer (which is coarser than the initial substrate), so the overall effect of an oscillation is a coarsening of the topmost part of the substrate. This coarsening acts as a regularization mechanism, which not only restores hyperbolicity but also dampens oscillations that arrive from upstream by limiting the source of fine material.
It is likely that, because of the regularization mechanism, computations do not crash. As a result the extent and likelihood of ellipticity may in practice be underestimated. Yet, the results are physically unrealistic and implementing an automated check of the eigenvalues would be good practice for software developers and users.

Sensitivity analysis
The previous section has shown that, due to the non-linearity of the system, an ill-posed simulation generates non-physical oscillations that propagate downstream and grow until a certain maximum amplitude at which the mathematical problem is at the brink of ill-posedness and well-posedness. In this section we run a sensitivity analysis to generalize those results.
To this end, we vary 4 physical parameters and 2 parameters related to the domain discretization, using Simulation 1 as a reference case. Table 5 summarizes the parameter values used in the sensitivity analysis.
As we are interested in studying the behavior of simulations under ill-posed conditions we exclude from the analysis those simulations in which the combination of parameters yield a wellposed model. As we have observed that the oscillations need space to grow until a maximum value ( Section 5.1 ), we exclude those simulations in which the domain is not long enough to develop an oscillation that travels with a constant amplitude. A set of 173 out of 256 simulations fulfills these two requirements. Table 4 Domain definition, boundary conditions, and numerical parameters. The symbols not defined in the text are: reach length ( L ), channel width ( B ), simulation time ( T ), lowering rate of the downstream water level (low. rate), horizontal discretization length ( x ), and time step ( t ).

L [m]
B   Table 3 for the parameters definition.
In Fig. 5 a we plot the maximum flow depth ( h max ) nondimensionalized with the normal flow depth ( h n ), as a function of the discriminant, Eq. (26) . For the sake of clarity we plot the results of the simulations with a horizontal discretization length ( x ) equal to 0.1 m and a thickness of the substrate layers equal to 0.01 m and 0.10 m (see Appendix C of the supplementary material for the results of all simulations). The parameters used to evaluate the discriminant are those at the start of the simulation (normal flow). The vertical black lines join simulations with the same physical parameters (i.e., they only differ regarding numerical parameters) and the color of each dot is related to the thickness of the substrate layer. The linear analysis has shown that the growth rate depends on the discriminant ( Section 2.5 ), however there is only mild correlation between the discriminant and the maximum amplitude of the oscillations.
A thinly discretized substrate is associated with a larger amplitude of the oscillations ( Fig. 5 a). This effect can be seen only empirically since it is not a parameter of the system of equations nor does it appear in the linear stability analysis. For all simulations we compute the flow depth that yields a value of the discriminant at the initial condition equal to 0 (i.e., at the brink between ellipticity and hyperbolicity). This is done numerically finding the root of Δ alS2 ( h ), Eq. (26) , considering the water discharge and volume of sediment in the active layer and at the substrate of the initial condition. We term this flow depth the hyperbolic flow depth ( h hyp ), which is independent of the numerical parameters of the simulation and depends on physical parameters only. In Fig. 5 b we com-pare the measured maximum flow depth and the hyperbolic flow depth. The grey line represents the situation in which h max = h hyp . We see that the hyperbolic flow depth can be used as a rough estimate of the maximum flow depth that will occur in an elliptic simulation. One important source of scatter is the fact that the hyperbolic flow depth depends on the initial condition only, whereas the maximum flow depth also depends on the evolution of the solution as the oscillations interact with each other.
In Fig. 5 c we plot the nondimensional maximum flow depth as a function of the streamwise location where the maximum flow depth occurs (nondimensionalized with the total length of the domain). For the sake of clarity we plot the results of the simulations with a thickness of the bookkeeping layers ( z ) equal to 0.01 m (see Appendix C of the supplementary material for the results of all simulations). In the thinly discretized simulations we find the maximum amplitude more upstream compared to the coarsely discretized simulations. The location where the maximum amplitude of the oscillations is found is related to its growth rate since a faster growing oscillation develops its maximum amplitude in less distance than a slower one. This result confirms the findings of linear stability analysis ( Section 2.5 ).
Also the maximum flow depth and the domain discretization are mildly correlated. Similar to the vertical discretization, smaller cells yield a larger maximum flow depth but this effect can be seen only empirically.
The discretization of the substrate affects also the duration of the elliptic behavior. Fig. 6 shows the longitudinal profile of four  simulations from the sensitivity analysis at the end of the run (Simulations 1, 5, 6, and 7, see Table 3 ). The substrate of the wellposed simulations at the end of the runs is unaltered, whereas it has coarsened in the ill-posed cases due to the oscillatory behavior, which acts as a regularization mechanism ( Section 5.1 ). A thinly discretized substrate enhances the regularization mechanism as one full layer is created with the grain size distribution of the (coarse) active layer during the aggrading phase of the oscillation. If the bed is discretized into thick layers the material transferred to the substrate will be averaged with the sediment already present in the top substrate layer. The resulting grain size distribution of the substrate may not be sufficiently coarse to prevent the model from being ill-posed.

Implications of considering more than two size fractions or an unsteady active layer thickness
To obtain an analytical expression of the eigenvalues we have restricted our analysis to mixtures of sediment composed of two size fractions and steady flow. In this section we explore the consequences of relaxing these assumptions.

Ill-posed domain of a three-size-fractions case
A model for three sediment size fractions is too complex to obtain analytical expressions of the eigenvalues. We therefore first attempt to provide insight addressing a specific case based on the results of the case for two size fractions.
The concept of a finer or coarser active layer relative to the substrate is unequivocally applicable in the case of two size fractions. However, this concept is not as straightforward for three size fractions, as it requires the definition of a mean grain size. As an extension of the results for the two size fractions case where the model can only be ill-posed in degradational conditions into a substrate finer than the active layer (assuming certain conditions on the closure relations, Section 3.1 ), we consider a situation with three size fractions which, regardless of the method to compute the mean grain size, is governed by degradation into a substrate coarser than the active layer. This happens, for instance, if the volume fraction contents in the active layer of the fine, medium, and coarse size fractions are 0.5, 0.5, and 0, respectively; and at the interface are 0.5, 0, and 0.5.
We consider a sediment mixture with the above volume fraction contents and characteristic grain sizes of the fine, medium,   ( Table 2 ). This situation is elliptic as two of the eigenvalues of the system matrix are complex. Thus, in a three size fractions case the mean grain size of the sediment in the active layer relative to that at the interface is not a valid discriminant of the mathematical character of the system of equations. Fig. 7 shows the results of a numerical simulation based on the above parameters. The solution presents oscillations as in the previous ill-posed cases. However, the amplitude of these oscillations is now significantly smaller (compare Fig. 7 a to Fig. 4 e). A relatively large oscillation appears after approximately 3 h which entrains coarse sediment from the substrate ( Fig. 7 b). During the aggradational phase of the oscillations fine sediment from the active layer is transferred to the substrate. Thus, at the end of the simulation the top part of the substrate is finer than initially ( Fig. 7 c). To illustrate the implications of this result we study the effects of discretizing the same sediment mixture into two or three sediment fractions. To discretize the sample into three sediment fractions we use characteristic grain sizes equal to 0.0 01 m,0.0 03 m and 0.0 05 m and to discretize it into two sediment fractions we use 0.002 m and 0.004 m. The volume fraction content in the medium size of the three fraction mixture is equally split between the fine and coarse bins of the two fraction mixture. We vary all volume fraction contents between 0 and 1 to obtain different sediment mixtures and the flow depth (keeping the water discharge per unit width constant) between 0.15 m and 1.5 m to obtain different flow conditions. All other parameters are equal to the reference case.  Fr ) and the difference between the mean grain size of the sediment in the active layer ( D ma ) and at the interface between the active layer and the substrate ( D mI ). When the mixture is discretized into three size fractions the model may be ill-posed under degradational conditions into a substrate coarser than the active layer. Fig. 8 a-b shows the elliptic domain when the mixture is discretized into two and three size fractions, respectively. On the vertical axis we plot the difference between the mean grain size of the sediment in the active layer ( D ma [m]) and at the interface between the active layer and the substrate ( D mI [m]). Note that some situations that are well-posed when the mixture is discretized into two fractions are ill-posed when it is discretized into three fractions. We cannot prove that the eigenvalues of the system matrix for three size fractions are always real under aggradational conditions due to the complexity of the expressions. Nevertheless, we have not obtained a single complex value in any of the aggradational tests we have conducted.

Effect of an unsteady active layer thickness in the ill-posed domain
In this section we analyze the implications of considering a variable active layer thickness with respect to the ill-posedness of the system of equations. The flow needs to be considered unsteady to study the variability of the active layer thickness if it is related to dune growth ( Section 3 ). We simplify the system assuming two sediment size fractions and negligible hiding ( Stecca et al., 2014 and Appendix D of the supplementary material).
We obtain the characteristic polynomial of the system matrix, Eq. (11) . We prove that in aggradational and degradational conditions into a substrate coarser than the active layer the characteristic polynomial has 5 real roots (Appendix D of the supplementary material). Therefore the model is well-posed regardless of the unsteady active layer thickness. Regarding degradational conditions into a substrate finer than the active layer we prove that if λ b > λ s 1 F a 1 / f I 1 , considering a variable active layer thickness increases the likelihood of the model becoming elliptic. Note that, assuming a similar order of magnitude of the bed and sorting celerities, in conditions prone to be elliptic (i.e., degradation into a substrate significantly finer than the active layer) a variable active layer thickness increases the domain in which the active layer model is elliptic. We numerically test several sets of parameters and we find no case where the model is hyperbolic if the active layer is unsteady but elliptic if it is constant (Appendix D of the supplementary material). This suggests that, although we do not provide a formal proof, an unsteady active layer thickness always increases the likelihood of the model being ill-posed.

Conclusions
We have assessed the well-posedness of the equations used to model mixed sediment river morphodynamics. In particular we have studied the system formed by the flow equations ( Saint-Venant, 1871 ) together with the active layer model ( Hirano, 1971 ) and a simplified vertically continuous model . Our findings are the following: • Considering two size fractions and the quasi-steady flow assumption we obtain an analytical expression for the discriminant that determines whether the active layer and continuous models are ill-posed. • Assuming (i) two size fractions, (ii) steady flow, (iii) no reverse mobility, and (iv) that under aggradational conditions the depositional flux of sediment to the substrate is entirely composed of active layer sediment, the active layer model can be ill-posed under degradational conditions into a substrate finer than the active layer only. • The use of a hiding factor increases the likelihood of illposedness. Strong hiding that causes reverse mobility may cause ill-posedness also in aggradational conditions. • Aggradational cases may be ill-posed if the depositional flux of sediment to the substrate includes bedload sediment ( Hoey and Ferguson, 1994 ). • The active layer model may be ill-posed in degradational conditions into a substrate coarser than the active layer if more than two size fractions are considered. • Considering a variable active layer thickness associated with dune growth increases the likelihood that the active layer model is ill-posed. • The simplified vertically continuous model can be ill-posed under both aggradational and degradational conditions. A small vertical gradient of the probability of bed elevation and volume fraction content decreases the likelihood of the model being illposed. • Ill-posedness results in non-physical oscillations that grow until a maximum amplitude is reached, at which the model recovers its hyperbolic character (and becomes well-posed). The non-physical oscillations itself act as a regularizing mechanism by coarsening the substrate.
The numerical solution of an ill-posed problem may be reliable if perturbations do not have space and/or time to grow or if the consequences of the perturbations are negligible compared to the accuracy of the problem data. However, the reliability of the solution becomes subjective. This implies that it is up to the modeller to decide whether a solution is representative of the physical phenomenon under consideration.
In a well-posed model a finer grid provides more accurate results. This is opposite in ill-posed models as the growth rate of oscillations decreases with grid size. Thus, if a model is ill-posed, one may be tempted to use a larger grid size such that oscillations do not have space to grow and numerical viscosity is sufficient to suppress the consequences of ill-posedness. We do not recommend to follow this strategy because of the subjectivity of the solution.
We do not recommend discarding the active layer and vertically continuous models for modeling mixed sediment river morphodynamics. The former has proven its validity over a large range of situations ( Section 1 ) and the latter is yet a simplified version of a continuous sediment conservation model. Moreover, both models are well-posed for a vast range of situations.
The ill-posedness of the system of equations is a fundamental mathematical problem independent of the numerical solver. It can only be solved by an improved set of equations that represents physical processes in a better way than existing models do. In this regard we are currently conducting laboratory experiments to investigate the physical mechanisms that are relevant under conditions in which the active layer model is ill-posed. One may want to introduce minimal changes to the active layer model to regularize it (i.e., to ensure that the model is always well-posed). In this case, the most straightforward solution is to check whether the model is ill-posed and to change the active layer thickness to a value that provides a well-posed model. One needs to be aware that this approach is not simply a numerical trick, as it implies a change of the time scale of the physical processes under consideration. Moreover, it implies a temporal change of the active layer thickness which may be relatively large and local. Preliminary simulations show that this solution is not always stable. Another possibility may be to artificially modify the celerities without changing the actual thickness of the active layer. A similar approach has been used by Zanotti et al. (2007) to regularize the ill-posed twolayer shallow water model. Current work by the authors builds on this idea.

Acknowledgments
This research is part of the research programme RiverCare, supported by the Applied and Engineering Sciences (TTW) division of the Netherlands Organization for Scientific Research (NWO), and which is partly funded by the Ministry of Economic Affairs under grant number P12-14 (Perspective Programme). Guglielmo Stecca is supported by a Marie Curie Outgoing Fellowship within the European Union's Seventh Framework Programme for research, technological development and demonstration under grant agreement no. PIOF-GA-2013-621886. We acknowledge the editorial work by Dr. D'Odorico and the comments of two anonymous reviewers, who helped in clarifying the message. The fruitful discussions and comments on the manuscript of Liselot Arkesteijn and Robert Jan Labeur are gratefully acknowledged.

Appendix A. Flow equations
The water phase is mathematically described by the Saint-Venant equations ( Saint-Venant, 1871 )  Considering steady flow, the conservation of water mass, Eq. (A.1) , reduces to a spatially constant discharge, and the conservation of momentum, Eq. (A.2) , to the backwater equation:

Appendix B. Active layer equations
The conservation of the total amount of sediment in the bed is represented by the Exner equation ( Exner, 1920 ): where q b [m 2 /s] is the sediment transport rate per unit width mul- is the bed porosity (i.e., the sediment transport rate q b accounts for pores). For simplicity, mechanisms such as subsidence and uplift, compaction and dilation of sediment are neglected in the above equation ( Paola and Voller, 2005 ). Of special relevance is the implicit assumption that the temporal change of the storage of sediment within the water column and its effects on bed elevation are negligible ( Park and Jain, 1987;Stevens, 1988;Correia et al., 1992;Morris and Williams, 1996 ). We consider that there is no lag between changes in bottom bed shear stress ( Bell and Sutherland, 1983;Jain, 1992 ). Worded differently, the sediment transport rate is at capacity and adapts instantaneously to the flow field. Thus, the sediment transport rate does not require a constitutive equation and is treated as a closure relation.
Assuming constant porosity and density, the active layer equation describes the conservation of mass of grain size fraction k in the active layer ( Hirano, 1971 ): where M ak [m] is the volume of sediment of size fraction k in the active layer per unit of bed area, f I k [-] is the volume fraction content of size fraction k at the interface between the active layer and the substrate ( f I k ∈ [0 , 1] ), L a [m] is the active layer thickness, and q bk [m 2 /s] is the sediment transport rate per unit width of size fraction k multiplied by 1 / (1 − p) . The addition of the sediment transport rate for each size fraction equals the total amount of sediment in transport including pores: Assuming constant porosity and density, mass conservation of sediment of size fraction k in the substrate is expressed by: The summation of N active layer equations yields the Exner equation ( Ribberink, 1987;Parker et al., 20 0 0 ), as the active layer equation, (B.2) , represents fractional mass conservation of sediment and the Exner equation, (B.1) , represent the conservation of the total amount of sediment. Thus, to consider N active layer equations is equivalent to considering N − 1 active layer equations and the Exner equation. We here choose for the second option, as in this way the conservation of sediment mass per size fraction can be considered as an extension of the unisize model.
The substrate equation, (B.4) , is a linear combination of the Exner equation, (B.1) , and the active layer equation, (B.2) , which means that the substrate equation does not play a role in the mathematical behavior of the system and can be treated in a decoupled manner.

Appendix C. Sediment transport closure relation
The sediment transport rate of size fraction k per unit width (including pores), q bk , is expressed as the product of the volume fraction content of size fraction k at the bed surface ( F bk [-]) and the sediment transport capacity Q bk [m 2 /s] which is the sediment transport we would obtain if the bed was formed by unisize sediment yet including hiding effects ( Deigaard and Fredsøe, 1978;Ribberink, 1987;Armanini, 1995 ): (C.1) The sediment transport capacity is the product of a nondimensional sediment transport rate ( q * bk [ -] ) and the parameter gRd 3 k ( Einstein, 1950 ): where we account for the volume of pores multiplying by 1 − p. R = ρ s /ρ w − 1 [-] is the submerged specific gravity, ρ s = 2650 kg/m 3 the sediment density, and ρ w = 10 0 0 kg/m 3 the water density. The volume fraction content of size fraction k at the bed surface is constrained by the condition: In the active layer model the volume fraction content of size fraction k at the bed surface is considered to be equal to the volume fraction content in the active layer ( F bk = F ak ). In the vertically continuous model developed by Viparelli et al. (2017) , it is considered equal to the integral of the volume fraction content of size fraction k in the bed sediment weighted by the elevation's exposure to the flow: F bk = + ∞ −∞ f k (y ) p e (y ) d y.

(C.4)
The sediment transport rate q * bk is related to the mean characteristics of the flow. Here we consider a generalized form of the Meyer-Peter and Müller (1948) transport relation, which estimates sediment transport as a power function of the excess bed shear ( Eq. 6 ). The nondimensional bed shear stress of size fraction k or Shields (1936) parameter is computed as θ k = C f u 2 / (gRd k ) [-].
A commonly used hiding relation is the one due to Egiazaroff (1965) : where D m [m] is a characteristic mean grain size of the mixture obtained as an average of the grains in movement and in the bed surface. In practical terms, D m is computed as the arithmetic mean (e.g, Wu et al., 20 0 0 ), geometric mean (e.g., Bettess and Frangipane, 2003 ) or the median grain size (e.g., Kleinhans et al., 2002 ) of the bed surface sediment. A simpler expression was developed by  : where the characteristic mean grain size is the median grain size ( D 50 ) of the subpavement sediment (below an armor layer)  or the geometric mean of the surface sediment ( Parker, 1990 ). If the nondimensional parameter b is equal to 0, there is no hiding effect and each grain size behaves independently of each other. If b = 1 , the sediment transport of each size fraction is independent of its grain size (for B = 1 . 5 ), thus only depends on its presence at the surface ( F bk ). Buffington and Montgomery (1997) made an inventory of values of b spanning between 0.32 and 1.25. A value of b > 1 implies reverse mobility ( Solari and Parker, 20 0 0 ).
In this paper we compute the characteristic mean grain size, D m , as the geometric mean: where d ref = 1 mm is a reference grain size that makes the grain size on φ-scale nondimensional.

Appendix D. System of equations of the steady vertically continuous model consisting of two size fractions
The vector of dependent variables ( Q vcS2 ) is:

Supplementary material
Supplementary material associated with this article can be found, in the online version, at 10.1016/j.advwatres.2018.02.011