2.1 Primary flow equations

The primary flow is described using the depth-averaged shallow water equations (e.g. Vreugdenhil Reference Vreugdenhil1994):

(2.1) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}q_{x}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}q_{y}}{\unicode[STIX]{x2202}y}=0, & & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}q_{x}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(q_{x}^{2}/h+gh^{2}/2)}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}\left(\displaystyle \frac{q_{x}q_{y}}{h}\right)}{\unicode[STIX]{x2202}y}+gh\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}x}-F_{sx}\nonumber\\ \displaystyle & & \displaystyle \quad =2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\unicode[STIX]{x1D708}h\frac{\unicode[STIX]{x2202}\left(\displaystyle \frac{q_{x}}{h}\right)}{\unicode[STIX]{x2202}x}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\unicode[STIX]{x1D708}h\left(\displaystyle \frac{\unicode[STIX]{x2202}\left(\displaystyle \frac{q_{x}}{h}\right)}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}\left(\displaystyle \frac{q_{y}}{h}\right)}{\unicode[STIX]{x2202}x}\right)\right)-ghS_{fx},\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}q_{y}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(q_{y}^{2}/h+gh^{2}/2)}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}\left(\displaystyle \frac{q_{x}q_{y}}{h}\right)}{\unicode[STIX]{x2202}x}+gh\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}y}-F_{sy}\nonumber\\ \displaystyle & & \displaystyle \quad =2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\unicode[STIX]{x1D708}h\frac{\unicode[STIX]{x2202}\left(\displaystyle \frac{q_{y}}{h}\right)}{\unicode[STIX]{x2202}y}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\unicode[STIX]{x1D708}h\left(\frac{\unicode[STIX]{x2202}\left(\displaystyle \frac{q_{y}}{h}\right)}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}\left(\displaystyle \frac{q_{x}}{h}\right)}{\unicode[STIX]{x2202}y}\right)\right)-ghS_{fy},\end{eqnarray}$$

where $(x,y)$ (m) are Cartesian coordinates and $t$ (s) is the time coordinate. The variables $(q_{x},q_{y})=(uh,vh)$ ( $\text{m}^{2}~\text{s}^{-1}$ ) are the specific water discharges in the $x$ and $y$ direction, respectively, where $h$ (m) is the flow depth and $u$ ( $\text{m}~\text{s}^{-1}$ ) and $v$ ( $\text{m}~\text{s}^{-1}$ ) are the depth-averaged flow velocities. The variable $\unicode[STIX]{x1D702}$ (m) is the bed elevation and $g$ ( $\text{m}~\text{s}^{-2}$ ) the acceleration due to gravity. The friction slopes are $(S_{fx},S_{fy})$ ( $-$ ) and the diffusion coefficient $\unicode[STIX]{x1D708}$ ( $\text{m}^{2}~\text{s}^{-1}$ ) is the horizontal eddy viscosity. The depth-averaging procedure of the equations of motion introduces terms that originate from the difference between the actual velocity at a certain elevation in the water column and the depth-averaged velocity. We separate the contributions due to turbulent motion and secondary flow caused by the flow curvature. The contribution due to turbulent motion is accounted for by the diffusion coefficient. Elder (Reference Elder1959) derived an expression for the diffusion coefficient that accounts for the effect of turbulent motion on the depth-averaged flow assuming a logarithmic profile for the primary flow and negligible effect of molecular viscosity:

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D708}_{E}={\textstyle \frac{1}{6}}\unicode[STIX]{x1D705}hu^{\ast }, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}=0.41$ ( $-$ ) is the von Kármán constant and $u^{\ast }=\sqrt{C_{f}}Q/h$ ( $\text{m}~\text{s}^{-1}$ ) is the friction velocity. Parameter $C_{f}$ ( $-$ ) is a non-dimensional friction coefficient, which we assume to be constant (Ikeda, Parker & Sawai Reference Ikeda, Parker and Sawai1981; Schielen, Doelman & De Swart Reference Schielen, Doelman and de Swart1993) and $Q=\sqrt{q_{x}^{2}+q_{y}^{2}}$ ( $\text{m}^{2}~\text{s}^{-1}$ ) is the module of the specific water discharge. In the numerical simulations we will assume the eddy viscosity to be a constant equal to the value given by $\unicode[STIX]{x1D708}_{E}$ in a reference state (e.g. Falconer Reference Falconer1980; Lien et al. Reference Lien, Hsieh, Yang and Yeh1999). Appendix A presents the limitations of the coefficient derived by Elder (Reference Elder1959).

The terms ( $F_{sx},F_{sy}$ ) ( $\text{m}^{2}~\text{s}^{-2}$ ) account for the effect of secondary flow. These terms are responsible for a transfer of momentum that shifts the maximum velocity to the outer bend (Kalkwijk & De Vriend Reference Kalkwijk and de Vriend1980), as well as for a sink of energy in the secondary circulation (Flokstra Reference Flokstra1977; Begnudelli, Valiani & Sanders Reference Begnudelli, Valiani and Sanders2010). We deal with these terms in § 2.2.

We assume a Chézy-type friction:

(2.5a,b ) $$\begin{eqnarray}\displaystyle S_{fx}=\frac{C_{f}q_{x}Q}{gh^{3}},\quad S_{fy}=\frac{C_{f}q_{y}Q}{gh^{3}}. & & \displaystyle\end{eqnarray}$$

One underlying assumption of the system of equations presented above is that the vertical length and velocity scales are negligible with respect to the horizontal ones. Another assumption is the fact that the concentration of sediment (the ratio between the solid and liquid discharge) is small (below $6\times 10^{-3}$ (Garegnani, Rosatti & Bonaventura Reference Garegnani, Rosatti and Bonaventura2011, Reference Garegnani, Rosatti and Bonaventura2013)), such that we apply the clear water approximation.

2.2 Secondary flow equations

This section describes the equations that model secondary flow (i.e. formulations for $F_{sx}$ and $F_{sy}$ in (2.2) and (2.3)). The secondary flow velocity profile $u^{s}$ ( $\text{m}~\text{s}^{-1}$ ) (i.e. the vertical profile of the velocity component perpendicular to the primary flow) is assumed to have a universal shape as a function of the relative elevation in the water column $\unicode[STIX]{x1D701}=(z-\unicode[STIX]{x1D702})/h$ ( $-$ ), where $z$ (m) is the vertical Cartesian coordinate perpendicular to $x$ and $y$ increasing in the upward direction (Rozovskii Reference Rozovskii1957; Engelund Reference Engelund1974; De Vriend Reference de Vriend1977, Reference de Vriend1981; Booij & Pennekamp Reference Booij and Pennekamp1984). Worded differently, the vertical profile of the secondary flow is parametrised by a single value representing the intensity of the secondary flow $I$ ( $\text{m}~\text{s}^{-1}$ ), such that $u^{s}=f(\unicode[STIX]{x1D701})I$ . The secondary flow intensity $I$ is the integral of the absolute value of the secondary flow velocity profile (De Vriend Reference de Vriend1981). Among others, Rozovskii (Reference Rozovskii1957), Engelund (Reference Engelund1974) and De Vriend (Reference de Vriend1977), derive equilibrium profiles of the secondary flow that differ in the description of the eddy viscosity, vertical profile of the primary flow and the boundary condition of the flow at the bed. Following De Vriend (Reference de Vriend1977), we assume a logarithmic profile for the primary flow (i.e. a parabolic distribution of the eddy viscosity) and vanishing velocity close to the bed at $\unicode[STIX]{x1D701}=\exp (-1-1/\unicode[STIX]{x1D6FC})$ where $\unicode[STIX]{x1D6FC}=\sqrt{C_{f}}/\unicode[STIX]{x1D705}<0.5$ .

The depth-averaging procedure yields the integral value (along $z$ ) of the force per unit mass that the secondary flow exerts on the primary flow (De Vriend Reference de Vriend1977; Kalkwijk & De Vriend Reference Kalkwijk and de Vriend1980):

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle F_{sx}=\frac{\unicode[STIX]{x2202}T_{xx}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}T_{xy}}{\unicode[STIX]{x2202}y}, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle F_{sy}=\frac{\unicode[STIX]{x2202}T_{yx}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}T_{yy}}{\unicode[STIX]{x2202}y}, & \displaystyle\end{eqnarray}$$

where $T_{lm}$ ( $\text{m}^{3}~\text{s}^{-2}$ ) is the integral shear stress per unit mass in the direction $l$ - $m$ . Assuming a large width-to-depth ratio (i.e. $B/h\gg 1$ , where $B$ (m) is the characteristic channel width) and a mild curvature (i.e. $h/R_{s}\ll 1$ , where $R_{s}$ (m) is the radius of curvature of the streamlines), the shear stress terms are:

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle T_{xx}=-2\frac{\unicode[STIX]{x1D6FD}^{\ast }I}{Q}q_{x}q_{y}, & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle T_{xy}=T_{yx}=\frac{\unicode[STIX]{x1D6FD}^{\ast }I}{Q}(q_{x}^{2}-q_{y}^{2}), & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle T_{yy}=T_{yy}=2\frac{\unicode[STIX]{x1D6FD}^{\ast }I}{Q}q_{x}q_{y}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}^{\ast }=5\unicode[STIX]{x1D6FC}-15.6\unicode[STIX]{x1D6FC}^{2}+37.5\unicode[STIX]{x1D6FC}^{3}$ .

The simplest strategy to account for secondary flow assumes that the secondary flow is fully developed. This is equivalent to saying that the secondary flow intensity is equal to the equilibrium value $I_{e}=Q/R_{s}$ ( $\text{m}~\text{s}^{-1}$ ) found in an infinitely long bend (Rozovskii Reference Rozovskii1957; Engelund Reference Engelund1974; De Vriend Reference de Vriend1977, Reference de Vriend1981; Booij & Pennekamp Reference Booij and Pennekamp1983). A change in channel curvature leads to the streamwise adaptation of secondary flow to the equilibrium value (De Vriend Reference de Vriend1981; Ikeda & Nishimura Reference Ikeda and Nishimura1986; Johannesson & Parker Reference Johannesson and Parker1989; Seminara & Tubino Reference Seminara, Tubino, Ikeda and Parker1989). Booij & Pennekamp (Reference Booij and Pennekamp1984) and Kalkwijk & Booij (Reference Kalkwijk and Booij1986) not only account for the spatial adaptation but also the temporal adaptation of the secondary flow associated with a variable discharge or tides. Here we adopt the latter strategy, which has been applied, for instance, in modelling the morphodynamics of braided rivers (Javernick et al. Reference Javernick, Hicks, Measures, Caruso and Brasington2016; Williams et al. Reference Williams, Measures, Hicks and Brasington2016; Javernick, Redolfi & Bertoldi Reference Javernick, Redolfi and Bertoldi2018). The spatial and temporal adaptation of secondary flow is expressed by (Jagers Reference Jagers2003):

(2.11) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}I}{\unicode[STIX]{x2202}t}+\frac{q_{x}}{h}\frac{\unicode[STIX]{x2202}I}{\unicode[STIX]{x2202}x}+\frac{q_{y}}{h}\frac{\unicode[STIX]{x2202}I}{\unicode[STIX]{x2202}y}-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}I}{\unicode[STIX]{x2202}x}\right)-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}I}{\unicode[STIX]{x2202}y}\right)=S_{s}, & & \displaystyle\end{eqnarray}$$

where $S_{s}$ ( $\text{m}~\text{s}^{-2}$ ) is a source term which depends on the difference between the local secondary flow intensity and its equilibrium value:

(2.12) $$\begin{eqnarray}\displaystyle S_{s}=-\frac{I-I_{e}}{T_{I}}, & & \displaystyle\end{eqnarray}$$

where $T_{I}$ (s) is the adaptation time scale of the secondary flow:

(2.13) $$\begin{eqnarray}\displaystyle T_{I}=\frac{L_{I}h}{Q}, & & \displaystyle\end{eqnarray}$$

where $L_{I}=L_{I}^{\ast }h$ (m) is the adaptation length scale of the secondary flow, which depends on the non-dimensional length scale $L_{I}^{\ast }=(1-2\unicode[STIX]{x1D6FC})/2\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D6FC}$ (Kalkwijk & Booij Reference Kalkwijk and Booij1986).

The radius of curvature of the streamlines is defined as (e.g. Legleiter & Kyriakidis Reference Legleiter and Kyriakidis2006):

(2.14) $$\begin{eqnarray}\displaystyle \frac{1}{R_{s}}=\frac{\displaystyle \frac{\text{d}\hspace{1.0pt}x}{\text{d}t}\frac{\text{d}^{2}y}{\text{d}t^{2}}-\frac{\text{d}y}{\text{d}t}\frac{\text{d}^{2}x}{\text{d}t^{2}}}{\left(\left(\displaystyle \frac{\text{d}\hspace{1.0pt}x}{\text{d}t}\right)^{2}+\left(\displaystyle \frac{\text{d}y}{\text{d}t}\right)^{2}\right)^{3/2}}, & & \displaystyle\end{eqnarray}$$

assuming steady flow and in terms of water discharge we obtain:

(2.15) $$\begin{eqnarray}\displaystyle \frac{1}{R_{s}}=\frac{-q_{x}q_{y}\displaystyle \frac{\unicode[STIX]{x2202}q_{x}}{\unicode[STIX]{x2202}x}+q_{x}^{2}\frac{\unicode[STIX]{x2202}q_{y}}{\unicode[STIX]{x2202}x}-q_{y}^{2}\frac{\unicode[STIX]{x2202}q_{x}}{\unicode[STIX]{x2202}y}+q_{x}q_{y}\frac{\unicode[STIX]{x2202}q_{y}}{\unicode[STIX]{x2202}y}}{(q_{x}^{2}+q_{y}^{2})^{3/2}}. & & \displaystyle\end{eqnarray}$$

The secondary flow model described in this section closes the primary flow model described in § 2.1 given a certain bed elevation. In the following section we describe the model equations that describe changes in bed elevation as a function of the primary and secondary flow.

2.3 Morphodynamic equations

We consider an alluvial bed composed of an arbitrary number $N$ of non-cohesive sediment fractions characterised by a grain size $d_{k}$ (m), where the subscript $k$ denotes the grain size fraction in increasing order (i.e. $d_{1}<d_{2}<\cdots <d_{N}$ ). Bed elevation change depends on the divergence of the sediment transport rate (Exner Reference Exner1920):

(2.16) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}q_{bx}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}q_{by}}{\unicode[STIX]{x2202}y}=0, & & \displaystyle\end{eqnarray}$$

where $q_{bx}=\sum _{k=1}^{N}q_{bxk}$ ( $\text{m}^{2}~\text{s}^{-1}$ ) and $q_{by}=\sum _{k=1}^{N}q_{byk}$ ( $\text{m}^{2}~\text{s}^{-1}$ ) are the total specific (i.e. per unit of differential length) sediment transport rates including pores in the $x$ and $y$ direction, respectively. The variables $q_{bxk}$ ( $\text{m}^{2}~\text{s}^{-1}$ ) and $q_{byk}$ ( $\text{m}^{2}~\text{s}^{-1}$ ) are the specific sediment transport rates of size fraction $k$ including pores. For simplicity we assume a constant porosity and density of the bed sediment. The sediment transport rate is assumed to be locally at capacity, which implies that we do not model the temporal and spatial adaptation of the sediment transport rate to capacity conditions (Bell & Sutherland Reference Bell and Sutherland1983; Phillips & Sutherland Reference Phillips and Sutherland1989; Jain Reference Jain1992).

Changes in the bed surface grain size distribution are accounted for using the active layer model (Hirano Reference Hirano1971). For simplicity, we assume a constant active layer thickness $L_{a}$ (m). Conservation of sediment mass of size fraction $k$ in the active layer reads:

(2.17) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}M_{ak}}{\unicode[STIX]{x2202}t}+f_{k}^{I}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}q_{bxk}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}q_{byk}}{\unicode[STIX]{x2202}y}=0\quad k\in \{1,N-1\}, & & \displaystyle\end{eqnarray}$$

and in the substrate (Chavarrías et al. Reference Chavarrías, Stecca and Blom2018):

(2.18) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}M_{sk}}{\unicode[STIX]{x2202}t}-f_{k}^{I}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}t}=0\quad k\in \{1,N-1\}, & & \displaystyle\end{eqnarray}$$

where $M_{ak}=F_{ak}L_{a}$ (m) and $M_{sk}=\int _{\unicode[STIX]{x1D702}_{0}}^{\unicode[STIX]{x1D702}_{0}+\unicode[STIX]{x1D702}-L_{a}}f_{sk}(z)\,\text{d}z$ (m) are the volume of sediment of size fraction $k$ per unit of bed area in the active layer and the substrate, respectively. Parameter $\unicode[STIX]{x1D702}_{0}$ (m) is a datum for bed elevation. Parameters $F_{ak}\in [0,1]$ , $f_{sk}\in [0,1]$ and $f_{k}^{I}\in [0,1]$ are the volume fraction content of sediment of size fraction $k$ in the active layer, substrate and at the interface between the active layer and the substrate, respectively. By definition, the sum of the volume fraction content over all size fractions equals 1:

(2.19a-c ) $$\begin{eqnarray}\displaystyle \mathop{\sum }_{k=1}^{N}F_{ak}=1,\quad \mathop{\sum }_{k=1}^{N}f_{sk}(z)=1,\quad \mathop{\sum }_{k=1}^{N}f_{k}^{I}=1. & & \displaystyle\end{eqnarray}$$

Under degradational conditions, the volume fraction content of size fraction $k$ at the interface between the active layer and the substrate is equal to that at the top part of the substrate ( $f_{k}^{I}=f_{sk}(z=\unicode[STIX]{x1D702}-L_{a})$ for $\unicode[STIX]{x2202}\unicode[STIX]{x1D702}/\unicode[STIX]{x2202}t<0$ ). This allows for modelling of arbitrarily abrupt changes in grain size due to erosion of previous deposits. Under aggradational conditions the sediment transferred to the substrate is a weighted mixture of the sediment in the active layer and the bed load (Parker Reference Parker1991; Hoey & Ferguson Reference Hoey and Ferguson1994; Toro-Escobar, Paola & Parker Reference Toro-Escobar, Paola and Parker1996). Here we simplify the analysis and we assume that the contribution of the bed load to the depositional flux is negligible (i.e. $f_{k}^{I}=F_{ak}$ for $\unicode[STIX]{x2202}\unicode[STIX]{x1D702}/\unicode[STIX]{x2202}t>0$ ) (Hirano Reference Hirano1971).

The magnitude of the sediment transport rate is assumed to be a function of the local bed shear stress. We apply the sediment transport relation by Engelund & Hansen (Reference Engelund and Hansen1967) in a fractional manner (Blom, Viparelli & Chavarrías Reference Blom, Viparelli and Chavarrías2016; Blom et al. Reference Blom, Arkesteijn, Chavarrías and Viparelli2017) as well as the one by Ashida & Michiue (Reference Ashida and Michiue1971) (appendix B).

The direction of the sediment transport ( $\unicode[STIX]{x1D711}_{sk}$ (rad)) is affected by the secondary flow and the bed slope (Van Bendegom Reference van Bendegom1947):

(2.20) $$\begin{eqnarray}\displaystyle \tan \unicode[STIX]{x1D711}_{sk}=\frac{\sin \unicode[STIX]{x1D711}_{\unicode[STIX]{x1D70F}}-\displaystyle \frac{1}{g_{sk}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}y}}{\cos \unicode[STIX]{x1D711}_{\unicode[STIX]{x1D70F}}-\displaystyle \frac{1}{g_{sk}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}{\unicode[STIX]{x2202}x}}\quad k\in \{1,N\}, & & \displaystyle\end{eqnarray}$$

where $g_{sk}$ ( $-$ ) is a function that accounts for the influence of the bed slope on the sediment transport direction and $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D70F}}$ (rad) is the direction of the sediment transport accounting for the secondary flow only:

(2.21) $$\begin{eqnarray}\displaystyle \tan \unicode[STIX]{x1D711}_{\unicode[STIX]{x1D70F}}=\frac{q_{y}-h\unicode[STIX]{x1D6FC}_{I}\displaystyle \frac{q_{x}}{Q}I}{q_{x}-h\unicode[STIX]{x1D6FC}_{I}\displaystyle \frac{q_{y}}{Q}I}. & & \displaystyle\end{eqnarray}$$

Assuming a mild curvature, uniform flow conditions and a logarithmic profile of the primary flow, the constant $\unicode[STIX]{x1D6FC}_{I}$ ( $-$ ) is (De Vriend Reference de Vriend1977):

(2.22) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}_{I}=\frac{2}{\unicode[STIX]{x1D705}^{2}}(1-\unicode[STIX]{x1D6FC}). & & \displaystyle\end{eqnarray}$$

The effect of the bed slope on the sediment transport direction depends on the grain size (Parker & Andrews Reference Parker and Andrews1985). We account for this effect setting:

(2.23) $$\begin{eqnarray}\displaystyle g_{sk}=A_{s}\unicode[STIX]{x1D703}_{k}^{B_{s}}\quad k\in \{1,N\}, & & \displaystyle\end{eqnarray}$$

where $A_{s}$ ( $-$ ) and $B_{s}$ ( $-$ ) are non-dimensional parameters and $\unicode[STIX]{x1D703}_{k}$ ( $-$ ) is the Shields (Reference Shields1936) stress (appendix B). Different values of the coefficients $A_{s}$ and $B_{s}$ have been proposed (for a recent review, see Baar et al. (Reference Baar, de Smit, Uijttewaal and Kleinhans2018)). We consider two possibilities: (i) $A_{s}=1$ , $B_{s}=0$ (Schielen et al. Reference Schielen, Doelman and de Swart1993) and (ii) $A_{s}=1.70$ and $B_{s}=0.5$ (Talmon et al. Reference Talmon, Struiksma and Mierlo1995). In the first and simpler case, the bed slope effect is independent of the bed shear stress (Engelund & Skovgaard Reference Engelund and Skovgaard1973; Engelund Reference Engelund1975). In the second, more complex, case, the bed slope effect is assumed to be dependent on the fluid drag force on the grains, which is assumed to depend on the Shields stress (Koch & Flokstra Reference Koch and Flokstra1981).

2.4 Linearised system of equations

The system of equations describing the flow, change of bed level and change of the bed surface texture is highly nonlinear. Here we linearise the system of equations to provide insight on the fundamental properties of the model and to study the stability of perturbations. To this end we consider a reference state that is a solution to the system of equations. The reference state is a steady uniform straight flow in the $x$ direction over an inclined plane bed composed of an arbitrary number of size fractions. Mathematically: $h_{0}=\text{ct.}$ , $q_{x0}=\text{ct.}$ , $q_{y0}=0$ , $I_{0}=0$ , $\unicode[STIX]{x2202}\unicode[STIX]{x1D702}/\unicode[STIX]{x2202}x=\text{ct.}=-C_{f}q_{x0}^{2}/gh_{0}^{3}$ , $\unicode[STIX]{x2202}\unicode[STIX]{x1D702}/\unicode[STIX]{x2202}y=0$ , $M_{ak0}=\text{ct.}$ $\forall k\in \{1,N-1\}$ , where ct. denotes a constant different from 0 and subscript $0$ indicates the reference solution.

We add a small perturbation to the reference solution denoted by $^{\prime }$ and we linearise the resulting system of equations. After substituting the reference solution we obtain a system of equations of the perturbed variables:

(2.24) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{Q}^{\prime }}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D63F}_{\boldsymbol{x}\mathbf{0}}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{Q}^{\prime }}{\unicode[STIX]{x2202}x^{2}}+\unicode[STIX]{x1D63F}_{\boldsymbol{y}\mathbf{0}}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{Q}^{\prime }}{\unicode[STIX]{x2202}y^{2}}+\unicode[STIX]{x1D63C}_{\boldsymbol{x}\mathbf{0}}\frac{\unicode[STIX]{x2202}\boldsymbol{Q}^{\prime }}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D63C}_{\boldsymbol{y}\mathbf{0}}\frac{\unicode[STIX]{x2202}\boldsymbol{Q}^{\prime }}{\unicode[STIX]{x2202}y}+\unicode[STIX]{x1D63D}_{\mathbf{0}}\boldsymbol{Q}^{\prime }=0, & & \displaystyle\end{eqnarray}$$

where the vector of dependent variables is:

(2.25) $$\begin{eqnarray}\displaystyle \boldsymbol{Q}^{\prime }=[h^{\prime },q_{x}^{\prime },q_{y}^{\prime },I^{\prime },\unicode[STIX]{x1D702}^{\prime },[M_{ak}^{\prime }]]^{\text{T}}, & & \displaystyle\end{eqnarray}$$

where the square bracket indicates the vector character.

The diffusive matrix in $x$ direction is:

(2.26)

where $\mathbb{0}$ denotes the zero matrix. The diffusive matrix in the $y$ direction is:

(2.27)

The advective matrix in the $x$ direction is:

(2.28)

The advective matrix in the $y$ direction is:

(2.29)

The matrix of linear terms is:

(2.30)

We assume that the perturbations can be represented as a Fourier series, which implies that they are piecewise smooth and bounded for $x=\pm \infty$ . Using this assumption the solution of the perturbed system is expressed in the form of normal modes:

(2.31) $$\begin{eqnarray}\displaystyle \boldsymbol{Q}^{\prime }=\text{Re}(\boldsymbol{V}\text{e}^{\text{i}(k_{wx}+k_{wy}-\unicode[STIX]{x1D714}t)}), & & \displaystyle\end{eqnarray}$$

where $\text{i}$ is the imaginary unit, $k_{wx}$ ( $\text{rad}~\text{m}^{-1}$ ) and $k_{wy}$ ( $\text{rad}~\text{m}^{-1}$ ) are the real wavenumbers in the $x$ and $y$ direction, respectively, $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D714}_{i}$ ( $\text{rad}~\text{s}^{-1}$ ) is the complex angular frequency, $\boldsymbol{V}$ is the complex amplitude vector and Re denotes the real part of the solution (which we will omit in the subsequent steps). The variable $\unicode[STIX]{x1D714}_{r}$ is the angular frequency and $\unicode[STIX]{x1D714}_{i}$ the attenuation coefficient. A value of $\unicode[STIX]{x1D714}_{i}>0$ implies growth of perturbations and $\unicode[STIX]{x1D714}_{i}<0$ decay. Substitution of (2.31) in (2.24) yields:

(2.32) $$\begin{eqnarray}\displaystyle [\unicode[STIX]{x1D648}_{\mathbf{0}}-\unicode[STIX]{x1D714}\mathbb{1}]\boldsymbol{V}=0, & & \displaystyle\end{eqnarray}$$

where:

(2.33) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D648}_{\mathbf{0}}=\unicode[STIX]{x1D63F}_{\boldsymbol{x}\mathbf{0}}k_{wx}^{2}\text{i}+\unicode[STIX]{x1D63F}_{\boldsymbol{y}\mathbf{0}}k_{wy}^{2}\text{i}+\unicode[STIX]{x1D63C}_{\boldsymbol{x}\mathbf{0}}k_{wx}+\unicode[STIX]{x1D63C}_{\boldsymbol{y}\mathbf{0}}k_{wy}-\unicode[STIX]{x1D63D}_{\mathbf{0}}\text{i}, & & \displaystyle\end{eqnarray}$$

and $\mathbb{1}$ denotes the unit matrix. Equation (2.32) is an eigenvalue problem in which the eigenvalues of $\unicode[STIX]{x1D648}_{\mathbf{0}}$ (as a function of the wavenumber) are the values of $\unicode[STIX]{x1D714}$ satisfying (2.32).

The solution of the linear model provides information regarding the development of small amplitude oscillations only, but for an arbitrary wavenumber. For this reason the linear model is convenient for studying the well posedness of the model, which we will assess in the following section.

4.1 Ill posedness due to secondary flow

In this section we study how the stability of the system of equations is affected by the secondary flow model. To gain insight we compare three cases. In the first case we omit secondary flow. In the second and third cases we include the secondary flow model with and without considering diffusion (table 4).

Table 4. Variations to the reference state (table 1) and results of the linear analysis with respect to secondary flow.

The first case is equivalent to I1 (table 2), yet the eddy viscosity is equal to the value derived by Elder ((2.4), $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{E}=0.0057~\text{m}^{2}~\text{s}^{-1}$ ). In figure 3(a,b) we plot the maximum growth rate of perturbations as a function of the wavenumber and the wavelength, respectively. Diffusion appears to significantly dampen perturbations (compare figure 1(a) in which diffusion is neglected to figure 3 a).

Figure 3. Growth rate of perturbations added to the reference case (tables 1 and 4) as a function of the wavenumber and the wavelength: (a,b) without secondary flow (Case S1, well posed), (c,d) accounting for secondary flow with diffusion (Case S2, well posed) and (e,f) accounting for secondary flow without diffusion (Case S3, ill posed). The panels in the two columns show the same information but highlight the behaviour for large wavenumbers (a,c,e) and for large wavelengths (b,d,f). Red and green indicate growth and decay of perturbations, respectively.

In the second case we repeat the analysis including the equation for advection and diffusion of the secondary flow intensity (i.e. the first four rows and columns of matrix $\unicode[STIX]{x1D648}_{\mathbf{0}}$ in (2.33), Case S2, table 4). We observe that accounting for secondary flow introduces an instability mechanism (figure 3 d). For the specific conditions of the case, a growth domain appears for wavelengths between 0.7 m and 39 m long and between 0.4 m and 19 m wide. The maximum growth corresponds to a wavelength in the $x$ and $y$ direction equal to 1.29 m and 0.74 m, respectively. This situation is well posed, as for large wavenumbers perturbations decay (figure 3 c). Yet, the model is unsuitable for reproducing such instability, as it predicts growth of perturbations with a length scale of the order of the flow depth and shorter, for which the shallow water equation model is not suited. Given the fact that we consider a depth-averaged formulation of the primary flow, processes that scale with the flow depth are not resolved by the model and consequently perturbations at that scale must decay to yield physically realistic results. Otherwise, scales of the order of the flow depth become relevant, which contradicts the assumptions of the depth-averaged formulation. To model processes that scale with the flow depth such as dune growth, it is necessary to account for non-depth-averaged flow formulations that consider, for instance, rotational flow (Colombini & Stocchino Reference Colombini and Stocchino2011, Reference Colombini and Stocchino2012), or non-hydrostatic pressure (Giri & Shimizu Reference Giri and Shimizu2006; Shimizu et al. Reference Shimizu, Giri, Yamaguchi and Nelson2009).

In the third case we test the secondary flow model without accounting for diffusion in the system of equations ( $\unicode[STIX]{x1D708}=0$ , Case S3, table 4). We observe that the instability domain extends to infinitely large wavenumbers (figure 3 e), which implies that this model is ill posed (§ 3). We now aim to prove that the shallow water equations in combination with the secondary flow model without diffusion always yields an ill-posed model. To this end we obtain the characteristic polynomial of matrix  $\unicode[STIX]{x1D648}_{\mathbf{0}}$ (2.33). We compute the discriminant of the fourth-order characteristic polynomial and we find that for $k_{wx}<k_{wy}$ the growth rate of perturbations is positive (appendix C). The model is ill posed, as there always exists a domain of growth extending to infinitely large wavenumbers in the transverse direction.

We assess how the length scale of the instability related to the secondary flow model depends on the flow parameters. For this purpose we compute the shortest wave with positive growth for a varying diffusion coefficient and flow conditions (figure 4). We observe that, independently of the flow conditions, the theoretical value of the diffusion coefficient derived by Elder (Reference Elder1959) (2.4) is insufficient for dampening oscillations scaling with the flow depth. We conclude that if the diffusion coefficient is realistic, the treatment of the secondary flow yields an unrealistic model. It is necessary to use an unrealistically large value of the diffusion coefficient to obtain a realistic depth-averaged model in which perturbations scaling with the flow depth decay.

Figure 4. Wavelength of the shortest perturbation with positive growth rate ( $l_{wm}$ ) relative to the flow depth ( $h$ ) as a function of the Froude number ( $Fr$ ) and the diffusion coefficient ( $\unicode[STIX]{x1D708}$ ) relative to the diffusion coefficient according to Elder (Reference Elder1959) ( $\unicode[STIX]{x1D708}_{E}$ ). Different flow conditions are studied varying the flow depth between 0.2 m and 1.5 m from the reference case (table 1). The cyan markers indicate the conditions of three numerical simulations with different values of the diffusion coefficient (§ 5.1). The arrow next to the diamond marker indicates that the value lies outside the figure. Red (green) colour indicates that the shortest wavelength with positive growth rate are smaller (larger) than the flow depth.

4.2 Ill posedness due to bed slope effect

In this section we study the influence of considering the effect of the bed slope on model well posedness. To gain insight we compare five cases in which we consider unisize and mixed-size sediment, various sediment transport relations and various bed slope functions (table 5). We neglect secondary flow and diffusion to reduce the complexity of the problem (Parker Reference Parker1976; Fredsøe Reference Fredsøe1978; Colombini et al. Reference Colombini, Seminara and Tubino1987; Schielen et al. Reference Schielen, Doelman and de Swart1993).

Our reference case is B1 (§ 3) which considers unisize sediment conditions, and the effect of the bed slope on the sediment transport direction is accounted for using the simplest formulation, $g_{s}=1$ . We have shown that this case is well posed. Neglecting the effect of the bed slope on the sediment transport direction (Case B2, table 5) makes the problem ill posed (figure 5 a). This illustrates that accounting for the effect of the bed slope is required for obtaining not only physically realistic but also mathematically well-posed results. We prove that the shallow water equations in combination with the Exner (Reference Exner1920) equation without considering the effect of the bed slope always yields an ill-posed model by studying the growth rate of perturbations in the limit as the wavenumber $k_{wy}$ tends to infinity (appendix D).

Figure 5. Growth rate of perturbations added to the reference case (tables 1 and 5) as a function of the wavenumber and the wavelength: (a,b) Case B2 (ill posed), (c,d) Case B3 (well posed), (e,f) Case B4 (ill posed) and (g,h) Case B5 (ill posed). The panels in the two columns show the same information but highlight the behaviour for large wavenumbers (a,c,e,g) and for large wavelengths (b,d,f,h). Red and green indicate growth and decay of perturbations, respectively.

Table 5. Variations to the reference state (table 1) and results of the linear analysis with respect to the effect of the bed slope on the sediment transport direction. EH and AM refer to the sediment transport relations by Engelund & Hansen (Reference Engelund and Hansen1967) and Ashida & Michiue (Reference Ashida and Michiue1971), respectively.

The fact that the bed slope effect dampens perturbations under unisize conditions is expected from the fact that the only diffusive term in the system of equations is $\unicode[STIX]{x2202}q_{by}/\unicode[STIX]{x2202}s_{y}$ (2.27), where $s_{y}=\unicode[STIX]{x2202}\unicode[STIX]{x1D702}/\unicode[STIX]{x2202}y$ . This term is negative and approximately equal to $-q_{b}/g_{s}$ for a small streamwise slope. When we consider more than one grain size, diffusive terms appear in each active layer equation. We find that these diffusive terms may be positive, which hints at the possibility of an antidiffusive behaviour, which may lead to ill posedness. To study this possibility we compute the growth rate of perturbations of a simplified case consisting of two sediment size fractions. In the limit for the wavenumbers tending to infinity, the maximum growth rate is:

(4.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{i}^{lim}=\unicode[STIX]{x1D6FC}_{1}(r_{y1}-d_{x1,1})^{2}+\unicode[STIX]{x1D6FC}_{2}(r_{y1}-d_{x1,1})+\unicode[STIX]{x1D6FC}_{3}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}_{i}$ for $i=1,2,3$ are parameters relating the flow and the sediment transport rate (appendix E). The parameter $r_{y1}$ explains how the sediment transport rate of the fine fraction is affected by changes in the transverse bed slope:

(4.2) $$\begin{eqnarray}\displaystyle r_{y1}=\frac{\unicode[STIX]{x2202}q_{by1}/\unicode[STIX]{x2202}s_{y}}{\unicode[STIX]{x2202}q_{by}/\unicode[STIX]{x2202}s_{y}}, & & \displaystyle\end{eqnarray}$$

and the parameter $d_{x1,1}$ relates changes in the sediment transport rate to changes in the volume of sediment in the active layer:

(4.3) $$\begin{eqnarray}\displaystyle d_{x1,1}=\frac{\unicode[STIX]{x2202}q_{bx1}/\unicode[STIX]{x2202}M_{a1}}{\unicode[STIX]{x2202}q_{bx}/\unicode[STIX]{x2202}M_{a1}}. & & \displaystyle\end{eqnarray}$$

As $\unicode[STIX]{x1D6FC}_{1}>0$ (appendix E), there exists an interval $(r_{y1}-d_{x1,1})^{-}<(r_{y1}-d_{x1,1})<(r_{y1}-d_{x1,1})^{+}$ in which $\unicode[STIX]{x1D714}_{i}^{lim}<0$ and the model is well posed. Outside the interval, $\unicode[STIX]{x1D714}_{i}^{lim}>0$ and the problem is ill posed.

The physical interpretation of the limit values for obtaining a well-posed model is as follows. The effect of the transverse bed slope ( $r_{y1}$ ) needs to be balanced with respect to the effect of changes in surface texture ( $d_{x1,1}$ ) to obtain a well-posed model. For a given $d_{x1,1}$ , if parameter $r_{y1}$ is too small (i.e. the bed slope effect is not sufficiently strong) perturbations in the transverse direction are not dampened and the model is ill posed. On the other hand, for a given $r_{y1}$ , if parameter $d_{x1,1}$ is too small (e.g. due to relatively strong hiding or in conditions close to incipient motion) perturbations in the streamwise direction do not decay and the model is also ill posed. The critical values $r_{y1}^{\pm }$ that allow for a well-posed model are shown in appendix E.

In Cases B3–B5 we illustrate the possibility of ill posedness due to the bed slope closure relation (table 5). In Case B3 the sediment mixture consists of two grain size fractions with characteristic grain sizes equal to 0.001 m and 0.004 m. The volume fraction content of the fine sediment in the active layer and at the interface between the active layer and the substrate is equal to 0.5. The sediment transport rate is computed using the relation developed by Ashida & Michiue (Reference Ashida and Michiue1971). The other parameters are equal to the reference case (table 1). The system is well posed when the effect of the bed slope does not depend on the bed shear stress (figure 5 c). The situation is ill posed if the effect of the bed slope depends on the bed shear stress (Case B4, table 5, figure 5 e). The cause of the ill posedness is not found in the closure relation for the bed slope effect only but in the combination of the closure relation with the flow and bed surface conditions. A case equal to B3 except for the fact that the coarse grain size is equal to 0.012 m is ill posed (Case B5, table 5, figure 5 g).

We assess how the domain of ill posedness due to the bed slope effect depends on the model parameters. For given sediment sizes, flow and bed surface texture, parameter $B_{s}$ (2.23) controls the effect of the bed slope by modifying $r_{y1}$ only. The parameter $A_{s}$ (2.23) cancels in $r_{y1}$ and does not play a role. For this reason we study how $g_{s1}/A_{s}$ ( $-$ ) affects the domain of ill posedness for varying sediment properties, flow and bed surface grain size distribution (figure 6). We consider Case B3 and we vary $B_{s}$ between 0 and 0.5 to vary the bed slope effect. The sediment size of the coarse fraction varies between $d_{1}$ and 0.020 m. The mean flow velocity varies between $1~\text{m}~\text{s}^{-1}$ and $2.8~\text{m}~\text{s}^{-1}$ . The volume fraction content of fine sediment at the bed surface varies between 0 and 1. We aim to isolate the problem of ill posedness due to bed slope effect from the problem of ill posedness due to a different grain size distribution at the bed surface and at the interface between the bed surface and the substrate (Chavarrías et al. Reference Chavarrías, Stecca and Blom2018). For this reason, in this analysis the volume fraction content of fine sediment at the interface is equal to the one at the bed surface. Under this condition the problem can be ill posed due to the effect of the bed slope only.

As we have shown analytically, under unisize conditions (i.e. $d_{1}/d_{2}=1$ or $F_{a1}=1$ or $F_{a1}=0$ ) the model is well posed (figure 6 a,c). For sufficiently different grain sizes ( $d_{1}/d_{2}\lessapprox 0.15$ ) the model is well posed regardless of the bed slope effect but for a small range of values ( $0.08\lessapprox d_{1}/d_{2}\lessapprox 0.1$ ). This small range of ill-posed values is associated with $r_{y1}$ increasing for decreasing values of $d_{1}/d_{2}$ and acquiring a value larger than $r_{y1}^{+}$ such that the model becomes ill posed for all values of the bed slope effect. A further decrease in $d_{1}/d_{2}$ increases the limit value $r_{y1}^{+}$ faster than $r_{y1}$ such that the model becomes well posed for all values of the bed slope effect.

An increase in the Froude number decreases the domain of well posedness, which is explained from the fact that an increase in Froude number decreases $r_{y1}$ while it does not modify $r_{y1}^{-}$ . We have assumed steady flow in deriving $\unicode[STIX]{x1D714}_{i}^{lim}$ to reduce the complexity of the model such that we can find an analytical solution (appendix E). This causes a physically unrealistic resonance phenomenon for $Fr\rightarrow 1$ (Colombini & Stocchino Reference Colombini and Stocchino2005). In reality we do not expect that for $Fr=1$ the model is always ill posed regardless of the bed slope effect. Apart from the limit values in which the problem becomes unisize, the surface volume fraction content does not significantly affect the domain of ill posedness (figure 6 c) as it rescales in more or less a similar way $r_{y1}^{\pm }$ and $r_{y1}$ .

While Case B4 is ill posed because the effect of the bed slope ( $r_{y1}$ ) is small, Case B5 is ill posed because parameter $d_{x1,1}$ is small. The different origin of ill posedness does not cause a significant difference in the growth rate of perturbations as a function of the wavenumber (figure 5 e–g). However, we will find out that the pattern resulting from the perturbations depends on the origin of the ill posedness.

Figure 6. Domain of ill posedness due to the bed slope effect under mixed-size sediment conditions: as a function of the ratio between fine and coarse sediment (a), the Froude number (b) and the volume fraction content of fine sediment in the active layer (c). The bed slope effect is measured by $g_{s1}/A_{s}$ and the range of parameters is obtained by varying $B_{s}$ (2.23). The range of values of $d_{1}/d_{2}$ is obtained by varying $d_{2}$ . The range of values of the Froude number is obtained by varying $u$ . The volume fraction content of fine sediment at the interface between the active layer and the substrate is kept equal to the volume fraction content of fine sediment in the active layer. The conditions represent unisize sediment when $d_{1}/d_{2}=1$ , $F_{a1}=0$ , or $F_{a1}=1$ .

5.1 Secondary flow

We run five numerical simulations with a fixed bed (i.e. only the flow is computed) to study the roles of the secondary flow model and the diffusion coefficient in the ill posedness of the system of equations. The first three simulations reproduce the conditions of Cases S1, S2 and S3 (table 4). The domain is rectangular, 100 m long and 10 m wide. We use square cells with size equal to 0.1 m. The time step is equal to 0.01 s and we simulate 10 min of flow. We set a constant water discharge and the downstream water level remains constant with time. The initial condition represents normal flow for the values in table 1 (i.e. equilibrium conditions).

The simulation not accounting for secondary flow does not present growth of perturbations as predicted by the linear analysis and remains stable with time (figure 7 a). We observe growth of perturbations when we account for secondary flow with the diffusion coefficient derived by Elder (Reference Elder1959) (figure 7 b). The results are physically unrealistic as the flow depth presents variations of up to 7 % of the normal flow depth and the length scale of the perturbations is smaller than the flow depth. We have not introduced any perturbation in the initial or boundary conditions which implies that perturbations grow from numerical truncation errors. This supports the fact that the simulation is physically unrealistic. The case with a diffusion coefficient equal to 0 is ill posed and the solution presents unreasonably large oscillations (figure 7 c). These numerical results confirm the results of the linear stability analysis.

In the fourth simulation we set a diffusion coefficient 100 times larger than the one derived by Elder (Reference Elder1959) (figure 7 d). Under this condition the linear analysis predicts all short waves to decay (diamond in figure 4). These numerical results confirm the linear theory. The last simulation is equal to Case S2 except for the fact that we use a coarser grid ( $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=1~\text{m}$ ). In this case the numerical grid is not sufficiently detailed to resolve the perturbations due to secondary flow and the simulation is stable (figure 7 e). This last case highlights an important limitation of a physically unrealistic model. Although Case S2 is mathematically well posed, the solution presents similarities with ill-posed cases in the sense that a refinement of the grid causes non-physical oscillations to appear.

Figure 7. Flow depth at the end of the simulations: (a) without accounting for secondary flow (Case S1), (b) setting $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{E}$ (Case S2), (c) setting $\unicode[STIX]{x1D708}=0$ (Case S3), (d) setting $\unicode[STIX]{x1D708}=100\unicode[STIX]{x1D708}_{E}$ and (e) setting $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{E}$ using a coarser numerical grid (Case S2). The colour map is adjusted for each case and centred on the initial and equilibrium values ( $h=1~\text{m}$ ).

5.2 Bed slope effect

In this section we focus on the consequences of accounting for the bed slope effect on the mathematical character of the model. To this end we run five more numerical simulations without accounting for secondary flow and updating the bed (i.e. accounting for morphodynamic change). The simulations reproduce Cases B1–B5 (table 5). We simulate 8 h using a time step $\unicode[STIX]{x0394}t=0.1~\text{s}$ .

We have proved that accounting for the effect of the bed slope makes a unisize simulation well posed (§ 4.2 and figure 1 c). The numerical solution of this case (B1, table 2) is stable and perturbations do not grow (figure 8 a). Moreover, no alternate bars appear as the channel width is below the critical value (§ 3). Perturbations grow when the effect of the bed slope is not taken into account (Case B2, figure 8 b), which confirms that this case is ill posed.

Figure 8. Flow depth at the end of the simulations of: (a) Case B1, (b) Case B2; and volume fraction content of fine sediment in the active layer: (c) Case B3, (d) Case B4, (e) Case B5. The colour map is adjusted for each case and centred on the initial and equilibrium values.

The mixed-size sediment conditions of Case B3 yield a well-posed model (figure 5 e) and the simulation is stable (figure 8 c). On the other hand, the ill-posed cases B4 and B5 present growth of unrealistic and non-physical perturbations (figure 8 d,e). While the growth of perturbations in Case B5 seems random, in Case B4 we observe a clear pattern. Moreover, oscillations have grown significantly faster in Case B5 than in Case B4. While after 8 h the changes in volume fraction content at the bed surface are of the order of $10^{-3}$ in Case B4, these are of order 1 in Case B5.

The fact that oscillations grow faster in Case B5 than in Case B4 is related to the different origin of ill posedness. While Case B4 is ill posed because the effect of the bed slope is not sufficiently strong (i.e. $r_{y1}<r_{y1}^{-}$ ), Case B5 is ill posed because changes in the sediment transport rate due to changes in the volume of fine sediment in the active layer are too small (i.e. $r_{y1}>r_{y1}^{+}$ ). The first case is closely linked to the lateral direction, in which sediment transport is initially zero. The fact that initially the lateral sediment transport rate is zero limits the rate at which lateral changes occur. In the second case perturbations are linked to the streamwise direction, in which the sediment transport rate initially is non-zero, which enhances the rate at which changes develop.