3.1. Hydraulic limit

We begin our analysis by briefly recapping the ‘purely hydraulic’ stability problem within our framework. That is, we address the limiting case of weak morphodynamic processes, by sending both $Q \to 0$ and $\varGamma \to 0$. In this case, perturbations in $\psi$ and $b$ can only be advected along the slope, since there are no morphodynamic feedbacks through which they may grow or decay. Equation (3.5) possesses the solutions $\sigma = -\mathrm {i} k$ and $\sigma = 0$, that respectively correspond to these modes of disturbance. The remaining two solutions are

(3.8)\begin{equation} \sigma ={-}\mathrm{i} k - \frac{\tau_{u_0}}{2} \pm \sqrt{ \tau_{u_0}^2 / 4 - k^2 / {\textit{Fr}}^2 + \mathrm{i} k (\tau_{h_0} - 1) }. \end{equation}

These branches correspond to disturbances in the hydraulic governing equations for $h$ and $u$, studied in the case of general drag by Trowbridge (Reference Trowbridge1987). When $k = 0$, they pass through $\sigma = -\tau _{u_0}$ and $0$. It can be shown straightforwardly that $\mathrm {Re}(\sigma )$ is a monotonic function with respect to $|k|$, meaning that the maximum growth for each branch must occur at either $k = 0$, or in the limit $|k|\to \infty$. Growth rate saturation at short wavelengths is a known property of the classical roll-wave instability that highlights the omission of physics (e.g. turbulent dissipation) that would otherwise damp out disturbances over short length scales. Evaluating the limit of (3.8) as $|k|\to \infty$ yields

(3.9)\begin{equation} \mathrm{Re}(\sigma) \to \frac{-\tau_{u_0} \pm |1 - \tau_{h_0}|{\textit{Fr}}}{2}. \end{equation}

If $\tau _{u_0} < 0$, then there is always unstable growth (i.e., at $k=0$). However, we consider the more physically reasonable situation where $\tau _{u_0} > 0$ (i.e., a drag parametrisation that increases resistance to flow at higher shear rates). Then, if $\tau _{h_0} = 1$, both branches are everywhere stable and asymptote to $\mathrm {Re}(\sigma ) = -\tau _{u_0}/2$. Otherwise, since the argument of the square root in (3.8) always has a non-zero imaginary part (away from $k = 0$), the growth rates are always distinct and in particular, the branch with positive root always dominates. This turns unstable when (3.9) exceeds zero, which occurs if

(3.10)\begin{equation} {\textit{Fr}} > \frac{\tau_{u_0}}{|1 - \tau_{h_0}|}. \end{equation}

This is the stability criterion due to Trowbridge (Reference Trowbridge1987), written in our dimensionless quantities. Inclusion of the absolute value in the denominator constitutes a minor correction to the original formula that accounts for the case where $\tau _{h_0} > 1$.

3.2. Suspended load model

We now reintroduce morphodynamics, by allowing for non-vanishing mass exchange with the bed ($\varGamma \neq 0$), but continuing to neglect bed load transport ($Q = 0$). This substantially complicates (3.5), which becomes a fully $4\times 4$ problem. Motivated by the above discussion, we divide our morphodynamic analysis into two tractable regimes: the long-wave (or global) limit $k = 0$ and the short-wave limit $k \gg 1$, and verify later that these limits control most of the important aspects of the problem.

3.2.1. Global modes: $k = 0$

A given steady morphodynamic flow is specified by four state variables $\tilde {h}_0$, $\tilde {u}_0$, $\tilde {\psi }_0$ and $\tilde {b}_0$, which are constrained by only two (3.1a,b). Therefore, the solution space is underdetermined and there is a two-dimensional linear family of possible steady states. In nature, selection of a particular flow from this family is assured via some boundary condition, such as the total flux of material through a flow cross-section. Moreover, transitions from one steady flow to another within this space can occur (e.g. through an increase in the total flux). Infinitesimal transitions between steady states are linear perturbations in the sense of (3.4a)–(3.4d), with $k=0$ and $\sigma = 0$ (neutral stability). Therefore, by (3.5) they satisfy $\boldsymbol{\mathsf{C}}\boldsymbol {q} = \boldsymbol {0}$. Solving for $\boldsymbol {q}$ reveals a two-dimensional space of neutral modes spanned by

(3.11a,b) \begin{align} \boldsymbol{v}_1 = \begin{pmatrix} (\tau_{\psi_0} - \Delta\rho) \varGamma_{u_0} - \tau_{u_0}\varGamma_{\psi_0} \\ (\Delta\rho - \tau_{\psi_0}) \varGamma_{h_0} + (\tau_{h_0} - 1)\varGamma_{\psi_0} \\ \tau_{u_0}\varGamma_{h_0} - (\tau_{h_0} - 1)\varGamma_{u_0} \\ 0 \end{pmatrix},\quad \boldsymbol{v}_2 = \boldsymbol{e}_4, \end{align}

where we adopt the convention of using $\boldsymbol {e}_j$ to denote the $j$th standard basis vector. The first of these, $\boldsymbol {v}_1$, may be interpreted in the following way. Written in our dimensionless variables, the equations for steady flows (3.1a,b) are the roots of the function $\boldsymbol {F}(h, u, \psi ) = (\tau - \rho h, \varGamma )^{\textrm {T}}$. It is straightforward to verify that $\boldsymbol {\nabla } \boldsymbol {F}(h_0,u_0,\psi _0) \boldsymbol {\cdot } \boldsymbol {v}_1 = \boldsymbol {0}$ and therefore $\boldsymbol {v}_1$ represents a shift along the curve of solutions, implicitly defined by $\boldsymbol {F} = \boldsymbol {0}$. The second neutral mode $\boldsymbol {v}_2$ accounts for invariance to arbitrary translations of the bed height.

The remaining two global modes have non-zero growth rate and therefore, by (3.5), they obey

(3.12)\begin{equation} \sigma\boldsymbol{\mathsf{A}}\boldsymbol{q} + \boldsymbol{\mathsf{C}}\boldsymbol{q} = \boldsymbol{0}. \end{equation}

After factoring out the neutral growth rates, the characteristic equation yields a quadratic from which the remaining two eigenvalues may be directly computed. The full set of eigenvalues of (3.12) is then

(3.13a,b)\begin{equation} \sigma = 0\ \mathrm{(repeated)}, \quad \frac{s_0 \pm \sqrt{s_c}}{2}, \end{equation}

where $s_0$, $s_c$ are placeholders for

(3.14a)\begin{gather} s_0 ={-}\tau_{u_0} + \varGamma_{h_0} + (\upsilon_0 - \rho_b) \varGamma_{u_0} + (1 - \psi_0)\varGamma_{\psi_0}, \end{gather}

(3.14b)\begin{gather}s_c = \varGamma_{u_0}^2 (\upsilon_0 - \rho_b)^2 + 2\varGamma_{u_0} \left\{ (\upsilon_0 - \rho_b) \left[ \varGamma_{h_0} + (1-\psi_0)\varGamma_{\psi_0} - \tau_{u_0} \right] \right. \nonumber\\ \qquad \qquad \qquad - \left. 2\left[\tau_{h_0} + \tau_{\psi_0}(1 - \psi_0) - \rho_b\right] \right\} + \left[ \tau_{u_0} + \varGamma_{h_0} + (1 - \psi_0)\varGamma_{\psi_0} \right]^2. \end{gather}

Here, we have made use of (3.7) with $\psi = \psi _0$ and $\psi = \psi _b = 1$, to eliminate $\Delta \rho$ in favour of the bed density $\rho _b = 1 + \Delta \rho (1 - \psi _0)$ in these expressions, which nevertheless depend on all nine independent quantities in the matrices $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{C}}$. Before moving on to the next section, we note two important special cases.

In the non-erosive limit $\varGamma \to 0$, (3.14a) and (3.14b) reduce to simply $s_0 = -\tau _{u_0}$ and $s_c = \tau _{u_0}^2$. Substituting these into (3.13b) leaves only one (typically negative) non-zero growth rate, $\sigma = -\tau _{u_0}$, consistent with the analysis in § 3.1.

If instead, $\varGamma$ is finite, but $|\varGamma _{u_0}|$ is sufficiently small, relative to the other components of (3.14a,b), so that it may be neglected, the non-zero eigenvalues become

(3.15a,b)\begin{equation} \sigma = \frac{s_0 - \sqrt{s_c}}{2} ={-}\tau_{u_0}\quad \mathrm{and}\quad \sigma= \frac{s_0 + \sqrt{s_c}}{2} = \varGamma_{h_0} + \varGamma_{\psi_0}(1 - \psi_0). \end{equation}

Since the latter eigenvalue (later referred to as $\sigma _a$) may be positive, there exists a route to a purely morphodynamic instability in this case, which depends on the signs and relative magnitudes of $\varGamma _{h_0}$ and $\varGamma _{\psi _0}$. Positive values for these derivatives imply positive morphodynamic feedbacks, amplifying the flow depth and concentration respectively. We return to this in § 4, where we demonstrate using some generic model closures that this mode can indeed be unstable.

3.2.2. Short wavelengths: $k \gg 1$

We now focus on short-wavelength perturbations. By analogy with the non-morphodynamic case of § 3.1, we anticipate that the limit $k\to \infty$ controls the onset and growth of instabilities by maximising $\mathrm {Re}[\sigma (k)]$. (We confirm that this is often the case for example model closures in § 4.) The form of (3.5) suggests the following asymptotic expansions for the four growth rates and their corresponding eigenmodes in this regime:

(3.16a,b)\begin{equation} \sigma ={-}\mathrm{i}\lambda_1 k + \lambda_0 + \lambda_{{-}1} k^{{-}1} + \ldots, \quad \boldsymbol{q} = \boldsymbol{q}_0 + \boldsymbol{q}_{{-}1} k^{{-}1} + \ldots. \end{equation}

Here, $\lambda _{1}$, $\lambda _0$, $\lambda _{-1}$ and $\boldsymbol {q}_0$, $\boldsymbol {q}_{-1}$, are unknown constants and vectors to be determined shortly. Substituting these expressions into (3.5) and retaining only the leading $O(k)$ terms leaves an eigenproblem for $\lambda _1$

(3.17)\begin{equation} \lambda_1\boldsymbol{\mathsf{A}}\boldsymbol{q}_0 = \boldsymbol{\mathsf{B}}\boldsymbol{q}_0. \end{equation}

This may be solved to obtain four distinct values

(3.18)\begin{equation} \lambda_1 = 1 \pm {\textit{Fr}}^{{-}1},\ 1,\ 0. \end{equation}

Since $c = -\mathrm {Im}(\sigma ) / k \to \lambda _1$ as $k\to \infty$, these are the wave speeds for disturbances in the short-wavelength regime (and also the characteristics of the governing equations in this context). The corresponding eigenvectors of (3.17) are

(3.19)\begin{equation} \boldsymbol{q}_0 = \begin{pmatrix} \pm {\textit{Fr}} \\ 1 \\ 0\\ 0 \end{pmatrix},\quad \begin{pmatrix} \Delta \rho/2 \\ 0 \\ -1 \\ 0 \end{pmatrix},\quad \begin{pmatrix} 1 \\ -1 \\ 0\\ {\textit{Fr}}^2 - 1 \end{pmatrix} . \end{equation}

Recalling the definition $\boldsymbol {q}=(h_1,u_1,\psi _1,b_1)^{\textrm {T}}$ and (3.4a)–(3.4d), the elements of these vectors are the leading-order amplitudes for each mode. Throughout the rest of the paper, we label these modes I–IV. Since these asymptotic vectors separate the four solution branches of the general linear problem (3.5), it will be convenient later to use the same labels to refer to quantities at finite $k$, though they may not necessarily share the properties of their asymptotic counterparts.

The first pair of modes (I,II) in (3.19) contain no morphodynamic content. Indeed, they are identical to the short-wavelength modes of the purely hydraulic problem (§ 3.1), which is guaranteed since $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$ do not depend on $\varGamma$. They describe disturbances in $h$ and $u$, propagating at speeds $c = 1 \pm {\textit {Fr}}^{-1}$. Mode III couples unit speed perturbations in $\psi$ with the flow free surface, while mode IV is stationary ($c=0$) and disturbs the bedform, as well as $h$ and $u$.

The second term in the expansion of $\sigma$ determines the leading-order real part of the growth rate. We substitute (3.16a,b) back into (3.5) and subtract away the $O(k)$ component, i.e., (3.17). Retaining only $O(1)$ terms in the remaining equation, leaves

(3.20)\begin{equation} \lambda_0 \boldsymbol{\mathsf{A}}\boldsymbol{q}_0 + \mathrm{i} (\boldsymbol{\mathsf{B}} - \lambda_1 \boldsymbol{\mathsf{A}})\boldsymbol{q}_{{-}1} = \boldsymbol{\mathsf{C}}\boldsymbol{q}_0. \end{equation}

The unknown vector $\boldsymbol {q}_{-1}$ can be eliminated by solving the eigenproblem adjoint to (3.17), which yields vectors $\boldsymbol {r}_0$ such that $\lambda _1 \boldsymbol {r}_0^{\textrm {T}} \boldsymbol{\mathsf{A}} = \boldsymbol {r}_0^{\textrm {T}} \boldsymbol{\mathsf{B}}$. Multiplying (3.20) on the left by $\boldsymbol {r}_0^{\textrm {T}}$ and rearranging gives the formula

(3.21)\begin{equation} \lambda_0 ={-}\frac{\boldsymbol{r_0}\boldsymbol{\cdot}\boldsymbol{\mathsf{C}}\boldsymbol{q}_0}{\boldsymbol{r}_0\boldsymbol{\cdot}\boldsymbol{\mathsf{A}}\boldsymbol{q}_0}. \end{equation}

Using this, the following four expressions for $\lambda _0$ are obtained, which we label $\lambda _{0,1}, \ldots , \lambda _{0,4}$ for later reference:

(3.22a,b)\begin{gather} \lambda_{0,1} = \frac{f_+({\textit{Fr}})}{{\textit{Fr}}({\textit{Fr}} + 1)},\quad \lambda_{0,2} = \frac{f_-({\textit{Fr}})}{{\textit{Fr}}({\textit{Fr}} - 1)}, \end{gather}

(3.22c,d)\begin{gather}\lambda_{0,3} = (1 - \psi_0)(\varGamma_{\psi_0} - \Delta\rho \varGamma_{h_0} / 2),\quad \lambda_{0,4} = \frac{\varGamma_{u_0} - \varGamma_{h_0}}{{\textit{Fr}}^2 - 1}. \end{gather}

By (3.16a), these asymptotic values dictate the limits of $\mathrm {Re}(\sigma )$ as $k\to \infty$ for modes I–IV. We list them in the same order as their respective wave speeds in (3.18) and the $O(1)$ eigenvectors in (3.19). The functions $f_\pm$ are third-order polynomials in ${\textit {Fr}}$ defined by

(3.23) \begin{align} f_\pm({\textit{Fr}}) &={\pm}\tfrac{1}{2}[ \varGamma_{h_0}(\upsilon_0 - \rho_b) + (1 - \tau_{h_0}) ] {\textit{Fr}}^3 \nonumber\\ &\quad +\tfrac{1}{2}[ \varGamma_{h_0} (2\upsilon_0 - \rho_b + 1) / 2 + \varGamma_{u_0} (\upsilon_0 - \rho_b) + 1 - \tau_{h_0} - \tau_{u_0} ] {\textit{Fr}}^2 \nonumber\\ &\quad \pm\tfrac{1}{4}[ \varGamma_{h_0}(\rho_b - 1) + \varGamma_{u_0} (2\upsilon_0 - \rho_b + 1) - 2\tau_{u_0} ] {\textit{Fr}} + \tfrac{1}{4} \varGamma_{u_0} (\rho_b - 1). \end{align}

It is easily confirmed that as $\varGamma \to 0$, $\lambda _{0,1}$ and $\lambda _{0,2}$ reduce to the high-$k$ growth rates of the non-erosive problem, given in (3.9), while $\lambda _{0,3}, \lambda _{0,4}\to 0$. For this reason, we will sometimes label $\lambda _{0,1}$, $\lambda _{0,2}$ and their corresponding modes (I,II) as ‘hydraulic’ and $\lambda _{0,3}$, $\lambda _{0,4}$ as ‘morphodynamic’ even though all of (3.22a–d) are coupled to the bed and sediment dynamics when $\varGamma$ is non-vanishing.

Just as in the non-erosive problem, the asymptotic growth rates in (3.22a–d) are non-zero, but typically finite. However, there is an extra complication. Since $f_\pm (0) = \varGamma _{u_0} (\rho _b - 1) / 4$ and $f_\pm (\mp 1) = (\varGamma _{h_0} - \varGamma _{u_0})/2$, the pairs $(\lambda _{0,1}, \lambda _{0,2})$ and $(\lambda _{0,2}, \lambda _{0,4})$ possess singularities at ${\textit {Fr}} = 0$ and ${\textit {Fr}} = 1$ respectively (provided $\varGamma _{h_0} \neq \varGamma _{u_0} \neq 0$).

When ${\textit {Fr}} = 0$, the steady flow velocity is zero. The singularities in the expressions for $\lambda _{0,1}$ and $\lambda _{0,2}$ are artefacts arising from the fact that the time scale chosen to non-dimensionalise (3.3a–d) vanishes in the limit $\tilde {u}_0 \to 0$. Referring back to (3.2b,h,l), we may rewrite (3.22a,b) in dimensional units and verify that these growth rates remain finite. Specifically, when $\tilde {u}_0 = 0$, the expressions are $\tilde {\lambda }_{0,j} = (3/4-j/2)\tilde {\varGamma }_{\tilde {u}_0}(\tilde {\rho }_b/\tilde {\rho }_0 - 1)(g \cos \phi / \tilde {h}_0)^{1/2}$, for $j = 1,2$.

However, the singularities in (3.22b,d) at unit Froude number cannot be removed by a choice of units. They occur when the wave speeds $\lambda _1 = 1 - {\textit {Fr}}^{-1}, 0$ for disturbances to the flow and bedform, coalesce, as do the corresponding $O(1)$ modes in (3.19). Since $\lambda _0$ cannot be $O(1)$ at this singular point, our expansions in (3.16a,b) are inappropriate here. Therefore, we propose instead that at ${\textit {Fr}} = 1$ (and for modes II and IV only), $\sigma$ and $\boldsymbol {q}$ take the asymptotic form

(3.24a,b)\begin{equation} \sigma = \lambda_{1/2} k^{1/2} + \lambda_0 + \ldots, \quad \boldsymbol{q} = \boldsymbol{q}_0 + \boldsymbol{q}_{{-}1/2}k^{{-}1/2} + \boldsymbol{q}_{{-}1}k^{{-}1} + \ldots, \end{equation}

where $\lambda _{1/2}$, $\lambda _0$ and $\boldsymbol {q}_0$, $\boldsymbol {q}_{-1/2}$, $\boldsymbol {q}_{-1}$ are to be determined. We proceed as before, by substituting these expressions into (3.5) and isolating its constituent parts at different orders in $k$. Retaining only $O(k)$ terms yields $\boldsymbol{\mathsf{B}}\boldsymbol {q}_0 = \boldsymbol{0}$, with only one solution, $\boldsymbol {q}_0 = (1, -1, 0, 0)^{\textrm {T}}$. As it must, this matches the coalescent modes in (3.19), when they are evaluated at ${\textit {Fr}} = 1$. At $O(k^{1/2})$, we have

(3.25)\begin{equation} \lambda_{1/2}\boldsymbol{\mathsf{A}}\boldsymbol{q}_0 + \mathrm{i}\boldsymbol{\mathsf{B}}\boldsymbol{q}_{{-}1/2} = \boldsymbol{0}. \end{equation}

On substituting $\boldsymbol {q}_0$ into this equation, a little algebra shows that $\boldsymbol {e}_4 \boldsymbol {\cdot } \boldsymbol {q}_{-1/2} = 2 \mathrm {i} \lambda _{1/2}$. To find $\lambda _{1/2}$, we use the $O(1)$ equation, which is

(3.26)\begin{equation} \lambda_0 \boldsymbol{\mathsf{A}}\boldsymbol{q}_0 + \lambda_{1/2} \boldsymbol{\mathsf{A}} \boldsymbol{q}_{{-}1/2} + \mathrm{i} \boldsymbol{\mathsf{B}}\boldsymbol{q}_{{-}1} + \boldsymbol{\mathsf{C}} \boldsymbol{q}_0 = \boldsymbol{0}. \end{equation}

Now we notice that $\boldsymbol {e}_4 \boldsymbol {\cdot } \boldsymbol{\mathsf{A}} \boldsymbol {q}_0 = \boldsymbol {e}_4 \boldsymbol {\cdot } \boldsymbol{\mathsf{B}} \boldsymbol {v} = 0$ for any vector $\boldsymbol {v}$. Therefore, projecting (3.26) onto $\boldsymbol {e}_4$ eliminates the unknowns $\lambda _0$ and $\boldsymbol {q}_{-1}$. On doing this, substituting our expressions for $\boldsymbol {q}_0$ and $\boldsymbol {e}_4 \boldsymbol {\cdot } \boldsymbol {q}_{-1/2}$ from above and rearranging, we find

(3.27)\begin{equation} \lambda_{1/2} ={\pm} \frac{1 + \mathrm{i}}{2} (\varGamma_{h_0} - \varGamma_{u_0})^{1/2}. \end{equation}

Except for the particular case $\varGamma _{h_0} = \varGamma _{u_0}$ (where there is no singularity in $\lambda _{0,2},\lambda _{0,4}$), these expressions have non-zero real part. Therefore, at ${\textit {Fr}} = 1$ and for $k\gg 1$, the second and fourth modes (3.19) of the linear stability problem (3.5) diverge with amplitudes $\sim \exp (\pm A \sqrt {k}t)$, where $A = |\mathrm {Re}(\lambda _{1/2})|$. Crucially, one of these amplitudes is strictly positive and unbounded in the limit $k \to \infty$. This implies that the morphodynamic governing equations (3.3a–d) at are ill posed as an initial value problem when ${\textit {Fr}} = 1$, because their solutions do not depend continuously on initial data (Joseph & Saut Reference Joseph and Saut1990). Mathematically speaking, this is a direct consequence of the model losing the property of strict hyperbolicity when its characteristic wave speeds (3.18) intersect. More intuitively, problems arise because over any finite time interval, there are short-wavelength disturbances that grow arbitrarily rapidly, making it impossible for the governing equations to behave in a physically consistent way. This fact has practical consequences beyond the theory of steady flows on constant slopes. Computer simulations of these models (in both one and two spatial dimensions) conducted on complex topographies almost inevitably feature locations where the conditions locally match our problem at unit Froude number and shortwave oscillations can grow catastrophically. Numerical ‘solutions’ in this case may nevertheless look physically reasonable, since spatial discretisation imposes an upper limit on $k$. However, they will not converge as the numerical resolution increases and cannot be relied upon to model real flows.

The essential issue of unbounded growth rates in these models was recognised by Balmforth & Vakil (Reference Balmforth and Vakil2012), who studied the stability of a similar, but non-equivalent system: uniform flows eroding at a constant positive rate in the Saint-Venant equations. In the limit of slow erosion, they also observed that the high-wavenumber growth rate of perturbations suffers a singularity at unit Froude number. Moreover, they were able to show that the inclusion of a diffusive term in the momentum dynamics was sufficient to regularise their system. Our more general setting adds dynamic coupling with the solid phase and an arbitrary basal drag parametrisation, thereby demonstrating that the same problem affects a far broader range of shallow morphodynamic flow models. Indeed, it suggests that in any situations close to, but not strictly covered by our framework, it is important to check carefully whether the governing equations are well posed and amend them if necessary. Therefore, we continue with an analysis of how this might be achieved.

3.3. Regularisation

We shall introduce a new term to (3.3c) in order to quash unbounded growth at small length scales. As noted by Joseph & Saut (Reference Joseph and Saut1990), ill posedness often signals that there are physical processes missing from a model. In our case, two possible culprits are the shallow-layer approximation and the omission of bed load from the current analysis. We assess the effect of bed load shortly, in § 3.4. A particular effect neglected by the assumption of shallow flow is the aggregate loss of horizontal momentum caused by turbulent eddies. This is usually acceptable, since it is only significant at length scales shorter than the flow depth. However, for short waves it is no longer strictly negligible. A simple and common way to include this missing physics is to try to capture it via diffusion-like process. We denote a characteristic eddy viscosity for the flow by $\tilde {\nu }$ and non-dimensionalise by setting $\nu = \tilde {\nu } / (\tilde {u}_0 \tilde {l}_0)$. This free parameter sets the scale of the diffusive term $({\partial }/{\partial x})(\nu \rho h ({\partial u}/{\partial x}))$, which we add to the right-hand side of (3.3c). We note that the extra term does not affect the steady uniform layer itself. Similar expressions have been employed elsewhere, as a regularisation term by Balmforth & Vakil (Reference Balmforth and Vakil2012) in their analysis and in the shallow-flow models of Simpson & Castelltort (Reference Simpson and Castelltort2006), Xia et al. (Reference Xia, Lin, Falconer and Wang2010) and Langendoen et al. (Reference Langendoen, Mendoza, Abad, Tassi, Wang, Ata, El kadi Abderrezzak and Hervouet2016). Later, in § 4.6 we briefly address the implications of adding a similar term to (3.3b) to encapsulate turbulent sediment diffusivity.

It is unclear a priori whether the eddy viscosity term is sufficient to regularise the ill-posed model equations on its own. Therefore, we must extend the high-wavenumber growth rate analysis of § 3.2.2. With the extra term, the linearised system of (3.5) generalises to

(3.28)\begin{equation} \sigma\boldsymbol{\mathsf{A}}\boldsymbol{q} + \mathrm{i} k \boldsymbol{\mathsf{B}}\boldsymbol{q} + \boldsymbol{\mathsf{C}}\boldsymbol{q} ={-}k^2 \boldsymbol{\mathsf{D}}\boldsymbol{q}, \end{equation}

where $\boldsymbol{\mathsf{D}} = ({\mathsf{D}}_{ij})$ is a $4\times 4$ matrix with entries ${\mathsf{D}}_{32} = \nu$ and ${\mathsf{D}}_{ij} = 0$ otherwise. At high wavenumber, the leading-order component in the linearised momentum equation is given by the new diffusive term itself. Suppose that there is at least one eigenvalue that balances this term. This motivates the following asymptotic expansions for $\sigma$ and $\boldsymbol {q}$ when $k \gg 1$

(3.29a,b)\begin{equation} \sigma = \lambda_2 k^2 + \lambda_1 k + \lambda_0 + \ldots, \quad \boldsymbol{q} = \boldsymbol{q}_0 + \boldsymbol{q}_{{-}1} k^{{-}1} + \boldsymbol{q}_{{-}2} k^{{-}2} + \ldots, \end{equation}

where $\lambda _2$, $\lambda _1$, $\lambda _0$ and $\boldsymbol {q}_0$, $\boldsymbol {q}_{-1}$, $\boldsymbol {q}_{-2}$ are to be determined. At $O(k^2)$, (3.28) reduces to the eigenproblem

(3.30)\begin{equation} \lambda_2 \boldsymbol{\mathsf{A}} \boldsymbol{q}_0 ={-}\boldsymbol{\mathsf{D}}\boldsymbol{q}_0, \end{equation}

with characteristic equation $-\lambda _2^3(\lambda _2 + \nu ) = 0$. When $\lambda _2 = -\nu$, it may be easily verified that $\boldsymbol {q}_0 = \boldsymbol {e}_2$. Therefore, the diffusion operator creates one stable eigenvalue ${\sigma = -\nu k^2 + O(k)}$ associated with viscous damping of $u$. The remaining three solutions are all $\lambda _2 = 0$ and so in these cases $\sigma$ is determined by a lower order balance. Using (3.30), the corresponding vector $\boldsymbol {q}_0$ is determined up to the eigenspace spanned by $\boldsymbol {e}_1$, $\boldsymbol {e}_3$ and $\boldsymbol {e}_4$.

We shall concentrate on the three eigenvalues with $\lambda _2 = 0$, since the mode with $\lambda _2 = -\nu$ is always stable at high $k$. Then, at $O(k)$, (3.28) becomes

(3.31)\begin{equation} (\lambda_1 \boldsymbol{\mathsf{A}} + \mathrm{i}\boldsymbol{\mathsf{B}})\boldsymbol{q}_0 ={-}\boldsymbol{\mathsf{D}} \boldsymbol{q}_{{-}1}. \end{equation}

Note that $\boldsymbol {e}_j^{\textrm {T}} \boldsymbol{\mathsf{D}} = \boldsymbol {0}$ for $j = 1,2,4$, since eddy viscosity appears only in the third row of (3.28). We can therefore eliminate the unknown $\boldsymbol {q}_{-1}$ in (3.31) as so

(3.32)\begin{equation} \boldsymbol{e}_j \boldsymbol{\cdot} (\lambda_1 \boldsymbol{\mathsf{A}} + \mathrm{i}\boldsymbol{\mathsf{B}})\boldsymbol{q}_0 = 0,\quad\text{for}\ j=1,2,4. \end{equation}

Likewise, we can write down the eigenproblem adjoint to (3.31) as $(\lambda _1\boldsymbol{\mathsf{A}} + \mathrm {i}\boldsymbol{\mathsf{B}})^{\textrm {T}}\boldsymbol {r}_0 = -\boldsymbol{\mathsf{D}}^{\textrm {T}}\boldsymbol {r}_{-1}$, where $\boldsymbol {r}_0$ and $\boldsymbol {r}_{-1}$ are unknown left eigenvectors. We will shortly need the constrained problem, with $\boldsymbol {r}_{-1}$ eliminated as so

(3.33)\begin{equation} \boldsymbol{e}_k \boldsymbol{\cdot} (\lambda_1 \boldsymbol{\mathsf{A}} + \mathrm{i}\boldsymbol{\mathsf{B}})^{\textrm{T}}\boldsymbol{r}_0 = 0,\quad\text{for}\ k=1,3,4. \end{equation}

Expanding $\boldsymbol {q}_0 = a_1\boldsymbol {e}_1 + a_3 \boldsymbol {e}_3 + a_4 \boldsymbol {e}_4$ in (3.32) results in a $3\times 3$ eigenvalue problem for the unknowns $a_1$, $a_3$, $a_4$ and $\lambda _1$. Its characteristic equation is

(3.34)\begin{equation} \lambda_1 (\lambda_1 + \mathrm{i})^2 = 0. \end{equation}

So $\lambda _1 = 0$ or $-\mathrm {i}$. In either case, we are forced to proceed further to determine the leading real part of $\sigma$. Therefore, we use the $O(1)$ part of (3.28), which is (for $\lambda _2 = 0$)

(3.35)\begin{equation} (\lambda_0 \boldsymbol{\mathsf{A}} + \boldsymbol{\mathsf{C}})\boldsymbol{q}_0 + (\lambda_1 \boldsymbol{\mathsf{A}} + \mathrm{i} \boldsymbol{\mathsf{B}})\boldsymbol{q}_{{-}1} ={-}\boldsymbol{\mathsf{D}} \boldsymbol{q}_{{-}2}. \end{equation}

We shall divide our pursuit of $\lambda _0$ according to the value of $\lambda _1$.

1. (i) Case: $\lambda _1 = 0$. Solving (3.32) for the eigenvector yields $\boldsymbol {q}_0 = \boldsymbol {e}_4$. To eliminate, the unknown vectors $\boldsymbol {q}_{-1}$ and $\boldsymbol {q}_{-2}$, from (3.35), we simply note that $\boldsymbol {e}_4^{\textrm {T}}\boldsymbol{\mathsf{B}} = \boldsymbol {0}$, since the last row of $\boldsymbol{\mathsf{B}}$ is all zeros (when $Q=0$). Therefore, we project (3.35) onto $\boldsymbol {e}_4$ and substitute in $\boldsymbol {q}_0$ to give

   (3.36)\begin{equation} \lambda_0 ={-}\frac{\boldsymbol{e}_4\boldsymbol{\cdot} \boldsymbol{\mathsf{C}} \boldsymbol{e}_4}{\boldsymbol{e}_4 \boldsymbol{\cdot} \boldsymbol{\mathsf{A}}\boldsymbol{e}_4} = 0. \end{equation}
   Referring back to our expansions (3.29a,b), this means that there always exists an eigenpair $(\sigma , \boldsymbol {q})$ with $\sigma \to 0$ and $\boldsymbol {q} \to \boldsymbol {e}_4$ as $k\to \infty$. Note that this situation corresponds to perturbations of the bedform.

2. (ii) Case: $\lambda _1 = -\mathrm {i}$. Since this is a repeated root, (3.32) only determines $\boldsymbol {q}_0$ within a two-dimensional subspace. Straightforward algebra gives this simply as $\boldsymbol {q}_0 = a_1\boldsymbol {e}_1 + a_3 \boldsymbol {e}_3$. Similarly, the corresponding adjoint eigenproblem (3.33), constrains the left eigenvector $\boldsymbol {r}_0$ to lie in the subspace spanned by $\boldsymbol {e}_1$ and $\boldsymbol {e}_2$. We project (3.35) onto these vectors, yielding

   (3.37a)\begin{gather} \boldsymbol{e}_1 \boldsymbol{\cdot} (\lambda_0 \boldsymbol{\mathsf{A}} + \boldsymbol{\mathsf{C}})\boldsymbol{q}_0 + \mathrm{i} \boldsymbol{e}_2 \boldsymbol{\cdot} \boldsymbol{q}_{{-}1} = 0, \end{gather}
   (3.37b)\begin{gather}\boldsymbol{e}_2 \boldsymbol{\cdot} (\lambda_0 \boldsymbol{\mathsf{A}} + \boldsymbol{\mathsf{C}})\boldsymbol{q}_0 + \mathrm{i} \psi_0\boldsymbol{e}_2 \boldsymbol{\cdot} \boldsymbol{q}_{{-}1} = 0. \end{gather}
   Note that only the second element of $\boldsymbol {q}_{-1}$ appears in these equations. To eliminate it, we return to the full $O(k)$ problem. Since $\boldsymbol {e}_3\boldsymbol {\cdot }\boldsymbol{\mathsf{D}}\boldsymbol {v} = \nu \boldsymbol {e}_2 \boldsymbol {\cdot } \boldsymbol {v}$ for any vector $\boldsymbol {v}$, we project (3.31) onto $\boldsymbol {e}_3$ and rearrange to give $\boldsymbol {e}_2 \boldsymbol {\cdot } \boldsymbol {q}_{-1} = -\mathrm {i} (2a_1 + \Delta \rho a_3) / (2\nu {\textit {Fr}}^2)$. Substituting this expression into (3.37a,b) yields a $2\times 2$ system for $a_1$, $a_3$ and $\lambda _0$, which we solve to find the two eigenvalues
   (3.38)\begin{equation} \lambda_0 = \lambda_\pm({\textit{Fr}}) \equiv \frac{R}{2} + \frac{\pm\sqrt{S} - 1}{2\nu{\textit{Fr}}^2}, \end{equation}
   where $R \equiv \varGamma _{h_0} + \varGamma _{\psi _0}(1-\psi _0)$ and $S \equiv (R\nu {\textit {Fr}}^2 + 1)^2 - 2\nu {\textit {Fr}}^2 \varGamma _{h_0}(\rho _b + 1)$. One of the pair, $\lambda _-$, possesses a singularity at ${\textit {Fr}} = 0$. However, as in § 3.2.2, this is merely an artefact of our choice of dimensionless units that may be removed by an appropriate rescaling. The corresponding eigenvectors are
   (3.39)\begin{equation} \boldsymbol{q}_0 = \boldsymbol{e}_1 + \frac{\nu {\textit{Fr}}^2 (R - 2\varGamma_{h_0}) + 1 \pm \sqrt{S}}{2\varGamma_{\psi_0}\nu{\textit{Fr}}^2 - \Delta\rho}\boldsymbol{e}_3. \end{equation}
   In the limit of vanishing eddy viscosity ($\nu \to 0$), $\boldsymbol {q}_0 \to \boldsymbol {e}_1 - (1\pm 1) \boldsymbol {e}_3 / \Delta \rho$. Therefore, since we anticipate small $\nu$, $\lambda _-$ corresponds largely to growth in $h$ only, whereas $\lambda _+$ corresponds to coupled growth in $h$ and $\psi$. By comparing with the modes of the unregularised problem in (3.19), $\lambda _-$ may be traced to the hydraulic modes and $\lambda _+$ to the third mode related to solid fraction perturbations. The former is responsible for very strong damping, since $\lambda _- \approx - 1/\nu {\textit {Fr}}^2$ for small $\nu$. Conversely, using l'Hôpital's rule, it may further be verified that $\lim _{\nu \to 0} \lambda _+ = \lambda _{0,3}$ from (3.22c). Therefore, in the limit of small $\nu$, this mode is not affected by the regularisation term.

To recap, using (3.29a,b)–(3.39), we have computed the growth rates $\sigma$ of perturbations for non-zero values of $\nu$ at leading order for large wavenumber. These expressions are valid for all ${\textit {Fr}}> 0$ and models of the form (3.28). They are

(3.40a–c)\begin{equation} \sigma ={-}\nu k^2 + O(k), \quad 0 + O(k^{{-}1}), \quad -\mathrm{i} k + \lambda_\pm({\textit{Fr}}) + O(k^{{-}1}), \end{equation}

where $\lambda _\pm ({\textit {Fr}})$ was defined in (3.38). The corresponding mode amplitudes are given by $\boldsymbol {q}_0 = \boldsymbol {e}_2$, $\boldsymbol {e}_4$ and the two expressions in (3.39). Since the growth rates are all bounded above, we conclude that the inclusion of a diffusive term in (3.3c) successfully regularises the singularities in (3.5) that are otherwise present at ${\textit {Fr}} = 1$, removing the problem of ill posedness. However, since the real parts $\lambda _\pm$ of the last two modes remain non-zero, they are still potentially unstable in the $k\gg 1$ regime (if $\lambda _\pm > 0$). We return to this point in § 4.

As $k\to 0$, the diffusive term vanishes in the linearised equations (3.28) to leading order. Therefore, the coefficient $\nu$ sets the effective length scale over which eddy viscosity damps out perturbations. For large-scale geophysical flows where $\nu$ is relatively small, we may thus anticipate linear growth rates for the most part matching those of the $\nu = 0$ problem, with eddy viscosity only affecting very short wavelengths. In this case, the linear stability may still be controlled by the values of the $\nu = 0$ asymptotic growth rates given in (3.22). The intuition here is that for sufficiently small $\nu$ there must be a scale separation between the ‘asymptotic’ ($k\gg 1$) regime of the $\nu = 0$ case and any damping (at still higher $k$) of $\sigma$ induced by turbulent momentum diffusion. We verify this for illustrative model closures in § 4.6.

3.4. Bed load

Returning to the original system (3.3a–d), with the eddy viscosity regularisation $\nu = 0$, we widen our perspective, to allow for non-zero bed flux $Q$. Unfortunately, in this case, analytical solutions of the linear system (3.5) become too complex to work with (even in the long- and short-wavelength regimes) and cease to be useful. Therefore, in this subsection we limit our scope to one important concern: how is the well-posedness of the model affected by $Q$? To address this, we compute the system characteristics, $\lambda _1$, as given by solutions to (3.17), since these are sufficient to determine whether (3.3a–d) may be correctly posed as an initial value problem. Specifically, if the characteristics are all real valued and distinct, the system is strictly hyperbolic and well posed. If instead, any of the characteristics have non-zero imaginary part, the system is not hyperbolic and ill posed (see e.g. Ivrii & Petkov Reference Ivrii and Petkov1974). Alternatively, if the characteristics are all real, but one or more are repeated, the equations are hyperbolic, but may still be ill posed, as we saw in § 3.2.2. In the case of bed load models ($\varGamma = 0$), strict hyperbolicity has previously been demonstrated for various cases. A number of earlier papers derived formulae for the system characteristics by expanding $\lambda _1$ in terms of parameters that are small when the bed dynamics is slow compared with the hydraulic variables (Lyn Reference Lyn1987; Zanré & Needham Reference Zanré and Needham1994; Lyn & Altinakar Reference Lyn and Altinakar2002; Lanzoni et al. Reference Lanzoni, Siviglia, Frascati and Seminara2006). From this perspective, the degeneracy of the system characteristics that underpins ill posedness in the $Q\to 0$ limit is already well appreciated, since it causes naive asymptotic formulae for $\lambda _1$ to break down near ${\textit {Fr}} = 1$ (Lyn Reference Lyn1987; Zanré & Needham Reference Zanré and Needham1994). Later, Cordier et al. (Reference Cordier, Le and Morales de Luna2011) derived general requirements for bed load models to be strictly hyperbolic. We extend their analysis to our setting, which allows for the bulk density variations that may arise with a suspended load. Indeed, since $\varGamma$ does not appear in (3.17), it does not affect the characteristics and so the following analysis applies equally well whether or not bulk entrainment is included.

Bed load terms are most often employed in dilute systems, where we would not expect the basal stresses that drive bed load transport to be sensitive to small changes in the bulk solid fraction. Therefore, we make the additional simplifying assumption that $Q_{\psi _0}$ is small enough that it may be neglected, so $Q_{\psi _0} = 0$ in (3.6b). Then, the characteristic equation resulting from (3.17) reduces to

(3.41)\begin{equation} (\lambda_1 - 1)p(\lambda_1) = 0, \end{equation}

where $p(\lambda _1) = {\textit {Fr}}^2\lambda _1^3 - 2{\textit {Fr}}^2\lambda _1^2 + ({\textit {Fr}}^2 - Q_{u_0} - 1)\lambda _1 + Q_{u_0} - Q_{h_0}$. The system is strictly hyperbolic if and only if (3.41) has four distinct real solutions. Note that these solutions only depend on $Q_{h_0}$, $Q_{u_0}$ and ${\textit {Fr}}$. One of them, arising from the solid mass transport equation (3.3b), is always $\lambda _1 = 1$. In the particular case where $Q_{h_0} = -1$, we also have $p(1)=0$, so this eigenvalue is degenerate. Otherwise, if $Q_{h_0} \neq -1$, then $p(1) = -Q_{h_0} - 1 \neq 0$, and the remaining solutions to (3.41) are never unity. Therefore, we only need to assess the roots of the cubic polynomial $p$ to see if all four characteristics are distinct.

On differentiating $p$ (with respect to $\lambda _1$), its turning points may be found at

(3.42)\begin{equation} \lambda_1 = \frac{2}{3} \pm \frac{1}{3{\textit{Fr}}}\sqrt{{\textit{Fr}}^2 + 3(Q_{u_0} + 1)}. \end{equation}

Hence, a necessary condition for strict hyperbolicity is $Q_{u_0} > -1 -{\textit {Fr}}^2 / 3$. Labelling the two turning points as $\lambda _1 = \ell _\pm$, it follows immediately from considering $p$ as $\lambda _1\to \pm \infty$ in this case, that $\ell _-$ is always a local maximum and $\ell _+$ a local minimum. Therefore, for $p$ to possess three real roots, it must additionally satisfy $p(\ell _-) > 0$ and $p(\ell _+) < 0$. By evaluating $p(\ell _\pm )$ and noting that $\partial p / \partial Q_{h_0} = -1$, it is straightforward to show that $p(\ell _-) > 0$ when $Q_{h_0} < G_+$ and $p(\ell _+) < 0$ when $Q_{h_0} > G_-$, where

(3.43)\begin{equation} G_\pm({\textit{Fr}}, Q_{u_0}) = \frac{1}{27{\textit{Fr}}}[2{\textit{Fr}}^3 \pm 2({\textit{Fr}}^2 + 3Q_{u_0} + 3)^{3/2} + 9 {\textit{Fr}}(Q_{u_0} - 2)]. \end{equation}

Therefore, the region of parameter space where the system is strictly hyperbolic satisfies $G_- < Q_{h_0} < G_+$. Conversely, when either $Q_{h_0} < G_-$ or $Q_{h_0} > G_+$ two of the characteristics, i.e., roots of (3.41), are complex conjugate and the system is non-hyperbolic. We now seek to identify constraints on $Q_{h_0}$ and $Q_{u_0}$ such that the system is strictly hyperbolic for all ${\textit {Fr}}$.

In the particular case when $Q_{h_0} = Q_{u_0}$, the solutions of (3.41) may be readily computed to be

(3.44)\begin{equation} \lambda_1 = 1 \pm {\textit{Fr}}^{{-}1} \sqrt{Q_{h_0} + 1}, 1, 0. \end{equation}

These expressions are commensurate with the case $Q\to 0$, whose characteristics were given in (3.18). Here, only the hydraulic modes are altered by the bed load term. All four values are distinct (so $G_- < Q_{h_0} < G_+$), unless ${\textit {Fr}} = \sqrt {Q_{h_0} + 1}$, where a hydraulic mode intersects with the bed characteristic $\lambda _1 = 0$, and $Q_{h_0} = Q_{u_0} = G_+$. On differentiating (3.43) with respect to ${\textit {Fr}}$, it can be shown that this point is a global minimum of $G_+({\textit {Fr}}, Q_{u_0})$ for any fixed $Q_{u_0}$. Therefore, $Q_{h_0} < G_+$ for all ${\textit {Fr}}$ only if $Q_{h_0} < Q_{u_0}$.

The lower limit $G_-$, is a strictly increasing function of ${\textit {Fr}}$ (for any $Q_{u_0}$), with $G_-({\textit {Fr}}, Q_{u_0}) \to -\infty$ as ${\textit {Fr}} \to 0$ and $G_-({\textit {Fr}}, Q_{u_0}) \to -1$ as ${\textit {Fr}} \to \infty$. Combining this information with the lower bound for $G_+$, we conclude that the system (with $Q_{\psi _0}=0$ assumed) is strictly hyperbolic over all Froude numbers if and only if

(3.45)\begin{equation} {-}1 < Q_{h_0} < Q_{u_0}. \end{equation}

We note, with reference to (3.42), that this stronger condition automatically satisfies the requirement that $p$ has two turning points. In figure 2(a), we plot examples of the bounds $G_\pm$, as (dotted) curves in $({\textit {Fr}},Q_{h_0})$-space for fixed $Q_{u_0}$, indicating the regions where the model fails to be hyperbolic. The axes have been chosen so as to encompass very general $Q$ closures. However, fortunately in applications sediment flux typically depends only very weakly, if at all on the flow depth, so $|Q_{h_0}| \ll 1$ is expected. In fact, it is common in fluvial models to have $Q_{u_0}$ strictly positive and $Q_{h_0} = 0$ or $-1 \ll Q_{h_0} < 0$ (e.g. in the latter case, if a Manning friction law is employed). Therefore, most studies that include bed load operate in a regime where (3.45) is satisfied.

Figure 2. Hyperbolicity of the morphodynamic model equations depends on the bed load function $Q$. (a) Regions of non-hyperbolicity as a function of ${\textit {Fr}}$ and $Q_{h_0}$, for fixed $Q_{u_0} = 0.1$ (purple shading), $1$ (blue shading) and $Q_{\psi _0} = 0$. Outside these regions, the model is strictly hyperbolic, save along the bounding curves $G_\pm$ (dotted lines) and the special case $Q_{h_0} = -1$ (dashed line), where one of the roots of $c$ intersects with the solid mass transport characteristic. (b) System characteristics as a function of ${\textit {Fr}}$ for $Q_{u_0} = 0.1$, $Q_{\psi _0} = 0$ and $Q_{h_0} = 0.095$ (dotted lines), $0.105$ (solid lines). We label the curves I–IV according to the ordering of the corresponding characteristics derived in (3.18) for the $Q=0$ case.

This analysis suggests an alternative to the regularisation strategy of § 3.3, since adding even a small bed load flux term can ensure that the model equations are well posed, provided that (3.45) is satisfied. We visualise the effect of bed load on the system characteristics in figure 2(b), either side of the threshold $Q_{h_0} = Q_{u_0}$ where the system loses hyperbolicity. We plot $\mathrm {Re}(\lambda _1)$ as a function of ${\textit {Fr}}$ for $Q_{u_0} = 0.1$ and $Q_{h_0} = 0.095$ (dotted lines), $0.105$ (solid lines). The four branches of $\lambda _1$ are labelled I–IV, according to the ordering of the corresponding modes in the $Q = 0$ case, adopted in § 3.2.2. When $Q_{h_0} = Q_{u_0} = 0.1$ (not shown), the mode II and IV characteristics intersect at a single point, ${\textit {Fr}} = \sqrt {1.1} \approx 1.05$ in this case, which is easily calculated using the expressions in (3.44). Decreasing $Q_{u_0}$ (so that $Q_{u_0} < Q_{h_0}$) causes these characteristics to coalesce into a complex conjugate pair, resulting in a region ($0.95 \lesssim {\textit {Fr}} \lesssim 1.15$) where the system is non-hyperbolic. Conversely, increasing $Q_{u_0}$ separates the intersecting characteristics so that four real values are present across all Froude numbers. (Note that these separated curves are labelled with both II and IV, since they originate from different modes at either end.) The other two characteristics (I and III) are essentially unaffected by small changes in the bed load.

4.1. Model closures

In order to capture a range of different sedimentary flows, from dilute suspensions to fully granular flow, we use a mixed drag formulation that depends on the bulk solid fraction, writing

(4.1)\begin{equation} \tilde{\tau} = (1 - \psi) \tilde{\tau}_f + \psi \tilde{\tau}_g, \end{equation}

where $\tilde {\tau }_f$ and $\tilde {\tau }_g$ are fluid and granular drag laws respectively. We assume that the bed is a saturated mixture of fluid and sediment containing the maximum possible sediment concentration $\tilde {\psi }_b = \tilde {\psi }^*$. The maximum solid fraction $\tilde {\psi }^*$ depends on how efficiently particles can be packed and is typically observed to be around 60–70 % (Santiso & Müller Reference Santiso and Müller2002; Farr & Groot Reference Farr and Groot2009). Since $\psi = \tilde {\psi } / \tilde {\psi }_b$, we have $0 \leq \psi \leq 1$ for all flows and therefore (4.1) contains all weighted combinations of fluid and granular drag. For the fluid law, we employ the common Chézy formula for turbulent shear stress, $\tilde {\tau }_f = C_d \tilde {\rho } \tilde {u}^2$, where $C_d$ is a drag coefficient, which we assume to be constant. For the granular drag, we use the frictional law due to Pouliquen & Forterre (Reference Pouliquen and Forterre2002), which sets $\tilde {\tau }_g = \mu (I) \tilde {\rho } g \cos (\phi ) \tilde {h}$. The phenomenological parameter $\mu$ is modelled as an increasing function of the dimensionless inertial number $I \equiv uh^{-3/2}d{\textit {Fr}}$, where $d = \tilde {d} / \tilde {h}_0$ and $\tilde {d}$ denotes the characteristic diameter of the sediment particles. It is constructed so as to vary smoothly between a lower (static) limit $\mu _1$ and an upper (dynamic) bound $\mu _2$, with $0 < \mu _1 < \mu _2$, and takes the form

(4.2)\begin{equation} \mu(I) = \mu_1 + \frac{\mu_2 - \mu_1}{1 + \beta I^{{-}1}}, \end{equation}

where $\beta$ is a positive constant that may be determined experimentally.

We suppose that mass transfer is governed by the competing processes of bed erosion at a rate $\tilde {E}$ and particle deposition at rate $\tilde {D}$, writing $\tilde {\varGamma } = \tilde {E} - \tilde {D}$. Since these processes take place at the scale of individual particles, we opt to non-dimensionalise these closure terms using the velocity $\tilde {u}_p = (g'_\perp \tilde {d})^{1/2}$, where $g_\perp ' = g \cos \phi ( \tilde {\rho }_s / \tilde {\rho }_f - 1 )$ is the reduced gravity for a particle in dilute suspension, resolved perpendicular to the slope. The dimensionless transfer rates are then $E_p = \tilde {E} / \tilde {u}_p$ and $D_p = \tilde {D} / \tilde {u}_p$. This rescaling allows us to fix appropriate dimensionless constants for these closures when considering the steady balance $E_p = D_p$ independently. However, note that care must be taken when reintroducing these terms into (3.3a–d), which uses a different velocity scale for $\partial b / \partial t$. Specifically, if $\varGamma _p \equiv E_p - D_p$, then (3.2h) implies that $\varGamma = \varGamma _p {\textit {Fr}} d^{1/2}(\rho _s/\rho _f-1)^{1/2}\cot \phi$.

Entrainment of particles into the flow is caused by turbulent shear stresses at the bed, which must overcome the static friction experienced by resting grains. Competition between $\sim \tilde {\tau } \tilde {d}^2$ drag forces and $\sim \tilde {\rho }_f g_\perp ' \tilde {d}^3$ frictional forces (assumed proportional to the submerged weight of individual grains) can be captured by their ratio, the dimensionless Shields number $\theta \equiv \tilde {\tau } / (\tilde {\rho }_f g_\perp ' \tilde {d})$. Experimental observations for dilute flows suggest that at sufficiently high drag, flow erosion obeys a power law of the form $\tilde {E} \propto (\theta - \theta _c)^{3/2}$, where $\theta _c$ is a critical Shields number below which there is no entrainment. Beyond this there is considerable disagreement concerning both the exact functional form for $\tilde {E}$ and its magnitude (Lajeunesse, Malverti & Charru Reference Lajeunesse, Malverti and Charru2010). Moreover, it is unclear whether this general erosion model applies for more concentrated suspensions, where effects such as the pore water pressure modify the force relationship encapsulated in the Shields number. Since our aim here is only to elucidate some general properties of solutions, we prefer simplicity here and suppose that $\tilde {E}$ depends linearly on $\tilde {u}_p (\theta - \theta _c)^{3/2}$. However, we shall make one important modification to this dilute erosion law. Since concentrated layers may be held static on shallow grades by their granular friction, we must not permit erosion to occur in situations where $\theta > \theta _c$, yet $\tilde {u} = 0$. Therefore, we set $\theta _c = \theta _c^* + \theta _0$, where $\theta _c^*$ denotes the usual critical Shields number (for dilute suspensions) and $\theta _0(\tilde {h}, \tilde {\psi }) = \theta |_{\tilde {u} = 0}$, i.e., the Shields number of a resting flow, which may become large as $\tilde {\psi }$ increases. For simplicity, we consider $\theta _c^*$ constant in this study, even though in principle it depends on flow properties such as the particle Reynolds number $Re_p \equiv \tilde {u}_p \tilde {d}/\tilde {\nu }_f$, where $\tilde {\nu }_f$ is the kinematic viscosity of the fluid (Soulsby Reference Soulsby1997). On dividing through by $\tilde {u}_p$, the dimensionless entrainment rate is then

(4.3)\begin{equation} E_p(h, u, \psi) = \begin{cases} \varepsilon \left[\theta(h,u,\psi) - \theta_c(h, \psi)\right]^{3/2}, & \text{if } \theta > \theta_c, \\ 0, & \text{otherwise}, \end{cases} \end{equation}

where $\varepsilon$ is a proportionality coefficient that characterises the erodibility of the bed.

We treat sediment deposition as being governed by a process of hindered settling. At low concentrations, particles settle independently, so the (monodisperse) sediment deposits at a rate ${\sim }\tilde {w}_s \psi$, where $\tilde {w}_s$ denotes the characteristic falling speed of the grains. As concentrations increase, pure settling becomes disrupted as particles increasingly interact, ultimately shutting off as $\psi \to 1$ and grains can no longer fall (in a time-averaged sense) under gravity. Therefore, we take the deposition term to be

(4.4)\begin{equation} D_p(\psi) = w_s \psi (1 - \psi), \end{equation}

where $w_s = \tilde {w}_s / \tilde {u}_p$. This a slight simplification of the widely used formula due to Richardson & Zaki (Reference Richardson and Zaki1954). More detailed and accurate expressions for $D_p$ typically involve empirical fits featuring the same essential form (e.g. Spearman & Manning Reference Spearman and Manning2017).

When considering a bed load, we use the following standard expression, which mirrors the entrainment rate of (4.3)

(4.5)\begin{equation} Q_p = \begin{cases} \gamma \left[\theta(h,u,\psi) - \theta_c(h, \psi)\right]^{3/2}, & \text{if } \theta > \theta_c, \\ 0, & \text{otherwise}, \end{cases} \end{equation}

where $Q_p$ is a particle-scale non-dimensionalisation such that $Q_p = \tilde {Q} / (\tilde {d} \tilde {u}_p)$ and $\gamma$ is a constant that dictates the flux strength. For example, $\gamma = 8$ sets the well-known Meyer-Peter & Müller formula (Meyer-Peter & Müller Reference Meyer-Peter and Müller1948). The corresponding flow-scale non-dimensionalisation, as used in (3.3d) and given by (3.2i), is $Q = Q_p \tilde {d} \tilde {u}_p / (\tilde {h}_0 \tilde {u}_0) = Q_p d^{3/2} (\rho _s / \rho _f - 1)^{1/2} / {\textit {Fr}}$.

Finally, the basal velocity closure dictates the characteristic downslope flow speed at the bed, during particle entrainment. Where needed, we assume that it can be modelled by a turbulent friction velocity of the form

(4.6)\begin{equation} \tilde{u}_b = \sqrt{\tilde{\tau} / \tilde{\rho}}. \end{equation}

A large number of free parameters are involved in specifying these various model closures. Therefore, we choose to fix some illustrative values, as listed in table 1, making it clear whenever we deviate from these defaults. An investigation of other parameter choices indicated that our observations below are qualitatively robust to variations in these values.

Table 1. Illustrative model parameters. Except where otherwise stated, all results in § 4 assume these values. Definitions for these parameters may be found in (2.2), (2.3b), (4.1)–(4.4) and the accompanying text in each case.

4.2. Existence of steady layers

We are now in a position to assess when uniform flowing layers can exist in equilibrium. This is dictated by the steady balances in (3.1a,b), which we recall enforce that drag balances the downslope component of gravitational forcing and that there is no net material transfer between the flow and the bed. Note that these conditions are independent of whether there is a bed load or not. To begin with, we concentrate on the stress balance and assume that mass transfer with the bulk is negligible ($\varGamma \to 0$), thereby automatically satisfying (3.1b). Substituting our mixed drag law (4.1) into (3.1a), non-dimensionalising and rearranging gives the following condition for existence of a steady layer of solid fraction $\psi _0$:

(4.7)\begin{equation} \tan\phi - (1-\psi_0)C_d{\textit{Fr}}^2 = \psi_0 \mu(d{\textit{Fr}}). \end{equation}

In the dilute limit, $\psi _0 \to 0$, where drag is purely fluid like, this equation selects a unique Froude number, ${\textit {Fr}} = \sqrt {\tan \phi /C_d}$. Such states become unstable when ${\textit {Fr}} > 2$ (Jeffreys Reference Jeffreys1925). Conversely, in the concentrated limit, $\psi _0 \to 1$, where drag is purely granular, steady layers adopt a unique inertial number, given by $I_0 \equiv d{\textit {Fr}} = \mu ^{-1}(\tan \phi )$. Since $\mu _1 < \mu (I) < \mu _2$, only a range of slope angles (between $\arctan \mu _1$ and $\arctan \mu _2$) are permitted. The threshold for linear instability in this case does not depend on $\mu$ and was computed by Forterre & Pouliquen (Reference Forterre and Pouliquen2003) to be ${\textit {Fr}} > 2/3$.

For intermediate values of $\psi _0$, since the left-hand hand side in (4.7) is a decreasing function of ${\textit {Fr}}$ and unbounded below, the steady drag balance may be satisfied as long as $\tan \phi \geq \psi _0 \mu _1$. The effect of the Chézy drag term thereby relaxes the limits on existence imposed by the granular law. Steady layers that are more dilute can exist in mobile equilibrium at shallower slope angles, i.e., down to $\arctan (\psi _0 \mu _1)$, while steep steady flows ($\tan \phi \to \infty$) may always be maintained at a suitably high ${\textit {Fr}}$, since the turbulent drag, $C_d{\textit {Fr}}^2$, is not bounded above. However, note that such solutions are not necessarily stable. Indeed, the stability threshold for these flows may be readily computed using Trowbridge's criterion (3.10). After non-dimensionalising (4.1) and differentiating, one sees that $\tau _{u_0} = [2(1-\psi _0)C_d{\textit {Fr}}^2 + \psi _0\mu ' I_{u_0}]\cot \phi$, where $\mu ' = \partial \mu / \partial I |_{I=I_0}$ and $I_{u_0} = \partial I / \partial u |_{h,u=1} = d{\textit {Fr}}$. Likewise, $\tau _{h_0} = \psi _0[\mu (d{\textit {Fr}}) + \mu '(d{\textit {Fr}})I_{h_0}]\cot \phi$, with $I_{h_0} = \partial I / \partial h |_{h,u=1} = -\frac {3}{2} d{\textit {Fr}}$. On substituting these expressions into (3.10) and rearranging, one sees that these states are unstable when

(4.8)\begin{equation} (1-\psi_0)C_d ({\textit{Fr}} - 2) + \psi_0 d \mu' (3/2 - 1 / {\textit{Fr}}) > 0. \end{equation}

Note that this criterion smoothly interpolates between the ${\textit {Fr}} = 2$ and ${\textit {Fr}} = 2/3$ thresholds for the special cases of purely fluid ($\psi _0 = 0$) and granular ($\psi _0 = 1$) flows.

We summarise the existence of non-erosive layers subject to the drag law (4.1) in figure 3. Dashed contours trace out the unidimensional family of steady layers for each slope angle, across $(d, {\textit {Fr}})$-parameter space, computed from (4.7), with the three panes corresponding to different fixed solid fractions, $\psi _0 = 0.2$, $0.5$ and $0.8$, from left to right. The minimum slope angles for steady flows (dashed red contours) follow the line ${\textit {Fr}} = 0$. Red and blue filled contours show the asymptotic growing and decaying growth rates of perturbations respectively, computed by substituting $\tau _{u_0}$ and $\tau _{h_0}$ from above into the limiting formula given earlier in (3.9). These are separated by the neutral stability boundary (4.8), which is displayed in solid black. When $d \ll 1$ (small particles, relative to the flow depth), the drag is dominated by the Chézy term. In this regime, which covers most physically realisable grain sizes, growth rates are essentially independent of $d$ and the stability boundary is ${\textit {Fr}} \approx 2$. Increasing $d$ outside this region leads to less stable flows and lowers the stability boundary. Increasing the solid fraction generally leads to less severe growth rates, but decreases the range of slope angles that permit stable steady flows (through increasing the minimum slope angle). We find the qualitative properties of this plot to be largely insensitive to our specific choices of $C_d$, $\mu _1$, $\mu _2$ and $\beta$ whose values were given in table 1.

Figure 3. Existence of steady layers for our mixed drag formulation (4.1), without morphodynamics. Dashed lines show the existence of steady flows at fixed slope angles as labelled, with the lowest red dashed line indicating the minimum slope angle $\phi = \arctan (\psi _0\mu _1)$. Filled contours show maximum values of the non-dimensional linear growth rates given in (3.8). Stable (blue) and unstable (red) regions are separated by the neutral stability boundary (4.8) plotted in solid black. Successive panes represent increasing solid fractions $\psi _0$ from left to right as follows: ($a$) 0.2, ($b$) 0.5 and ($c$) 0.8.

When the morphodynamics is non-negligible, steady flows must also satisfy (3.1b). That is, erosion and deposition must be everywhere in balance, $E_p(h,u,\psi ) = D_p(\psi )$. This condition dictates the solid fractions where mass transfer can be in equilibrium. Despite our efforts to obtain simple closures in (4.3) and (4.4), it is a complex nonlinear equation that depends on many physical parameters. Nevertheless, we can determine some generic properties of solutions.

We first consider the system parameters to be fixed (but arbitrary) and suppose that the flow is in uniform steady balance with its drag (so $h = u = 1$), leaving $D_p$ and $E_p$ functions of $\psi$ only. We also assume that $\theta > \theta _c$, so that there is some particle entrainment. Then, in particular, $E_p(0) > 0$ and $E_p(1) > 0$. Conversely, the deposition rate curve (4.4) always obeys $D_p(0) = D_p(1) = 0$. Therefore, since $E_p - D_p > 0$ for both $\psi = 0$ and $1$, either: (a) there exist an even number of coexistent steady flows with different solid fractions, (b) erosion balances deposition exactly at a turning point in $E_p-D_p$ or (c) erosion always exceeds deposition. We visualise these three cases in figure 4, where we plot $D_p(\psi )$ and $E_p(\psi )$ at different values of ${\textit {Fr}}$, using the illustrative model parameters of table 1. Aside from where explicitly stated otherwise, these parameters are fixed for the remainder of this section. Note that our choice of $D_p$ does not depend on ${\textit {Fr}}$, while the Froude number dependence of $E_p$ enters through the basal drag term in the Shields number.

Figure 4. Example deposition and erosion closures, as functions of the solid fraction $\psi$. Steady flows occur where $D_p$ and $E_p$ intersect. There are three cases: ($a$) two steady states, one dilute and one concentrated (${\textit {Fr}} = 0.75$); ($b$) a single steady state where $E_p$ is tangent to $D_p$ (${\textit {Fr}} \approx 2.63$); ($c$) erosion always exceeds deposition (no steady states) (${\textit {Fr}} = 3.25$).

Case ($a$), where there are multiple steady flows, occurs at lower ${\textit {Fr}}$ numbers. Here, erosion increases with solids concentration, intersecting the deposition curve at two points. This leads to two corresponding steady flows: one dilute and one relatively concentrated. Additional intersections (leading to three or more steady flows) are a possibility, but would require a very particular erosion curve. Physical intuition suggests that the dilute solution is stable, since perturbations in either direction cause negative feedback: decreasing $\psi$ away from this state leads to $E_p > D_p$, while increasing $\psi$ leads to $E_p < D_p$. Likewise, the concentrated solution (at $\psi \approx 0.94$ in figure 4) invites either runaway deposition (if $\psi$ decreases) or runaway erosion (if $\psi$ increases). This process suggests a possible mechanism underlying sediment distribution in debris flows, which commonly feature an unsteady highly concentrated front trailed by a steady dilute layer (Pierson Reference Pierson1986; Hungr Reference Hungr2000; Ancey Reference Ancey2001). As ${\textit {Fr}}$ increases in figure 4, the pair of steady flows coalesces at a single point; this is case ($b$). Beyond this point, no steady solutions exist, case ($c$). Here, erosion everywhere exceeds deposition. Uniform layers in this case can never be truly steady, since they can only satisfy one of (3.1a,b). If the drag is ever in equilibrium with gravitational forces, then net entrainment injects material into the flow.

4.3. Global morphodynamic modes

The instability mechanism identified in figure 4 is purely morphodynamic and depends on a straightforward criterion: states are unstable to this mode if $\varGamma _{\psi _0} > 0$. The process is fundamentally an instability to uniform perturbations in flow concentration, though since there is no intrinsic dependence upon spatial gradients we might anticipate that it manifests as a destabilising feature at all wavenumbers. However, this picture is a simplification, since it omits feedbacks from the other flow fields. The full situation for uniform disturbances is contained within our analysis of zero-wavenumber perturbations in § 3.2.1, where we computed the general growth rates of the four different linear stability modes for $k=0$, two of which are always neutrally stable. As discussed earlier, if the approximation of small $\varGamma _{u_0}$ can be made, the two remaining modes may be understood simply: one is inherited from the hydrodynamic stability problem, the other contains morphodynamic feedbacks. Their growth rates are given in (3.15a,b). The morphodynamic rate in (3.15b) may by understood as a competition between two mechanisms for growth in $h$ and $\psi$. The process for $\psi$ (when $\varGamma _{\psi _0} > 0$) has already been outlined. If $\varGamma _{h_0} > 0$, then a small increase in the steady layer height leads to net entrainment which, via (3.3a), enhances growth in $h$ in turn. Likewise, a small decrease in $h$ would cause the depth to decrease away from its steady value. Conversely, if $\varGamma _{h_0} < 0$, then this term is stabilising. If $\varGamma _{h_0}$ is negligible with respect to $\varGamma _{\psi _0}(1-\psi _0)$, then the morphodynamic mode has approximate growth rate $\sigma ^*_m$, defined by

(4.9)\begin{equation} \sigma_m^* = \varGamma_{\psi_0}(1- \psi_0). \end{equation}

In this case, stability is only governed by the mechanism for concentration growth encapsulated by figure 4.

The accuracy of the above physical picture depends on the reliability of the approximations made in reaching (3.15a,b) and (4.9). These estimates are plausible (especially at lower $\psi _0$), since we might expect the relative steepness of the hindered settling curve (plotted in figure 4) to be more significant than gradients in $E_p$, which are solely responsible for dictating $\varGamma _{h_0}$, $\varGamma _{u_0}$ and depend on the small parameters $\varepsilon$ and $C_d$.

For a more detailed analysis, in figures 5(a) and 5(b) we show both non-zero branches of the exact growth rates $\sigma$ (solid lines) for global disturbances (which are purely real), as given in (3.13a)–(3.14b), for (a) dilute and (b) concentrated states. On the same axes, we plot both approximations to the rates: the more general formulae from (3.15a,b), which we denote $\sigma _a$ (dashed lines), and the cruder approximation to the morphodynamic mode rate $\sigma _m^*$ (yellow circles), made above in (4.9). For the dilute states, we confirm that $\sigma < 0$ across the range where steady layers exist ($0.22 \lesssim {\textit {Fr}} \lesssim 2.62$). Moreover, both branches of $\sigma$ are well approximated by $\sigma _a$ for ${\textit {Fr}} \lesssim 1.5$, and $\sigma _m^*$ lies very close to its corresponding branch of $\sigma _a$, indicating that $\varGamma _{\psi _0}$ is indeed primarily responsible for setting the sign of $\sigma$ in this case. At higher ${\textit {Fr}}$, the approximating curves are less accurate. However, this is to be expected. Whereas at low ${\textit {Fr}}$, the dilute solutions have $\psi _0 \approx 0$, which is where the hindered settling curve is at its steepest (and consequently where $|\varGamma _{\psi _0}|$ can be considered large compared with other derivatives of $\varGamma$), for higher ${\textit {Fr}}$ states have higher $\psi _0$ and the approximations of negligible $|\varGamma _{h_0}|$ and $|\varGamma _{u_0}|$ gradually break down as $\psi _0$ approaches the turning point in $D_p$ (see figure 4). However, we note that throughout, $\sigma _m^*$ lies close to $\sigma _a$ since $|\varGamma _{h_0}|$ is small for dilute states. Furthermore, we need only be strictly concerned with ${\textit {Fr}} \lesssim 1$, since outside this regime layers are susceptible to high-wavenumber instabilities.

Figure 5. Growth rates for uniform ($k=0$) perturbations as a function of ${\textit {Fr}}$, for the (a) dilute and (b) concentrated steady solution families identified in figure 4. Solid curves show the two non-zero branches of the exact growth rate $\sigma$, computed from the formulae in (3.13b)–(3.14b). The signs given in the legend indicate the branch, according to the sign of $\pm \sqrt {s_c}$ in (3.13b). Also plotted are two approximations to the exact rates, $\sigma _a$ (dashed lines), as defined in (3.15), and $\sigma _m^*$ (yellow circles), as defined in (4.9).

By contrast, the approximations to the growth rates for the concentrated solutions, in figure 5(b), are not especially good, as highlighted by the figure inset. Therefore, $\varGamma _{u_0}$ cannot be neglected here. Unfortunately, while perturbations in $h$ and $\psi$ and their respective feedbacks may be understood in simple physical terms, this is not easy to do in general for $u$, due to the many interacting contributions to momentum present in the governing equations. However, in this particular case, we note that the intuition that concentrated steady states are typically unstable is borne out, meaning that the feedbacks omitted in making the approximation $\sigma _a$ are not stabilising on aggregate. In fact, this is guaranteed, since the pair of steady states in figure 4 arises in a saddle-node bifurcation as ${\textit {Fr}}$ is decreased from infinity. This implies that, since the dilute flow is stable, the concentrated flow must have at least one unstable direction.

4.4. Linear growth rates for general wavenumbers

Instability to uniform disturbances is not the only feature introduced by the presence of morphodynamics, as our earlier analysis in § 3.2 indicates, since states are also vulnerable to the high-wavenumber instability near unit Froude number. Therefore, we broaden the discussion to incorporate the full linear stability problem. We begin by numerically solving the eigenproblem in (3.5) using our chosen model closures, over a range of finite wavenumbers. Each of the four eigenmodes in the problem possesses a linear growth rate continuously parametrised by $k$. As in § 3.2.2, we label the modes I–IV according to the ordering of their asymptotes used in (3.22a–d). Recall that in the high-$k$ regime, modes I and II are analogues to the modes of the purely hydraulic stability problem, whereas III and IV are additional morphodynamic modes that, involve perturbations in the solid fraction and bed surface respectively.

We denote the growth rates of each mode, indexed by wavenumber as $\sigma _n(k)$ for $n = 1,\ldots ,4$. In figure 6(a) we plot $\max _n \mathrm {Re}(\sigma _{n})$ versus $k$, for states on the dilute solution branch. Additionally, we plot their limiting high-$k$ growth rates, taken from the maxima of the four expressions derived in (3.22a–d), in dotted grey. Each curve passes through the origin, since the maximum growth at $k=0$ is given by the neutral modes for these states. The ${\textit {Fr}} = 0.5$ and $0.75$ solutions are stable to all perturbations, just as they would be in the non-erosive case, with the latter solution being more strongly damped. The ${\textit {Fr}} = 1.25$ curve is stable to long-wavelength disturbances and becomes unstable at $k \approx 6$. Its maximum over all $k$ is given by its asymptotic value, to which it converges more slowly than the other curves. This solution would be stable in the non-erosive case: given the same solid fraction ($\psi _0 \approx 0.02$), according to (4.8), non-erosive states turn unstable at ${\textit {Fr}} \approx 1.98$. The ${\textit {Fr}} = 2$ state is stable for a narrower range of wavelengths, becoming unstable at $k \approx 3$. This state (which has $\psi _0 \approx 0.1$) would also be unstable in the non-erosive situation. Its asymptotic growth rate is lower than the ${\textit {Fr}} = 1.25$ curve, which we shall see shortly is because it lies further from the ${\textit {Fr}} = 1$ singularities present in (3.22b,d). Aside from a narrow interval ($k \lesssim 1.5$) at small $k$ in the ${\textit {Fr}} = 2$ case where mode I dominates (not visible at this scale), each of these maximum growth rates for dilute states are everywhere given by the growth of the morphodynamic mode IV.

Figure 6. Perturbation growth rates for morphodynamic states as functions of wavenumber in the unregularised problem ($\nu = 0$). The curves were computed by taking the maximum real part of the four normal mode growth rates arising from (3.5), across a range of wavenumbers $k$, for states on the (a) dilute and (b) concentrated solution branches identified in figure 4. The Froude numbers are ${\textit {Fr}} = 0.5$ (orange), $0.75$ (blue), $1.25$ (olive), $2$ (teal). High-$k$ asymptotes, computed from maxima of the four expressions (3.22a–d), are plotted in dotted grey.

In figure 6(b), we plot the corresponding maximum growth rate curves for the concentrated solution branch. These match the intuition from § 4.3, that concentrated states should be everywhere unstable (since we expect the basic global instability mechanism to persist regardless of $k$). Growth at $k=0$ is a local maximum for each ${\textit {Fr}}$ and the corresponding instability is due to mode III, which dominates over all $k$ for ${\textit {Fr}} = 0.5$, $0.75$ and ${\textit {Fr}} = 2$. (The two lower ${\textit {Fr}}$ growth rate curves turn sharply at $k \approx 0.025$ and $0.075$ respectively, but are nevertheless smooth, as indicated on the inset.) Conversely, the ${\textit {Fr}} = 1.25$ curve is formed by a crossing of growth rates for modes III and IV, the latter of which is neutral at $k = 0$. The crossing point (at $k \approx 1.6$) is shown in an inset, with the subdominant portions of the mode III and IV curves also included (dashed lines).

The variations of the modes as functions of Froude number are encapsulated in figure 7. Here, we plot the asymptotic ($k \gg 1$) growth rates $\lambda _{0,i}$, given in (3.22a–d), for $i=1,\ldots 4$, with dashed lines. Their maximum for each ${\textit {Fr}}$ is overlaid as a solid line. With filled circles, we plot the maximum growth rate over all $k$. Therefore, any discrepancy between the solid lines and circles indicates that maximal growth is attained at some finite wavenumber. Figure 7(a) presents the data for the dilute solutions. At low ${\textit {Fr}}$, solutions are stable and the maximum possible growth is zero, due to the ($k=0$) neutral modes. The asymptotic growth rate is dictated by the mode IV curve, which is briefly surpassed by mode I at ${\textit {Fr}} \approx 0.85$, before growth is dominated by the singular behaviour of modes II (for ${\textit {Fr}} < 1$) and IV (for ${\textit {Fr}} > 1$), which diverge at ${\textit {Fr}} = 1$ with opposite sign. This induces instability at ${\textit {Fr}} \approx 0.9$. For all higher Froude numbers, the solutions remain unstable, even though they would be sufficiently dilute to remain stable well past ${\textit {Fr}} = 1$ if morphodynamic effects were neglected. We also note that the asymptotic rate correctly identifies the onset of instability and matches the maximum rate thereafter, thereby justifying the focus on short wavelengths in our analysis.

Figure 7. The maximum of the high-wavenumber asymptotic growth rates, $\max _i\lambda _{0,i}$, plotted with solid curves as a function of ${\textit {Fr}}$, for the (a) dilute and (b) concentrated solution branches. With dashed lines we plot the individual $\lambda _{0,i}({\textit {Fr}})$ curves, as labelled for $i=1,\ldots ,4$. Also plotted (filled circles) are the maximum growth rates of the corresponding solutions, over all wavenumbers, i.e. $\max _{k,n}\mathrm {Re}[\sigma _n(k)]$.

The corresponding curves for the concentrated solution branch are shown in figure 7(b). As expected, this also features a singularity at ${\textit {Fr}} = 1$ and is everywhere unstable. The maximum growth rates (filled circles) exactly match the corresponding asymptotic limits, except at low Froude number (${\textit {Fr}} \lesssim 0.85$), where growth at $k=0$ is slightly larger than in the $k\gg 1$ regime, as we saw in figure 6(b). Mode III is always unstable and dominates the other modes over a large region. Near the singularity it is surpassed by modes II (for ${\textit {Fr}} < 1$) and IV (for ${\textit {Fr}} > 1$) and it is briefly surpassed again by mode IV near ${\textit {Fr}} \gtrsim 2.6$, where the dilute and concentrated solution branches coalesce.

4.5. Effect of bed erodibility

In addition to the singularity introduced at unit Froude number, a further effect of the morphodynamics on dilute states may by identified in figure 7(a). In the limit of weak sediment entrainment, $\varepsilon \to 0$, mode I must turn unstable at ${\textit {Fr}} = 2$, due to the classical roll-wave instability for Chézy bottom drag (Jeffreys Reference Jeffreys1925). Conversely, in the regime of figure 7 ($\varepsilon =0.01$), this mode remains stable throughout the range of ${\textit {Fr}}$ where steady states exist ($0.22\lesssim {\textit {Fr}} \lesssim 2.62$). Unstable growth only occurs for mode IV.

We observe the effect of increasing $\varepsilon$ from a weakly erodible regime by reproducing the figure 7(a) plot for dilute states at various values of $\varepsilon$, in figure 8. The other model parameters remain as stated in table 1. Beginning with the case of small $\varepsilon = 3\times 10^{-4}$, in figure 8(a), the signature of hydraulic roll-wave instability is clear. Mode I becomes unstable at a Froude number just above $2$ and immediately dominates over the morphodynamic mode. Outside the singular point ${\textit {Fr}} = 1$, the high-$k$ growth rate of the latter mode (IV) is $O(\varepsilon )$, which may be confirmed by consulting the formula derived in (3.22d) and the closure for entrainment in (4.3). Consequently, there is a narrowing of the effective influence of the singularity, relative to the picture in figure 7(a), and away from this region, mode IV is only very weakly unstable. In figure 8(b), we plot $\mathrm {Re}(\sigma _n)$ as a function of $k$ for modes I and IV at ${\textit {Fr}} = 2.6$ (green) and $3.2$ (orange). Both curves for mode I are close to the corresponding unstable growth rates in the limiting case $\varepsilon \to 0$ with purely Chézy drag, which asymptote to exactly $0.3$ and $0.6$ respectively, by (3.9).

Figure 8. Effect of the bed erodibility $\varepsilon$ on the growth rates of the unstable modes, showing the suppression of mode I (associated with roll-wave instability) as $\varepsilon$ increases. Panels (a,c,e) show maximum asymptotic growth rates ($\max _i \lambda _{0,i}$) as a function of ${\textit {Fr}}$ for dilute states (solid blue) and individual $\lambda _{0,i}({\textit {Fr}})$ curves (grey dashed), as labelled. Maximum growth rates over all wavenumbers are plotted with blue circles (as in figure 7, which shows the $\varepsilon = 0.01$ case). Panels (b,d,f) show the growth rates $\sigma _n(k)$ of individual modes in dashed lines, as labelled, for ${\textit {Fr}} = 2.6$ (green) and ${\textit {Fr}} = 3.2$ (orange). The maxima of these curves are overlaid as solid lines. The erodibility values for each row are (a,b) $\varepsilon = 3\times 10^{-4}$, (c,d) $\varepsilon = 3\times 10^{-3}$ and (e,f ) $\varepsilon = 6\times 10^{-3}$.

We increase $\varepsilon$ by an order of magnitude in figure 8(c–f). When $\varepsilon = 3\times 10^{-3}$, growth in the $1 < {\textit {Fr}} \lesssim 2.3$ range is now dominated by the morphodynamic mode and $\lambda _{0,1}$ has substantially decreased, remaining negative until ${\textit {Fr}} \approx 2.56$. However, $\mathrm {Re}(\sigma _1)$ no longer attains its maximum in the high-$k$ limit. Indeed, when ${\textit {Fr}} \gtrsim 2.3$, the most significant growth comes from mode I at finite $k$, as evidenced by figure 8(c). By comparing the $\mathrm {Re}(\sigma _1)$ curves in figure 8(d) with their counterparts in figure 8(b) we see that growth of mode I is suppressed in the short-wavelength regime as $\varepsilon$ increases. On proceeding to $\varepsilon = 6\times 10^{-3}$, these trends continue. The plots in figure 8(e,f) show that growth of the morphodynamic mode dominates over the range $1 < {\textit {Fr}} \lesssim 2.8$, before being surpassed by mode I growth at low $k$. Interestingly, mode I is now only unstable across a limited band of wavelengths. Numerical inspection at indicates that at finite $k$, away from their asymptotic limiting expressions in (3.19), each eigenmode contains components of all the flow variables $h$, $u$, $\psi$ and $b$. Therefore, at high enough $\varepsilon$, there are no ‘purely hydraulic’ instabilities, including classical roll waves. However, this does not prohibit the existence of roll waves that are intrinsically coupled with the bed and concentration dynamics.

Due to the complexity of the analytic expression for $\lambda _{0,1}$, given in (3.22a) and (3.23) it is difficult to appreciate directly why mode I ultimately becomes suppressed when morphodynamic effects are significant. Instead, we can look at a simplified case, where the drag function is purely fluid like, $\tilde {\tau } = \tilde {\tau }_f = C_d \tilde {\rho } \tilde {u}^2$, and the bottom friction velocity $\tilde {u}_b$ is neglected. Then, $\tau _{h_0} = \varGamma _{h_0} = \upsilon _0 = 0$, $\tau _{u_0} = 2$ and (3.22a) simplifies to

(4.10)\begin{equation} \lambda_{0,1} = \frac{{\textit{Fr}} - 2}{2} + \frac{\varGamma_{u_0}}{4{\textit{Fr}}({\textit{Fr}} + 1)} \left[ (1 - \rho_b)({\textit{Fr}} - 1) - 2\rho_b{\textit{Fr}} \right]. \end{equation}

The first term in this expression stems from the hydraulic limit and yields the correct stability threshold in that case (${\textit {Fr}} = 2$). Provided that both ${\textit {Fr}} > 1$ and $\rho _b > 1$ (i.e., the bed density exceeds the bulk density), the second term is strictly negative. Moreover, on referring back to the mass transfer closures (4.3) and (4.4), we see that it is proportional to $\varepsilon$. Therefore, greater erodibility decreases $\lambda _{0,1}$ and correspondingly increases the stability threshold for mode I.

In the next two subsections, we investigate the effect of eddy viscosity and bed load on linear growth rates in the morphodynamically dominated regime. We do not independently investigate varying $\varepsilon$ in these cases. However, our numerical observations indicate that, away from ${\textit {Fr}} = 1$, its essential effect is preserved. Namely, that as $\varepsilon$ increases from a negligible value, the $O(\varepsilon )$ bed mode (IV) growth rates become larger and there is a corresponding suppression of mode I.

4.6. Effect of eddy viscosity

Figure 9 demonstrates the effect of the eddy viscosity term studied in § 3.3. The size of the parameter $\nu$ sets the scale over which diffusive effects are important. If chosen sufficiently small, the term only influences the high-$k$ regime. In figure 9(a), we plot maximum normal mode growth rates for dilute solutions with ${\textit {Fr}} = 0.5, 1.5$, $\varepsilon = 0.01$, $\nu = 10^{-4}$ (solid lines) and compare these with the corresponding rates for the unregularised system, $\nu = 0$ (dashed lines). At low $k$ (including $0 \leq k < 10$, not shown), the growth rates are essentially unchanged by the presence of the dissipative term. Then, at higher $k$, both rates for the regularised problem converge to exactly zero, where they remain in the high-$k$ limit. These may be compared with the analytical formulae in (3.40a–c). Computing them for our particular parameters confirms that both cases are dominated by the asymptote $\sigma = 0 + O(k^{-1})$, corresponding to perturbations of the bed (3.40b). For the ${\textit {Fr}} = 0.5$ case, the growth rate increases when it approaches this limit at high $k$. However, note that since it approaches zero from below, this does not affect flow stability. Therefore, we conclude that, away from the singularity at ${\textit {Fr}} = 1$, the addition of the diffusive term succeeds in damping out short-wavelength disturbances without affecting the system outside the asymptotic regime. The behaviour at ${\textit {Fr}} = 1$ is shown in figure 9(b) and confirms the successful regularisation of the growth rate singularity. While the unregularised rate diverges like ${\sim }k^{1/2}$ (as shown in § 3.2.2), when $\nu = 10^{-4}$ it decays to zero within $0 \leq k \lesssim 10^3$. It reaches a maximum growth rate of approximately $2$, around 3–6 times the magnitude of unregularised growth rates either side of the singularity at ${\textit {Fr}} = 0.75$ and $1.25$, plotted in figure 6(a).

Figure 9. Effect of the eddy viscosity regularisation term on maximum growth rates. Curves for the dilute solution branch are shown with $\nu = 0$ (dashed) and $10^{-4}$ (solid). The Froude numbers are: (a) ${\textit {Fr}} = 0.5$ (orange) and $1.5$ (red); (b) ${\textit {Fr}} = 1$ (purple).

Figure 10 shows the effect of regularisation on concentrated states. In these plots we include additional detail, plotting the growth rates for all four modes. Dashed lines show curves with $\nu = 0$; dotted lines show the corresponding rates with $\nu = 10^{-4}$. We label these I–IV by numerically computing the asymptotic rates (3.22a–d) and (3.27) for the unregularised system, which match the corresponding regularised rates at lower $k$. The case of ${\textit {Fr}} = 0.5$, away from the singularity is given in figure 10(a). Both hydrodynamic modes (I and II) are severely damped at high $k$. On checking these curves against the asymptotic rates in (3.40a–c), we find that mode I corresponds to (3.40a), scaling as $\sim -k^2$ asymptotically, while mode II eventually converges to $\lambda _- \approx -4 \times 10^4 = -1/\nu {\textit {Fr}}^2$ in this case (outside the range of the figure axes). The mode IV rate increases as $k$ increases, eventually asymptoting to $0$. Finally, as argued in § 3.3, mode III, which asymptotes to $\lambda _+$ ($\approx 0.015$), is not much affected by the eddy viscosity term, provided that $\nu$ is small relative to the other terms in (3.38). The situation at the ${\textit {Fr}} = 1$ singularity for concentrated states is shown in figure 10(b). When $\nu = 0$, the mode II and IV growth rates can be seen diverging to $-\infty$ and $+\infty$ respectively as $k\to \infty$. As established in § 3.2.2, these scale like ${\sim }\pm k^{1/2}$, asymptotically – while modes I and III converge to $\lambda _{0,1}$ ($\approx 0.04$) and $\lambda _{0,3}$ ($\approx 0.07$) respectively. With the addition of regularisation, the modes qualitatively mirror the ${\textit {Fr}} = 0.5$ case: the hydraulic modes are strongly damped, while mode IV no longer diverges and decays to zero, being eventually surpassed by mode III which is essentially unaffected by the eddy viscosity.

Figure 10. Individual growth rates as a function of $k$ for concentrated steady states. The curves are labelled I–IV according to their corresponding linear stability modes in the $\nu = 0$ analysis of § 3.2.2. Unregularised rates are plotted with dashed grey lines; regularised rates ($\nu = 10^{-4}$) with dotted lines. The maximum over all curves in the $\nu = 10^{-4}$ case is overlaid as a solid line. The Froude numbers are (a) 0.5 and ($\textit {b}$) 1.

We have briefly experimented with introducing an additional diffusive term to the solid mass transport equation. This takes the form $({\partial }/{\partial x})(\kappa h ({\partial \psi }/{\partial x}))$ and is added to the right-hand side of (3.3b), thereby contributing an additional non-zero term ${\mathsf{D}}_{23} = \kappa$ to the matrix $\boldsymbol{\mathsf{D}}$ in the linear stability eigenproblem (3.28). The free parameter $\kappa = \tilde {\kappa }/ (\tilde {u}_0 \tilde {l}_0)$, where $\tilde {\kappa }$ is a (constant) characteristic sediment diffusivity sometimes included in models (e.g. Balmforth & Vakil Reference Balmforth and Vakil2012; Bohorquez & Ancey Reference Bohorquez and Ancey2015). Its inclusion is consistent with the basic physics of sediment transport. However, the effect on the growth rates is modest, for the values explored ($10^{-3} < \kappa < 10^{-5}$), effectively serving only to damp mode III at high wavenumbers. Consequently, figure 9 is unaffected, while the maximum unstable growth in figure 10 ultimately rolls off at high $k$, with the roll-off being more severe for larger values of $\kappa$. Therefore, in the context of simple shallow-flow modelling, sediment diffusion, combined with momentum diffusion via eddy viscosity, could serve as an unobtrusive way to prevent instabilities from developing over unphysically short length scales.

We now examine the effect of varying $\nu$. Figure 11 shows the reduction of growth rate at the ${\textit {Fr}} = 1$ singularity as $\nu$ is increased over three orders of magnitude for both the (a) dilute and (b) concentrated solution branches. At $\nu = 10^{-4}$, the maximum growth is curtailed only in a relatively small neighbourhood of ${\textit {Fr}} = 1$. Consequently, growth still peaks sharply in this region and elsewhere follows the asymptotic rate for the $\nu = 0$ case (solid curves). Increasing $\nu$ smooths over the signature of the singularity: its (diminished) influence is clear at $\nu = 10^{-3}$, but by $\nu = 10^{-2}$ there is only a residual trace of it. For the larger values of $\nu$, growth largely fails to reach the asymptotic rates at all. We note also that increasing $\nu$ increases both the Froude number at which maximum growth occurs and, in the case of dilute states, the onset of instability.

Figure 11. Effect of eddy viscosity size on the severity of morphodynamic instability. Solid curves show the maximum asymptotic growth rates for the unregularised system (as shown previously in figure 7), near the ${\textit {Fr}} = 1$ singularity, for (a) dilute and (b) concentrated steady states. Dotted curves show the maximum growth rates over all wavenumbers and all four modes, for various $\nu$, as labelled.

A constant eddy viscosity is a crude parametrisation of the effects of turbulent dissipation. Consequently, studies of non-erosive shallow layers have typically treated selection of $\nu$ (when included) as a means to an end – either to smooth over hydraulic jumps in flow profiles, or to constrain instabilities to within a bounded spectrum (Needham & Merkin Reference Needham and Merkin1984; Hwang & Chang Reference Hwang and Chang1987; Balmforth & Mandre Reference Balmforth and Mandre2004). Provided $\nu$ is sufficiently small, this is reasonable since roll-wave onset and development tends to be largely insensitive to the exact magnitude of the dissipation term (Chang, Demekhin & Kalaidin Reference Chang, Demekhin and Kalaidin2000). However, figure 11 suggests this is not the case when the evolution of the bed is accounted for. Here, the size of $\nu$ necessarily dictates the severity of the morphodynamic instability. A plausible range for $\nu$ can be ascertained as follows. Suppose that dissipation acts over length and time scales set by shearing within the turbulent boundary layer. Then we may anticipate $\tilde {\nu } \sim \tilde {u}_* \tilde {h}$, where $\tilde {u}_*$ is the basal friction velocity, equivalent to the closure used for $\tilde {u}_b$, given in (4.6). For a dilute steady state flow, $\tilde {u}_* \approx \tilde {u}_0 \sqrt {C_d}$. In our dimensionless units, $\nu = \tilde {\nu } g \sin \phi / \tilde {u}_0^3 = \tilde {\nu } \tan \phi / (\tilde {h}_0 \tilde {u}_0 {\textit {Fr}}^2)$ and (since $\tan \phi \approx C_d {\textit {Fr}}^2$ in the dilute regime) therefore $\nu \sim C_d^{3/2}$. With our chosen drag coefficient, this yields $\nu \sim 10^{-3}$. After accounting for potential deviations from this rough order-of-magnitude estimate and especially the ways in which suspended sediment may suppress turbulent fluctuations, it is difficult to rule out any of the scenarios in figure 11 with confidence. However, it is at least reasonable to conclude that diffusive effects are neither negligible, nor likely to completely quash unstable growth near unit Froude number in morphodynamic shallow-flow models.

4.7. Effect of bed load

In § 3.4, we showed that the inclusion of a bed load flux $Q$ also prevents ill posedness, provided that the constraint (3.45) is satisfied. It does so by altering the system characteristics, thereby avoiding a resonance between modes II and IV at ${\textit {Fr}} = 1$ that otherwise leads to unbounded linear growth. In fluvial settings, the physical basis for such terms is well established (see e.g. Gomez Reference Gomez1991) and similar analyses to ours have noted the role it plays in ensuring strict hyperbolicity of models (Cordier et al. Reference Cordier, Le and Morales de Luna2011; Stecca et al. Reference Stecca, Siviglia and Blom2014; Chavarrías et al. Reference Chavarrías, Stecca and Blom2018, Reference Chavarrías, Schielen, Ottevanger and Blom2019). In more severely concentrated situations, such as debris and purely granular flows, the case for a bed load is less clear since it may not be possible to distinguish between grains crawling along the bed surface and grains in suspension. Therefore, bed load fluxes are not typically included in numerical models of these flows. However, since bed load dictates the characteristic wave speed of the bed surface, it may nonetheless have a role to play in these settings that is not currently appreciated.

We shall consider the effect of bed load separately from eddy viscosity by resetting $\nu =0$ and employing the flux closure in (4.5). Although the effect on concentrated steady states was briefly investigated, the conclusions did not differ substantially from the relatively dilute case, whose results we report below. In figure 12(a), we plot the maximum growth rates over all modes and wavenumbers as a function of ${\textit {Fr}}$ for various values of the bed load flux strength $\gamma$, increasing from zero (solid blue) to $\gamma = 1$, $5$ and $10$ (dotted). The plot may be compared with figure 11(a) and shares the same essential features. This should not be surprising: relatively small values of $\gamma$ do not displace the system characteristics $\lambda _1$ far from their resonant values, while larger values push them far apart, leaving little trace of the ${\textit {Fr}} = 1$ singularity. We observe this directly in figure 12(b), which plots $\lambda _1({\textit {Fr}})$ for each $\gamma$, as labelled (solid and dotted lines), using the same horizontal axes. These curves are obtained by numerically solving (3.17). The solid lines are the characteristics $\lambda _1 = 1 - {\textit {Fr}}^{-1}$ (mode II) and $\lambda _1 = 0$ (mode IV) of the system without bed load, as derived in (3.18). As $\gamma$ increases from zero, the curves diverge from the singular intersection point of these two characteristics. Typical values of $\gamma$ from the fluvial literature are $O(10)$ (see e.g. Cordier et al. Reference Cordier, Le and Morales de Luna2011), around the upper end of the range considered here. In this regime, growth rates are relatively modest and the onset of instability notably increases with $\gamma$, reaching ${\textit {Fr}} \approx 1.15$, when $\gamma = 10$. Where solutions are unstable, we also plot, in figure 12(b), points along $c_{max}({\textit {Fr}})$, which we define here to be the wave speed $c = -\mathrm {Im}(\sigma ) / k$ of the dominant unstable mode at its most unstable wavenumber, i.e., the mode (for each $\gamma$) whose corresponding growth rates are in figure 12(a). Each set of points lies exactly on its corresponding characteristic curve. This is because, in this case, the maximum growth at each ${\textit {Fr}}$ occurs in the high-$k$ regime and $c\to \lambda _1$ as $k\to \infty$. Furthermore, we note that since $c_{max} < 0$ for $\gamma > 0$, dominant perturbations travel upslope (and do so more rapidly if the bed load strength is larger). When $\gamma = 10$, the picture in figure 12(a,b) shares traits with some fluvial systems – there is neither ill posedness, nor severely accentuated growth near unit Froude number, and the morphodynamics drives slowly upstream-migrating bedforms (see e.g. Colombini & Stocchino Reference Colombini and Stocchino2008; Seminara Reference Seminara2010). It is tempting to think of the morphodynamic processes in this regime as being essentially ‘bed load dominant’. However, pure bed load formulations do not feature instabilities near ${\textit {Fr}} = 1$ (Lanzoni et al. Reference Lanzoni, Siviglia, Frascati and Seminara2006). Indeed, we have checked that decreasing $\varepsilon$ to $10^{-4}$ (and retaining $\gamma = 10$), removes the morphodynamic instability in our model, which then remains stable until the threshold for roll waves, near ${\textit {Fr}} = 2$. Therefore, the bulk mass transfer term (i.e., suspended load) plays a role in sustaining the instabilities of figure 12(a).

Figure 12. Effect of a bed load flux term on linear growth rates and the corresponding instabilities. Panel (a) shows curves of maximum growth rate over all modes and wavenumbers as a function of ${\textit {Fr}}$ (dotted), for $\varepsilon = 0.01$, $d = 0.01$ and various $\gamma = 1$, $5$ and $10$ as labelled. Also shown for reference, is the analytical curve of maximum growth rate at high wavenumber (solid blue), for $\gamma = 0$. In (b) we plot the system characteristics $\lambda _1({\textit {Fr}})$ with a solid blue line for $\gamma = 0$ and dotted lines for the $\gamma > 0$ cases. In the unstable regions, we also plot the corresponding wave speed $c_{max}({\textit {Fr}})$ of the dominant unstable mode at its most unstable wavenumber (filled circles). Panel (c) shows maximum growth rate curves, as in (a), for $\varepsilon = 6\times 10^{-3}$, $\gamma = 8$ and various $d = 10^{-3}$ (green dotted), $5\times 10^{-3}$ (orange dotted), $10^{-2}$ (black dotted). The reference curve $\max _i \lambda _{0,i}$ (solid blue) is plotted for $\gamma = 0$ and $d = 10^{-3}$. In (d), the characteristics (solid and dotted lines for $\gamma = 0$ and $\gamma = 8$ respectively) and $c_{max}({\textit {Fr}})$ (filled circles) are plotted using the same colour scheme as (c), with labels I–IV to indicate the corresponding eigenmodes.

Furthermore, even if $\gamma$ is large enough that the characteristics are well separated, it is possible to see the influence of the ${\textit {Fr}} = 1$ singularity if we move to a regime where suspended load is enhanced. In figure 12(c), we fix $\gamma = 8$, $\varepsilon = 6 \times 10^{-3}$ and plot maximum growth rate curves, in the vein of figure 12(a), for different $d = 0.01$, $5\times 10^{-3}$ and $10^{-3}$. Also shown for reference is the singular high-$k$ growth curve (solid blue), for $d = 10^{-3}$ and no bed load. Note that these curves may only be plotted for ${\textit {Fr}}$ where steady flows exist and this range shrinks as $d$ decreases. (This is because smaller particles are more easily eroded, meaning that the unsteady scenario of figure 4, where erosion always exceeds deposition, occurs at lower ${\textit {Fr}}$.) The trend of figure 12(c) shows that decreasing particle size leads to more severe growth near ${\textit {Fr}} = 1$, with the $d = 10^{-3}$, $\gamma = 8$ case (dotted green) inheriting a severe instability in this region. The corresponding characteristics, plotted in figure 12(d) below, are not greatly affected by changes in $d$. Therefore, enhanced growth around unit Froude number cannot be due to near intersection of characteristics in this case. Without simple analytical expressions for the growth rates when $Q \neq 0$, it is difficult to pin down exactly why small $d$ has this effect. However, we note that since our closures for erosion and bed load, given in (4.3) and (4.5) respectively, have similar dependencies on the excess shear stress, it is straightforward to see for any given ${\textit {Fr}}$, that $Q / E \propto d\gamma / \varepsilon$. Therefore, smaller $d$ implies that the magnitude of bed load is diminished, relative to the suspended load dynamics.

Finally, just as we saw in figure 12(b), we note that for the $d = 10^{-3}$ and $d = 5\times 10^{-3}$ case, $c_{max}$ (filled circles) lies exactly on the characteristic curves corresponding to slowly upstream-propagating disturbances. However, for the $d = 0.01$ branch, this mode is only dominant for Froude numbers up to approximately $2.4$, where the rapidly downstream-propagating perturbations of mode I become the most unstable. This is marked in figure 12(c) by a steepening of the maximum growth rate curve and may be compared to the curves in figure 8(a–d), which depict an analogous transition between modes IV and I in the $Q=0$ setting. The latter mode is related to hydraulic roll-wave instability, as discussed in § 4.5. Note that in this region (${\textit {Fr}} \gtrsim 2.4$), $c_{max}$ does not precisely follow the mode I characteristic. This is because, as in figure 8, maximum growth occurs at finite $k$, rather than in the asymptotic regime where disturbance wave speeds are given exactly by the characteristics.

4.8. Summary

Finally, we return to the unregularised suspended load model and explore the $(d,{\textit {Fr}})$-parameter space more broadly, by computing the steady states that can exist and their linear stability, when the solid diameter $d$ is varied over three decades. We assume here that the non-dimensional settling speed is constant and unity (as in table 1) throughout, even though its dimensional counterpart $\tilde {w}_s$ varies considerably with the particle size. This is reasonable for sufficiently large particles, since $\tilde {w}_s \approx \tilde {u}_p$, when the particle Reynolds number is high. (This is straightforward to confirm from typical empirical formulae for $\tilde {w}_s$, see e.g. Cheng Reference Cheng1997; Soulsby Reference Soulsby1997). Figure 13 summarises the existence of (a) dilute and (b) concentrated steady flows across parameter space. Where steady solutions exist, we contour them according to $\psi _0$. Elsewhere, we leave the region blank. Overlaid are contours showing lines of constant slope angle. Care must be taken when interpreting this plot, since its values depend on the choice of parameters. However, its qualitative characteristics are robust. The most striking observation is the separation of the two states, which are either highly dilute or highly granular. This is clear from considering the picture in figure 4. Solutions exist predominantly for higher $d$, where erosion rates are typically smaller and may therefore balance deposition at higher ${\textit {Fr}}$. As ${\textit {Fr}}$ increases, keeping $d$ fixed, eventually $E_p > D_p$ and states can no longer be steady. Bounding the region of existence from above is the unidimensional family of states of type (b) in figure 4, where the dilute and concentrated branches coalesce. Consequently, the dilute and concentrated states respectively possess greater and lesser solid fractions close to this boundary than they do in the bulk of parameter space. On the dilute contour map, figure 13(a), we also plot the neutral stability curve (solid black), below which states are stable to both hydraulic and morphodynamic modes of disturbance. It is determined (for our chosen closures) by the asymptotic growth rate $\lambda _{0,2}$ (which diverges to $+\infty$ at ${\textit {Fr}} = 1$) crossing zero. Therefore, we compute it by numerically solving $\,f_-({\textit {Fr}}) = 0$, where $f_-$ was given in (3.23). As indicated by figure 7(a), the upper stability limit lies just below the line ${\textit {Fr}} = 1$ for larger $d$ values. At smaller $d$, the neutral line dips as it approaches the existence boundary where solutions are more concentrated. For the parameters used herein, the region of stability includes only extremely shallow grades (typically less than $1^\circ$). However, this is not unexpected, since fluid drag predominates and purely Chézy layers turn supercritical (in this case) when $\phi = \arctan (C_d) \approx 0.6^\circ$. Beneath the region where stable dilute erosive flows exist is a region where flows are not sufficiently energetic to entrain material, $\theta < \theta _c$, indicated by a dotted line. Here, only steady flows with zero solid fraction exist. In the concentrated case, this region is only very narrow. Note that for any $d$, the limiting solutions as ${\textit {Fr}} \to 0$ on this branch are static granular layers resting at the neutral slope angle $\arctan (\mu _1) \approx 5.7^\circ$ [see (4.2) for the definition of $\mu _1$]. Such states typically have high Shields numbers in excess of the constant part $\theta _c^*$ of the critical value and consequently any increase in ${\textit {Fr}}$ from zero leads to entrainment. The remainder of parameter space in the concentrated case features steady flows at a range of more severe slope angles, all of which are unstable.

Figure 13. Existence and stability of uniform steady layers in the morphodynamic case. Plotted on two axes are filled contours of scaled solid fraction $\psi$ (shaded regions) for the (a) dilute and (b) concentrated solution branches, over a representative range of $d$ and ${\textit {Fr}}$ values. The region where steady erosive flows exist is outlined by dotted black lines. Where no steady flows exist, the plot is left blank. The dilute states possesses a region of stability; the corresponding neutral curve is shown as a solid black line in (a). Overlaid dashed contours indicate lines of constant slope angle.