Sorting and bed waves in unidirectional shallow-water flows

Abstract

The stability of a uniform flow above an erodible bed composed of a bimodal mixture of sediments is investigated by means of linear analysis. Results show that, for any given set of the flow and sediment parameters, two distinct modes of instability arise, each one characterized by its own wave speed, growth rate and longitudinal wavelength, each one involving spatial variations of both grain size density and bed elevation. Although at a linear level no information on the amplitude of the perturbations is gathered, the analysis of the eigenvectors associated with the two modes of instability allows for an easy classification in terms of the relative amplitudes of the perturbations of bed elevation and size density. One eigenvalue is shown to be associated with the modifications of bed forms induced by the presence of the heterogeneous mixture, such as the local accumulation of finer and coarser material along the unit wavelength, the other with the formation of the low-amplitude sorting waves known as bedload sheets. In the present unidirectional shallow-water framework, only the sorting wave is found to be unstable, since dunes and antidunes, the relevant bed forms for this case, require a more refined rotational flow model in order to become unstable. On the other hand, the simple flow model adopted allows for the formulation of an algebraic eigenvalue problem that can be solved analytically, allowing for a deep insight into the mechanisms that drive both instabilities.

1 Introduction

Sediments in nature can seldom be considered as homogeneous. Despite this, among the many simplifying assumptions usually adopted in the modelling of sediment transport, the one of well-sorted material is probably the most common. Following this approach, the dynamics of the sediment mixture is modelled as if it were composed of a single grain size, typically chosen as the median diameter of the mixture, i.e. the average particle diameter of the sample by mass. If the sediment is moving as bedload, a single sediment mass conservation equation (known as the Exner equation) is then usually introduced to represent the bed dynamics. Moreover, sediment discharge is assumed to equal the transport capacity of the flow. Once coupled with the set of equations of a suitable flow model, the Exner equation provides the simplest possible tool to study morphodynamic problems.

Figure 1. Side and plane view of bedload sheets as reported by Whiting et al. (1988). Flow is from left to right.

However, simplification always comes at a price, which in this case is related to the fact that sorting effects are completely neglected in the analysis. This can be acceptable when, as in the study of bed forms, the focus is on the instability of the bed interface as a whole and sorting is passively driven by changes in bed elevation, simply producing an accumulation of finer (coarser) material on crests (troughs) or vice versa. On the other hand, this becomes inappropriate when the formation of the bed form itself is inherently associated with, or dominated by, the effect of sorting, the heterogeneity of the sediment being the crucial mechanism driving the instability (Seminara 1995; Livesey et al. 1998). This is the case of the sorting waves named ‘bedload sheets’ after Whiting et al. (1988), which appear as rhythmic alternations of finer and coarser bands of the bed material aligned across the flow (see figure 1), characterized by downstream migration and negligible amplitude (Venditti, Nelson & Dietrich 2008). Bedload sheets have been observed both in the laboratory (Iseya & Ikeda 1987; Kuhnle & Southard 1988; Bennett & Bridge 1995; Recking et al. 2009) and in the field (Whiting et al. 1988; Cudden & Hoey 2003) sometimes superimposed on the stoss side of dunes and bars (Whiting et al. 1988; Bennett & Bridge 1995; Livesey et al. 1998; Rice et al. 2009). Moreover, their propagation has been associated with the rhythmic fluctuations of the bedload transport rate that are frequently observed in gravel bed rivers (Reid, Frostick & Layman 1985; Iseya & Ikeda 1987; Kuhnle & Southard 1988).

Several authors (Livesey et al. 1998; Cudden & Hoey 2003) have highlighted how difficult it is, from an experimental point of view, to distinguish among the different spatial and temporal scales associated with bed and sorting waves. This led Livesey et al. (1998) to conclude that ‘the interactions between bed topography, flow structure and surface sorting patterns remains poorly understood despite their central role in controlling mixed-size entrainment dynamics and fractional bedload transport rates’. In particular, the distinction between dunes with sorting and bedload sheets is difficult to make (Carling 1999), so that often the latter have been considered as some sort of dune ‘precursors’ (e.g. Bridge 1993), adding confusion to the interpretation. In fact, observations have shown that both the wavelength and the wave speed of dunes (with sorting) and bedload sheets can be remarkably different from one another, the former being typically longer and slower than the latter, as shown by Kuhnle et al. (2006), who explicitly state that ‘It is not entirely clear how bed load sheets found in sand and gravel mixtures relate to the dunes in gravel found in other studies’. The most distinctive character among the two patterns remains indeed the height of the bed form, which is almost negligible for bedload sheets as compared to dunes.

From a theoretical point of view, linear stability analyses are known to be the ideal tool to shed light on the mechanisms responsible for the formation of a wide variety of patterns, so it is not surprising that in the nineties several attempts were made to include sorting into existing studies on bed forms, initially developed for a well-sorted sediment (see Colombini & Parker (1995), Seminara, Colombini & Parker (1996) and Lanzoni & Tubino (1999), among others). These analyses were mainly focused on the formulation of coupled sediment transport and flow models valid for heterogeneous mixtures, thus including the effect of hiding and of the local surface roughness. The appearance of sorting waves was in fact acknowledged, but the distinct nature of topographic- and roughness-driven instabilities did not receive the attention it deserved, mainly because the focus was more on the effect of sorting on bed forms than on the formation of sorting waves. The only notable exception is the work of Seminara et al. (1996), where the formation of bedload sheets was investigated in detail. However, in order for the bedload-sheet mode to emerge, in Seminara et al. (1996) the dune mode was artificially ruled out from the analysis, by assuming that the perturbation of bed elevation scales with the (small) standard deviation of the grain size distribution.

More recently, in the scientific community dealing with coastal morphodynamics, the problem of bed forms in heterogeneous sediment has received considerable attention. In particular, Murray & Thieler (2004) described the existence in the inner continental shelf of ‘self-organized sorted bed forms’, which can be considered as the coastal counterpart of bedload sheets. Sorted bed forms, although the name chosen by Murray & Thieler (2004) can be misleading in this regard, are in fact sorting patterns generated by roughness variations, coupled to small undulations in the bed level (Van Oyen, de Swart & Blondeaux 2010), the occurrence of which has been explained by both the numerical solution of Coco et al. (2007) and the linear stability analysis of Van Oyen, de Swart & Blondeaux (2011). In particular, in the latter work the presence of two distinct modes of instability was highlighted, a ‘topography-driven’ mode primarily associated with changes in bed elevation and a ‘roughness-driven’ mode, characterized by shorter wavelength and by smaller amplitudes than the former, associated with the appearance of the ‘sorted bed forms’. With the present contribution, which can be considered as a revisitation of Seminara et al. (1996) with a simpler flow model, those concepts are brought back to the riverine applications where they were originally conceived. By doing so, we hope to gain a deeper insight into the mutual interactions between sorting and bed forms.

In the following, the terms ‘bed wave’ and ‘sorting wave’ have been adopted to indicate the topography-driven and the roughness-driven instabilities, respectively. In fact, these terms have been borrowed from Stecca, Siviglia & Blom (2014), where a linearized analysis of the dynamics of small waves in the case of heterogeneous sediment was carried out. Note that each of these modes of instability displays a topographic as well as a sorting pattern. It is only the relative importance of the amplitude of the perturbations of bed elevation and size density that, eventually, draws the distinction between bed waves (dunes) and sorting waves (bedload sheets).

The present analysis might also be of relevance for the investigation of the well posedness of the Saint-Venant–Hirano model, which has been recently questioned in several studies in terms of either loss of hyperbolicity of the system (Ribberink 1987; Sieben 1997; Stecca et al. 2014; Chavarrías, Stecca & Blom 2018) or unboundedness of the growth rate in the short-wave range (Chavarrías et al. 2019). Although the focus is herein on the stability of the governing differential system and therein on its mathematical properties, the concepts of stability, unboundedness and well posedness (and their own opposites) are in fact strongly tied to one another and all the above analyses ultimately rely on a linearization in terms of small-amplitude hydro-morphodynamic waves. Indeed, if the growth rate associated with an eigenvalue has a finite upper bound, independent of the wavenumber, and so is well posed in the sense of Chavarrías et al. (2019), the eigensystem is defined ‘regular’ (Birkhoff 1954). In turn, regularity implies hyperbolicity and so if the solution is regular in the short-wave range it is also mathematically well posed in the sense of Stecca et al. (2014).

Note that it is the unboundedness of the growth rate that makes the system irregular and ill posed, not the mere instability (i.e. a positive growth rate) in the short-wave range, though the numerical solution of a system where the maximum growth rate is found in the limit of infinitely short waves can become problematic. Examples of regular short-wave instability are indeed present in the literature (Hooper & Boyd 1983). It must be stressed, however, that the absence of a cutoff in the short-wave range makes the solution much less appealing from a physical point of view, because no wavelength of maximum amplification is selected by means of a normal mode analysis. As Joseph & Saut (1990) pointed out, that ‘there is always a cutoff’ is in fact an axiom in physics which can be considered as fully equivalent to the medical aphorism ‘the bleeding always stops’. Moreover, the absence of a cutoff in the short-wave range is quite often the symptom of the lack of an essential, stabilizing ingredient in the analysis, as is the case of the interfacial tension in Kelvin–Helmholtz and Rayleigh–Taylor instabilities, which transform an ill-posed (and irregular) problem into a well-posed one (Truzzolillo & Cipelletti 2017).

2 Formulation of the problem

The starting point of our analysis is the one-dimensional form of the governing equations of morphodynamics in a straight channel with an erodible bed. The triplet composed by the fluid mass density $\unicode[STIX]{x1D70C}$ and by the uniform, depth-averaged flow velocity $U^{\ast }$ and depth $D^{\ast }$ is used for non-dimensionalization. Here and in the following, an asterisk denotes dimensional variables.

Figure 2. Sketch of a longitudinal section of the flow configuration.

It is convenient to write the equations in the Cartesian coordinate system sketched in figure 2, sloping with slope $S$ . Three different interfaces are displayed in the picture. The curve $z=B(x,t)$ is the actual bed level and the curve $z=B(x,t)+R(x,t)$ represents the reference level, $R(x,t)$ being the roughness height, i.e. the distance above the bed at which the logarithmic vertical profile of velocity conventionally vanishes. Moreover, we stipulate that the reference level sets the lower boundary of the flow, so that the free surface is represented by the curve $z=B(x,t)+R(x,t)+D(x,t)$ , where $D$ is the non-dimensional local flow depth.

2.1 The flow model

Under the shallow-water approximation, whereby a hydrostatic pressure distribution along the vertical is assumed, the non-dimensional forms of the Saint-Venant and continuity equations read

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle U_{,t}+UU_{,x}-\frac{S}{Fr^{2}}+\frac{(B+R+D)_{,x}}{Fr^{2}}+\frac{T_{t}}{D}-\frac{(T_{n}D)_{,x}}{D}=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle D_{,t}+(UD)_{,x}=0, & \displaystyle\end{eqnarray}$$

where $U$ is the local value of the depth-averaged velocity, $Fr=U^{\ast }/\sqrt{gD^{\ast }}$ is the Froude number of the undisturbed uniform flow, $g$ is the gravitational acceleration, $T_{t}$ is the bed shear stress, $T_{n}$ is the depth-averaged value of the streamwise normal Reynolds stress and the comma notation is used to indicate partial derivatives. The closure for the stresses is found by assuming a self-similar vertical profile of velocity and mixing length of the kind

(2.3a,b ) $$\begin{eqnarray}u=\frac{u_{f}}{\unicode[STIX]{x1D705}}\ln \left(\frac{R+\unicode[STIX]{x1D701}D}{R}\right),\quad l=\unicode[STIX]{x1D705}(R+\unicode[STIX]{x1D701}D)(1-\unicode[STIX]{x1D701})^{1/2},\end{eqnarray}$$

where $u_{f}$ is the friction velocity, $\unicode[STIX]{x1D705}$ is the von Kármán constant, taken as 0.4 and the following transformation is employed, which maps the domain of figure 2 in a rectangular domain:

(2.4a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D701}=\frac{z-B-R}{D},\quad \unicode[STIX]{x1D709}=x.\end{eqnarray}$$

Hence

(2.5a,b ) $$\begin{eqnarray}T_{t}=u_{f}^{2}=\unicode[STIX]{x1D708}_{t}u_{,z}|_{\unicode[STIX]{x1D701}=0},\quad T_{n}=2\int _{0}^{1}\unicode[STIX]{x1D708}_{t}u_{,x}\,\text{d}\unicode[STIX]{x1D701},\quad \unicode[STIX]{x1D708}_{t}=l^{2}u_{,z}.\end{eqnarray}$$

Moreover, the friction velocity is proportional to the local value of the depth-averaged velocity through a non-dimensional Chézy coefficient $C$ , which, in turn, is assumed to depend on the average bed roughness $r^{\ast }$ through the Keulegan relationship (ASCE 1963), valid for a uniform flow in the rough regime

(2.6) $$\begin{eqnarray}C=\frac{U}{u_{f}}=\frac{1}{\unicode[STIX]{x1D705}}\ln \left(\frac{11.09D^{\ast }}{r^{\ast }}\right).\end{eqnarray}$$

In the case of homogeneous sediment, the roughness height $R$ is typically set equal to one thirtieth of the non-dimensional bed roughness, consistently with (2.6). In the present heterogeneous case, the same line of reasoning is followed, whereby the local roughness height $R$ is set to one thirtieth of the local bed roughness, which in turn depends on the local bed composition. As a consequence, a local increase of the percentage of the coarser (finer) fraction implies a larger (smaller) roughness height, thus creating the necessary feedback between the sorting process and the flow.

2.2 The sediment transport model

To treat the bedload transport of heterogeneous sediment, a simple three-layer model is used, which implements the concept of the ‘active layer’ of Hirano (1971). The active layer, depicted as a grey area in figure 2, is defined as the layer of sediment close to the bed interface which is available for entrainment into bedload. It is convenient at this point to briefly summarize some of the salient concepts of the formulation. For a more complete analysis, the reader is referred to Colombini & Parker (1995) and Seminara et al. (1996).

To characterize the size distribution of the mixture composing the bed, we introduce the sedimentological $\unicode[STIX]{x1D719}$ -scale as

(2.7) $$\begin{eqnarray}d^{\ast }(\unicode[STIX]{x1D719})=d_{50}^{\ast }2^{-\unicode[STIX]{x1D719}+\unicode[STIX]{x1D719}_{50}},\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{50}$ is the opposite of the binary logarithm of the median of the grain-size distribution $d_{50}^{\ast }$ , expressed in millimetres.

The probability density function for the grain size $\unicode[STIX]{x1D719}$ in the active layer (in short, the size density) is denoted as $F(\unicode[STIX]{x1D719};x,t)$ , where, by definition

(2.8) $$\begin{eqnarray}\int _{-\infty }^{\infty }F(\unicode[STIX]{x1D719})\,\text{d}\unicode[STIX]{x1D719}=1.\end{eqnarray}$$

The first and second moments of $F$ provide the mean $\unicode[STIX]{x1D719}_{m}$ and standard deviation $\unicode[STIX]{x1D70E}_{m}$ (on the $\unicode[STIX]{x1D719}$ -scale) of the grain-size distribution as

(2.9a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{m}=\int _{-\infty }^{\infty }\unicode[STIX]{x1D719}F(\unicode[STIX]{x1D719})\,\text{d}\unicode[STIX]{x1D719},\quad \unicode[STIX]{x1D70E}^{2}=\int _{-\infty }^{\infty }(\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{m})^{2}F(\unicode[STIX]{x1D719})\,\text{d}\unicode[STIX]{x1D719},\end{eqnarray}$$

from which the characteristic diameters $d_{g}^{\ast }$ and $d_{\unicode[STIX]{x1D70E}}^{\ast }$ can be calculated

(2.10a,b ) $$\begin{eqnarray}d_{g}^{\ast }=d_{50}^{\ast }2^{-\unicode[STIX]{x1D719}_{m}+\unicode[STIX]{x1D719}_{50}},\quad d_{\unicode[STIX]{x1D70E}}^{\ast }=d_{50}^{\ast }2^{-\unicode[STIX]{x1D719}_{m}+\unicode[STIX]{x1D719}_{50}+\unicode[STIX]{x1D70E}}.\end{eqnarray}$$

It is just worth noting that in the case of a log-normally distributed sediment

(2.11a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{m}=\unicode[STIX]{x1D719}_{50},\quad d_{g}^{\ast }=d_{50}^{\ast },\quad d_{\unicode[STIX]{x1D70E}}^{\ast }=d_{50}^{\ast }2^{\unicode[STIX]{x1D70E}}=d_{84}^{\ast },\end{eqnarray}$$

where $d_{84}^{\ast }$ is the size such that 84 % of the mass of a sample is finer. This grain size is close to $d_{90}^{\ast }$ , which is often considered as representative of both the roughness of the bed and the thickness of the active layer. However, in the present contribution we want to investigate the mechanisms which drive and are driven by sorting and, to this end, we are particularly interested in well-sorted mixtures which exhibit only a small degree of heterogeneity. Hence, we can safely make use of $d_{g}^{\ast }$ instead of $d_{\unicode[STIX]{x1D70E}}^{\ast }$ in the determination of the roughness and of the active layer thickness. Moreover, following Colombini & Parker (1995) and Seminara et al. (1996), we assume the active layer to be sufficiently thin to allow the neglect of a vertical structure. Any interaction with the substrate is neglected as well, assuming that only a minimal amount of aggradation and degradation takes place. This implies that the size density $F(\unicode[STIX]{x1D719})$ in the active layer coincides with the areal concentration of the sediment in the size range $\unicode[STIX]{x1D719}$ , $\unicode[STIX]{x1D719}+\text{d}\unicode[STIX]{x1D719}$ . The rationale for these hypotheses is that in the framework of a linear analysis the base uniform state, for which the bed neither aggrades nor degrades, is only slightly perturbed. Indeed, it must be pointed out that these assumptions completely rule out aggradational and degradational cases and, with them, the troublesome, possibly elliptic case of degradation into a finer substrate (Stecca et al. 2014).

The bed shear stress and the volumetric grain-specific sediment discharge per unit width are made non-dimensional using $d_{50}^{\ast }$ , as it is customary when treating a sediment mixture

(2.12a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D703}=\frac{T_{t}^{\ast }}{\unicode[STIX]{x1D70C}(s-1)gd_{50}^{\ast }},\quad q_{s}=\frac{q_{s}^{\ast }}{\sqrt{(s-1)gd_{50}^{\ast }}d_{50}^{\ast }}=q_{u}F,\end{eqnarray}$$

where $\unicode[STIX]{x1D703}$ is the Shields stress, $s$ is the relative mass density of the sediment and $q_{u}(\unicode[STIX]{x1D719})$ is the bedload transport density per unit size density in the active layer (in short, the transport density).

Imposing sediment mass conservation, a grain-size specific non-dimensional Exner equation can eventually be obtained that reads

(2.13) $$\begin{eqnarray}FB_{,t}+L_{a}F_{,t}+\unicode[STIX]{x1D6FE}(q_{u}F)_{,x}=0,\end{eqnarray}$$

where $L_{a}$ is the non-dimensional active layer thickness, which in the following is taken as twice the dimensionless median diameter $d_{50}=d_{50}^{\ast }/D^{\ast }$ (Seminara et al. 1996).

The parameter $\unicode[STIX]{x1D6FE}$ appearing in (2.13) represents the order of magnitude of the ratio between the fluxes per unit width of sediment and flow, corrected for the sediment porosity $p_{s}$ . This small parameter can also be interpreted as the ratio of the characteristic time scales of flow and sediment transport and reads

(2.14) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}=\frac{\sqrt{(s-1)gd_{50}^{\ast }}d_{50}^{\ast }}{(1-p_{s})U^{\ast }D^{\ast }}=O\left(\frac{d_{50}}{C\sqrt{\unicode[STIX]{x1D703}}}\right),\end{eqnarray}$$

where $O(x)$ is used in the following to indicate ‘of the order of $x$ ’.

Integration of (2.13) on the $\unicode[STIX]{x1D719}$ -space yields the classic Exner equation

(2.15) $$\begin{eqnarray}B_{,t}+\unicode[STIX]{x1D6FE}q_{t,x}=0,\end{eqnarray}$$

where (2.8) is used and

(2.16) $$\begin{eqnarray}q_{t}=\int _{-\infty }^{\infty }q_{u}(\unicode[STIX]{x1D719})F(\unicode[STIX]{x1D719})\,\text{d}\unicode[STIX]{x1D719}\end{eqnarray}$$

is the total sediment discharge per unit width.

A closure equation is needed for $q_{u}$ , which is assumed to depend on the fraction grain size through one of the many empirical relationships available in the literature for homogeneous sediment. In particular, we define (Parker, Klingeman & McLean 1982)

(2.17) $$\begin{eqnarray}q_{u}(\unicode[STIX]{x1D719})=\unicode[STIX]{x1D703}^{3/2}G[\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D719})],\quad \unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D719})=\frac{\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D703}_{c}}\frac{d_{50}^{\ast }}{d_{g}^{\ast }}\left(\frac{d_{g}^{\ast }}{d^{\ast }(\unicode[STIX]{x1D719})}\right)^{b},\end{eqnarray}$$

where $\unicode[STIX]{x1D703}_{c}$ is the critical Shields stress for sediment motion.

In table 1 some of the most popular relationships for the function $G(\unicode[STIX]{x1D6F7})$ are summarized, together with the corresponding values of the critical Shields stress. Note that the perturbation approach followed herein requires $G$ to possess continuous derivatives, so that we have assumed the function to be valid in the whole interval $\unicode[STIX]{x1D6F7}\geqslant 1$ , accepting a small error in the limit $\unicode[STIX]{x1D703}\rightarrow \unicode[STIX]{x1D703}_{c}$ , where the functions presented in table 1 are known to become less accurate.

Table 1. The function $G(\unicode[STIX]{x1D6F7})$ and the corresponding value of the critical Shields stress $\unicode[STIX]{x1D703}_{ch}$ for several empirical relationships.

The exponent $b$ in (2.17) controls the hiding effect. Two limiting cases are of relevance: when $b$ vanishes, the ‘equal mobility’ case is recovered, whereby the transport densities of the fractions are all equal to each other and depend on the local value of the Shields stress built upon the geometric mean size. On the contrary, when $b$ is set to unity, the transport density of each fraction depends on the local value of the Shields stress built upon the fraction grain size, as if the bed were composed by that sediment alone. We set $b$ to the value of 0.095, as in Parker (1990), so that even if the same function $G$ is used for all grain sizes, the mobility of each fraction is not identical. Indeed, a small positive value of $b$ corresponds to finer material being slightly more mobile than the coarser one. Hence, under uniform conditions the bed surface is coarser than the bedload and selective transport of surface grains can take place whenever uniform conditions are perturbed. As an example, using the classic Meyer-Peter & Müller (1948) relationship, we obtain

(2.18) $$\begin{eqnarray}q_{ud}=8(\unicode[STIX]{x1D703}_{d}-\unicode[STIX]{x1D709}_{d}\unicode[STIX]{x1D703}_{c})^{3/2},\quad \unicode[STIX]{x1D709}_{d}=\left(\frac{d_{g}^{\ast }}{d^{\ast }}\right)^{1-b},\end{eqnarray}$$

where $q_{ud}$ and $\unicode[STIX]{x1D703}_{d}$ are the transport density and Shields stress built upon the fraction grain size $d^{\ast }$ . The coefficient $\unicode[STIX]{x1D709}_{d}$ in (2.18), which corrects the critical Shields stress $\unicode[STIX]{x1D703}_{c}$ is shown to play the role of the hiding-exposure coefficient of Egiaziaroff (1965), to which it successfully compares.

Finally, the effect of gravity, whereby grains move more easily downhill than uphill, is included in the analysis by reducing the critical Shields stress of an amount proportional to the local slope

(2.19) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{c}=\unicode[STIX]{x1D703}_{ch}-\unicode[STIX]{x1D707}(S-B_{,x}),\end{eqnarray}$$

where $\unicode[STIX]{x1D703}_{ch}$ is the constant value appearing in table 1 and the constant $\unicode[STIX]{x1D707}$ has been set equal to 0.1 after Fredsøe (1974).

3 Linear analysis

A normal mode analysis is developed by expanding the main variables as follows

(3.1) $$\begin{eqnarray}(U,D,B,R)=(1,1,0,R_{0})+\unicode[STIX]{x1D716}\{(U_{1}.D_{1},B_{1},R_{0}R_{1})\exp [\text{i}k(x-wt)]+\text{c.c.}\},\end{eqnarray}$$

where $\unicode[STIX]{x1D716}$ is a small parameter, $k$ and $w$ are the wavenumber and complex wave speed of the perturbation and c.c. stands for complex conjugate. Similarly, any generic derived quantity $f(\unicode[STIX]{x1D719};x,t)$ is expanded as

(3.2) $$\begin{eqnarray}f(\unicode[STIX]{x1D719};x,t)=f_{0}(\unicode[STIX]{x1D719})+\unicode[STIX]{x1D716}\{\,f_{1}(\unicode[STIX]{x1D719})\exp [\text{i}k(x-wt)]+\text{c.c.}\}.\end{eqnarray}$$

Collecting terms with the same power of $\unicode[STIX]{x1D716}$ , the following problems arise.

3.1 Base state: $O(\unicode[STIX]{x1D716}^{0})$

The base state consists, as already stated, of a uniform flow, so that, by definition, we have

(3.3a,b ) $$\begin{eqnarray}\displaystyle T_{t0}=\frac{S}{Fr^{2}}=\frac{1}{C^{2}},\quad T_{n0}=0, & & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle R_{0}=\exp (-\unicode[STIX]{x1D705}C-1)=\frac{r_{0}}{30},\quad r_{0}=n_{r}d_{g0},~d_{g0}=d_{50}, & & \displaystyle\end{eqnarray}$$

where $n_{r}$ is taken as $2.5$ after Engelund & Hansen (1967). Note that the roughness height $R_{0}$ , the bed roughness $r_{0}$ , the conductance coefficient $C$ and the median diameter of the mixture $d_{g0}$ are interrelated by (3.4), so that either of them can be used to fix the grain-to-depth ratio $d_{50}$ .

The Exner equation (2.13) does not provide any additional information, since, under uniform flow conditions, the bed neither experiences aggradation nor degradation. The following relationships hold

(3.5a,b ) $$\begin{eqnarray}\displaystyle q_{u0}(\unicode[STIX]{x1D719})=\unicode[STIX]{x1D703}_{0}^{3/2}G[\unicode[STIX]{x1D6F7}_{0}(\unicode[STIX]{x1D719})],\quad \unicode[STIX]{x1D6F7}_{0}(\unicode[STIX]{x1D719})=\frac{\unicode[STIX]{x1D703}_{0}}{\unicode[STIX]{x1D703}_{c0}}2^{b(\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{m0})}, & & \displaystyle\end{eqnarray}$$

(3.6a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}_{0}=\frac{Fr^{2}}{C^{2}(s-1)d_{50}},\quad \unicode[STIX]{x1D703}_{c0}=\unicode[STIX]{x1D703}_{ch}-\unicode[STIX]{x1D707}\frac{Fr^{2}}{C^{2}}, & & \displaystyle\end{eqnarray}$$

(3.7a-c ) $$\begin{eqnarray}\displaystyle \int _{-\infty }^{\infty }F_{0}(\unicode[STIX]{x1D719})\,\text{d}\unicode[STIX]{x1D719}=1,\quad \unicode[STIX]{x1D719}_{m0}=\int _{-\infty }^{\infty }\unicode[STIX]{x1D719}F_{0}(\unicode[STIX]{x1D719})\,\text{d}\unicode[STIX]{x1D719},\quad \unicode[STIX]{x1D70E}=\int _{-\infty }^{\infty }(\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{m0})^{2}F_{0}(\unicode[STIX]{x1D719})\,\text{d}\unicode[STIX]{x1D719}. & & \displaystyle\end{eqnarray}$$

3.2 Linear level: $O(\unicode[STIX]{x1D716}^{1})$

Substituting the expansion (3.2) into (2.1)–(2.2) and equating terms of $O(\unicode[STIX]{x1D716})$ , the following set of algebraic equations is obtained

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle U_{1}(1-w)+\frac{B_{1}+R_{0}R_{1}+D_{1}}{Fr^{2}}-\frac{\text{i}}{kC^{2}}\left(\frac{T_{t1}}{T_{t0}}-D_{1}\right)-T_{n1}=0, & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle U_{1}+D_{1}(1-w)=0, & \displaystyle\end{eqnarray}$$

where (see appendix A)

(3.10a,b ) $$\begin{eqnarray}\frac{T_{t1}}{T_{t0}}=2U_{1}\quad T_{n1}=\frac{\text{i}k}{C^{2}}\left(NU_{1}-B_{1}-\frac{1}{3}R_{1}\right),\quad N=\frac{1}{3}\left(\unicode[STIX]{x1D705}C+\frac{1}{6}\right).\end{eqnarray}$$

By substituting (3.10) in (3.8) and (3.9), the following algebraic homogeneous system is obtained

(3.11) $$\begin{eqnarray}\left(\begin{array}{@{}cccc@{}}a_{11}-w & a_{12} & a_{13} & a_{14}\\ 1 & 1-w & 0 & 0\end{array}\right)\cdot \left(\begin{array}{@{}c@{}}U_{1}\\ D_{1}\\ B_{1}\\ R_{1}\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\end{array}\right),\end{eqnarray}$$

where

(3.12a,b ) $$\begin{eqnarray}\displaystyle a_{11}=1-\frac{\text{i}}{kC^{2}}(2+k^{2}N),\quad a_{12}=\frac{1}{Fr^{2}}+\frac{\text{i}}{kC^{2}}, & & \displaystyle\end{eqnarray}$$

(3.13a,b ) $$\begin{eqnarray}\displaystyle a_{13}=\frac{1}{Fr^{2}}+\frac{\text{i}k}{C^{2}},\quad a_{14}=\frac{R_{0}}{Fr^{2}}+\frac{\text{i}k}{3C^{2}}. & & \displaystyle\end{eqnarray}$$

Focusing now on the sediment transport model, at the linear level the Exner equation (2.13) yields

(3.14) $$\begin{eqnarray}-wF_{0}B_{1}-wL_{a}F_{1}+\unicode[STIX]{x1D6FE}(F_{0}q_{u1}+q_{u0}F_{1})=0,\end{eqnarray}$$

where

(3.15) $$\begin{eqnarray}\displaystyle q_{u1}=\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D719})\frac{T_{t1}}{T_{t0}}+\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D719})\left[\frac{T_{t1}}{T_{t0}}-\frac{\unicode[STIX]{x1D703}_{c1}}{\unicode[STIX]{x1D703}_{c0}}-(1-b)\frac{d_{g1}}{d_{g0}}\right], & & \displaystyle\end{eqnarray}$$

(3.16a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D719})={\textstyle \frac{3}{2}}\unicode[STIX]{x1D703}_{0}^{3/2}G[\unicode[STIX]{x1D6F7}_{0}(\unicode[STIX]{x1D719})],\quad \unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D719})=\unicode[STIX]{x1D703}_{0}^{3/2}\unicode[STIX]{x1D6F7}_{0}(\unicode[STIX]{x1D719})G^{\prime }[\unicode[STIX]{x1D6F7}_{0}(\unicode[STIX]{x1D719})], & & \displaystyle\end{eqnarray}$$

(3.17a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}_{c1}=\text{i}k\unicode[STIX]{x1D707}B_{1},\quad \frac{d_{g1}}{d_{g0}}=\frac{r_{1}}{r_{0}}=R_{1}=-\ln (2)\unicode[STIX]{x1D719}_{m1}, & & \displaystyle\end{eqnarray}$$

(3.18a,b ) $$\begin{eqnarray}\displaystyle \int _{-\infty }^{\infty }F_{1}(\unicode[STIX]{x1D719})\,\text{d}\unicode[STIX]{x1D719}=0,\quad \unicode[STIX]{x1D719}_{m1}=\int _{-\infty }^{\infty }\unicode[STIX]{x1D719}F_{1}(\unicode[STIX]{x1D719})\,\text{d}\unicode[STIX]{x1D719}, & & \displaystyle\end{eqnarray}$$

and $G^{\prime }$ denotes the derivative of $G$ with respect to $\unicode[STIX]{x1D6F7}$ .

Equation (3.14) can be integrated in $\unicode[STIX]{x1D719}$ leading to

(3.19) $$\begin{eqnarray}-wB_{1}+\unicode[STIX]{x1D6FE}q_{t1}=0,\end{eqnarray}$$

where

(3.20) $$\begin{eqnarray}q_{t1}=\int _{-\infty }^{\infty }[F_{0}(\unicode[STIX]{x1D719})q_{u1}(\unicode[STIX]{x1D719})+q_{u0}(\unicode[STIX]{x1D719})F_{1}(\unicode[STIX]{x1D719})]\,\text{d}\unicode[STIX]{x1D719}\end{eqnarray}$$

is the perturbation of the total sediment discharge per unit width.

Moreover, equation (3.19) can be substituted in (3.14) leading to

(3.21) $$\begin{eqnarray}-wL_{a}F_{1}(\unicode[STIX]{x1D719})+\unicode[STIX]{x1D6FE}[F_{0}(\unicode[STIX]{x1D719})q_{u1}(\unicode[STIX]{x1D719})-F_{0}(\unicode[STIX]{x1D719})q_{t1}+q_{u0}(\unicode[STIX]{x1D719})F_{1}(\unicode[STIX]{x1D719})]=0,\end{eqnarray}$$

which, together with (3.19), defines the linear system that controls the bed evolution, in terms of bed elevation and composition. Assuming that the size distribution of the mixture can be approximated by means of $n$ fractions, the continuous probability density function $F(\unicode[STIX]{x1D6F7};x,t)$ is formally replaced by its discrete counterpart: the probability mass function (in short, the size fraction) $F(\unicode[STIX]{x1D6F7}^{i};x,t)$ . Note that, in this case, equations (3.19) and (3.21) provide a total of $n$ equations in $n$ unknown variables: the amplitude of the bed perturbation $B_{1}$ plus the $n-1$ perturbations of the size fraction $F_{1}(\unicode[STIX]{x1D719}^{i})$ , the last size fraction $F_{1}(\unicode[STIX]{x1D719}^{n})$ being determined by the condition that the summation of all $F_{1}$ must vanish.

Let us restrict our attention, from now on, to the simplest possible case of heterogeneous sediment: a bimodal mixture composed by two species ( $a$ and $b$ ) in equal parts. Although most of the following considerations also apply to the more general case of $n$ fractions, this choice will allow for discussing the behaviour of a single sorting wave instead of that of a family of $n-1$ eigenvalues, which indeed behave quite similarly (Stecca et al. 2014). We set

(3.22a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}^{a}=\unicode[STIX]{x1D719}_{m0}+\unicode[STIX]{x1D70E},\quad \unicode[STIX]{x1D719}^{b}=\unicode[STIX]{x1D719}_{m0}-\unicode[STIX]{x1D70E},\quad F_{0}^{a}=F_{0}^{b}={\textstyle \frac{1}{2}}, & & \displaystyle\end{eqnarray}$$

(3.23a-c ) $$\begin{eqnarray}\displaystyle F_{1}^{a}+F_{1}^{b}=0,\quad \unicode[STIX]{x1D719}_{m1}=2\unicode[STIX]{x1D70E}F_{1}^{a},\quad R_{1}=-2\unicode[STIX]{x1D70E}\ln (2)F_{1}^{a}, & & \displaystyle\end{eqnarray}$$

where the superscripts $^{a}$ and $^{b}$ refer to the finer and the coarser fraction, respectively. Note that (3.23) shows that the amplitude of the perturbations of the roughness height $R_{1}$ is proportional to that of the coarser size fraction $F_{1}^{b}=-F_{1}^{a}$ , so that the former is maximum where the areal concentration of the coarser sediment is higher.

Using (3.19) and (3.21), the following algebraic homogeneous system is obtained:

(3.24) $$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\displaystyle -wB_{1}+\unicode[STIX]{x1D6FE}(q_{u1}^{+}+2q_{u0}^{-}F_{1}^{a})=0,\\ \displaystyle -wF_{1}^{a}-\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2L_{a}}(q_{u1}^{-}+2q_{u0}^{+}F_{1}^{a})=0,\end{array}\right\}\end{eqnarray}$$

where

(3.25a,b ) $$\begin{eqnarray}q_{u0}^{\pm }={\textstyle \frac{1}{2}}(q_{u0}^{a}\pm q_{u0}^{b}),\quad q_{u1}^{\pm }={\textstyle \frac{1}{2}}(q_{u1}^{a}\pm q_{u1}^{b}).\end{eqnarray}$$

Moreover, making use of (3.15)–(3.17) and (3.23) the system (3.24) can be rewritten as

(3.26) $$\begin{eqnarray}\left(\begin{array}{@{}cccc@{}}\unicode[STIX]{x1D6FE}a_{31} & 0 & \unicode[STIX]{x1D6FE}\text{i}ka_{33}-w & \unicode[STIX]{x1D6FE}a_{34}\\ \unicode[STIX]{x1D6E4}a_{41} & 0 & \unicode[STIX]{x1D6E4}\text{i}ka_{43} & \unicode[STIX]{x1D6E4}a_{44}-w\end{array}\right)\cdot \left(\begin{array}{@{}c@{}}U_{1}\\ D_{1}\\ B_{1}\\ R_{1}\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\end{array}\right),\end{eqnarray}$$

where

(3.27) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}a_{31}=2(\unicode[STIX]{x1D6FC}^{+}+\unicode[STIX]{x1D6FD}^{+}),\quad a_{33}=-\unicode[STIX]{x1D6FD}^{+}\displaystyle \frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D703}_{c0}},\quad a_{34}=-\displaystyle \frac{2}{3}\frac{\unicode[STIX]{x1D6FC}^{-}}{\unicode[STIX]{x1D70E}\ln (2)}-\unicode[STIX]{x1D6FD}^{+}(1-b),\\ a_{41}=-2\unicode[STIX]{x1D70E}\left(\unicode[STIX]{x1D6FC}^{-}+\unicode[STIX]{x1D6FD}^{-}\right),\quad a_{43}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FD}^{-}\displaystyle \frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D703}_{c0}},\quad a_{44}=\displaystyle \frac{2}{3}\frac{\unicode[STIX]{x1D6FC}^{+}}{\ln (2)}+\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FD}^{-}(1-b),\end{array}\right\}\end{eqnarray}$$

and, as in (3.25), $\unicode[STIX]{x1D6FC}^{\pm }$ and $\unicode[STIX]{x1D6FD}^{\pm }$ are the halved sums and differences of the values of the functions $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ evaluated at $\unicode[STIX]{x1D719}^{a}$ and $\unicode[STIX]{x1D719}^{b}$ .

The parameter $\unicode[STIX]{x1D6E4}$ , which appears in the last row of (3.26), is equal to

(3.28) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}=\frac{\unicode[STIX]{x1D6FE}\ln (2)}{L_{a}}=O\left(\frac{1}{C\sqrt{\unicode[STIX]{x1D703}_{0}}}\right)\end{eqnarray}$$

and, in strict analogy with $\unicode[STIX]{x1D6FE}$ , it can be considered as the ratio of the characteristic time scales of flow and sorting, an assumption the validity of which will be confirmed later on. Although large with respect to $\unicode[STIX]{x1D6FE}$ , since $L_{a}$ is proportional to the (small) grain-to-depth ratio $d_{50}$ , $\unicode[STIX]{x1D6E4}$ can still be considered as a small parameter. Indeed, the smallness of $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6E4}$ ultimately ensures that both the bed and the sorting processes do not evolve faster than the flow, a condition that would lead to quite unrealistic and unphysical results. Hence, we have

(3.29a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}\ll \unicode[STIX]{x1D6E4}\ll 1,\quad \unicode[STIX]{x1D6FF}=\frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D6E4}}\ll 1,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}$ represents the ratio of the characteristic time scales of topographic and sorting waves.

Combining together (3.11) and (3.26), the fully coupled morphodynamic system for a two-grain sediment mixture is obtained, in the form

(3.30) $$\begin{eqnarray}\left(\begin{array}{@{}cccc@{}}a_{11}-w & a_{12} & a_{13} & a_{14}\\ 1 & 1-w & 0 & 0\\ \unicode[STIX]{x1D6FE}a_{31} & 0 & \unicode[STIX]{x1D6FE}\text{i}ka_{33}-w & \unicode[STIX]{x1D6FE}a_{34}\\ \unicode[STIX]{x1D6E4}a_{41} & 0 & \unicode[STIX]{x1D6E4}\text{i}ka_{43} & \unicode[STIX]{x1D6E4}a_{44}-w\end{array}\right)\cdot \left(\begin{array}{@{}c@{}}U_{1}\\ D_{1}\\ B_{1}\\ R_{1}\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\\ 0\\ 0\end{array}\right)\end{eqnarray}$$

or, in compact form,

(3.31) $$\begin{eqnarray}(\boldsymbol{A}-w\unicode[STIX]{x1D644})\boldsymbol{\cdot }\boldsymbol{x}=\{0\},\end{eqnarray}$$

where $\boldsymbol{A}$ and $\boldsymbol{x}$ have obvious definitions and $\unicode[STIX]{x1D644}$ is the identity matrix. This equation reveals itself as a classic eigenvalue problem, whereby the system admits a non-trivial solution only for those specific values of $w$ , the eigenvalues, for which

(3.32) $$\begin{eqnarray}\text{det}(\boldsymbol{A}-w\unicode[STIX]{x1D644})=0,\end{eqnarray}$$

where $\text{det}(\cdot )$ stands for the determinant of the matrix.

The analysis of the eigenvalues of (3.32) and of the corresponding eigenvectors provides the required information on the stability of the system. In particular, for any given set of values of the parameters of the problem, the eigenvalue determines the growth rate and the celerity (i.e. the wave speed) of the perturbations, which read, respectively

(3.33a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}=kw^{i},\quad \unicode[STIX]{x1D714}=w^{r},\end{eqnarray}$$

where the superscripts $^{r,i}$ stand for the real and the imaginary part of the complex quantity, respectively.

A positive (negative) growth rate implies instability (stability) of the base state, whereas a positive (negative) celerity implies downstream (upstream) propagation.

The system (3.30) represents the core of the present linear stability analysis and, therefore, its solutions will be discussed in detail in the following. It is worth noting that this algebraic eigenvalue problem provides four separate eigenvalues, each of them describing a different physical process and each of them being in principle able to induce an instability of the base uniform flow. In the next section, the nature of these eigenvalues will be thoroughly investigated and it will be shown how they can be associated with different patterns. Indeed, the eigenvectors associated with each eigenvalue involve perturbations of the same quantities (velocity, flow depth, topography and roughness) but, as in cooking, the same ingredients, combined in different ways, can provide quite different results. In particular, it will be shown that one of the eigenvalues can be associated with an instability of the bedload-sheet kind: a sorting wave which manifests itself as a streamwise periodic perturbation of the surface composition of the bed, travelling downstream with only a minimal change in bed elevation.

4 Flow, bed and sorting eigenvalues

In this section we will examine the characteristics of the four eigenvalues of the complete system by isolating simpler eigenproblems, already contained in (3.30), which will provide some insight on the physical processes described by each of the eigenvalues of the original problem. Moreover, we will show how the character (and thus the labelling) of the eigenvalues remains valid as we zoom out from the simpler towards the more complex problem. When appropriate, we will also take advantage of the presence of the small parameters introduced in the previous section (namely $\unicode[STIX]{x1D6FE}$ , $\unicode[STIX]{x1D6E4}$ and $\unicode[STIX]{x1D6FF}$ ), by suitably expanding the eigenvalues. The quasi-steady approximation and the limiting case of weak sorting will also be considered and discussed.

Let us first rewrite (3.30) as

(4.1) $$\begin{eqnarray}\left(\begin{array}{@{}cccc@{}}a_{11}-\unicode[STIX]{x1D6FE}W & a_{12} & a_{13} & a_{14}\\ 1 & 1-\unicode[STIX]{x1D6FE}W & 0 & 0\\ a_{31} & 0 & \text{i}ka_{33}-W & a_{34}\\ a_{41} & 0 & \text{i}ka_{43} & a_{44}-\unicode[STIX]{x1D6FF}W\end{array}\right)\cdot \left(\begin{array}{@{}c@{}}U_{1}\\ D_{1}\\ B_{1}\\ R_{1}\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\\ 0\\ 0\end{array}\right),\end{eqnarray}$$

where a new complex wave speed $W=w/\unicode[STIX]{x1D6FE}$ has been introduced, thus adopting the slow time scale of sediment transport instead of the fast time scale of the flow.

4.1 The flow eigenvalues

We start from the simple hydrodynamic problem, whereby the bed is assumed to be flat, uniformly rough and unerodible. The latter conditions can be easily imposed by setting $B_{1}$ and $R_{1}$ to zero (flat bed of uniform roughness) and by dropping the Exner equations (unerodible bed).

The system reduces to the rank-2 submatrix in the upper-left corner of (4.1) and the corresponding characteristic polynomial takes the form

(4.2) $$\begin{eqnarray}(\unicode[STIX]{x1D6FE}W_{f})^{2}-a_{0}(\unicode[STIX]{x1D6FE}W_{f})+b_{0}=0,\end{eqnarray}$$

where

(4.3a,b ) $$\begin{eqnarray}a_{0}=1+a_{11},\quad b_{0}=a_{11}-a_{12}.\end{eqnarray}$$

Two complex eigenvalues $w_{f}^{\pm }=\unicode[STIX]{x1D6FE}W_{f}^{\pm }$ are obtained as the roots of the quadratic eigenrelation (4.2) and the associated growth rates and celerities are readily obtained using (3.33). In particular, in the short-wave limit ( $k\rightarrow \infty$ ) and neglecting the normal stress term $T_{n1}$ , we recover the classic result

(4.4a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{f}^{\pm }C^{2}\rightarrow -1\pm \frac{Fr}{2},\quad \unicode[STIX]{x1D714}_{f}^{\pm }\rightarrow 1\pm \frac{1}{Fr},\end{eqnarray}$$

where the ‘slow’ eigenvalue $w_{f}^{-}$ (i.e. the one characterized by the smallest celerity) is always stable and propagates upstream (downstream) in the subcritical (supercritical) regime, whereas the ‘fast’ eigenvalue $w_{f}^{+}$ always propagates downstream and it becomes unstable for $Fr>2$ , leading to the formation of roll waves.

The growth rate of the fast eigenvalue being bounded with a maximum for the shortest waves, the solution described by (4.4) is regular in the sense of Birkhoff (1954) and thus the problem is mathematically well posed. However, the physical interpretation of the stability results is not particularly appealing, because of the absence of a suitable cutoff mechanism. As soon as the Froude number exceeds the critical threshold, all modes in the short-wave range are equally unstable. In this regard, the normal stress perturbation $T_{n1}$ plays a fundamental role, since it provides a diffusive damping of the growth rate in the short-wave range (Needham & Merkin 1984). Including $T_{n1}$ , the limits of the growth rate and celerity of the fast eigenvalue in the short-wave range become, respectively

(4.5a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{f}^{+}\rightarrow -\frac{C^{2}}{NFr^{2}},\quad \unicode[STIX]{x1D714}_{f}^{+}\rightarrow 1.\end{eqnarray}$$

Note that the normal stress term only inhibits the growth of the shortest modes, whereas some intermediate modes are still unstable if the threshold is exceeded. This, in turn, provides a wavelength selection mechanism.

4.2 The bed eigenvalue

We now consider the rank-3 submatrix in the upper-left corner of (4.1), which defines the classic morphodynamic problem that controls bed form instability for uniform sediment, thus implying constant roughness. Hence, we can set $R_{1}$ to zero, dropping only the Exner equation for the fraction. In this case

(4.6a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}^{+}=\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D719}_{m0})=\unicode[STIX]{x1D6FC}_{m},\quad \unicode[STIX]{x1D6FD}^{+}=\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D719}_{m0})=\unicode[STIX]{x1D6FD}_{m}.\end{eqnarray}$$

The characteristic polynomial takes the form

(4.7) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}^{2}W^{3}-\unicode[STIX]{x1D6FE}W^{2}(a_{0}+\unicode[STIX]{x1D6FE}a_{1})+W(b_{0}+\unicode[STIX]{x1D6FE}b_{1})-c_{0}=0,\end{eqnarray}$$

where

(4.8a-c ) $$\begin{eqnarray}a_{1}=\text{i}ka_{33},\quad b_{1}=a_{0}a_{1}-a_{13}a_{31},\quad c_{0}=a_{1}b_{0}-a_{13}a_{31}.\end{eqnarray}$$

The presence of the small parameter $\unicode[STIX]{x1D6FE}$ in (4.7) suggests the use of perturbation theory to find approximate solutions for the roots of the cubic polynomial. However, for vanishing $\unicode[STIX]{x1D6FE}$ , the degree of the polynomial is reduced to unity so that two roots are missing. Such an abrupt change in the character of the solution is a symptom that we are dealing with a singular perturbation problem, as also hinted by the fact that the small parameter multiplies the higher degree term of the polynomial.

We then expand the roots of (4.7) as

(4.9) $$\begin{eqnarray}W=\unicode[STIX]{x1D6FE}^{-1}W_{-1}+W_{0}+O(\unicode[STIX]{x1D6FE}),\end{eqnarray}$$

where the first term of the expansion takes care of the singularity so that the approximate roots converge in the limit $\unicode[STIX]{x1D6FE}\rightarrow 0$ . Note that $O(x)$ is used in the following to indicate both ‘of the order of’ $x$ and to summarize all the least-significant terms in a series, as above. When $x$ is a complex number, the order of its modulus is considered.

Substituting (4.9) in (4.7) and collecting terms at $O(\unicode[STIX]{x1D6FE}^{-1})$ we obtain

(4.10a,b ) $$\begin{eqnarray}W_{-1}(W_{-1}^{2}-a_{0}W_{-1}+b_{0})=0,\quad W_{-1}=w_{f}^{+},w_{f}^{-},0\end{eqnarray}$$

having recognized that the roots of the quadratic polynomial between parentheses are none other than the flow eigenvalues over a flat bed $w_{f}^{\pm }$ . Moreover, at $O(\unicode[STIX]{x1D6FE}^{0})$ we get:

(4.11) $$\begin{eqnarray}W_{0}=\frac{a_{1}W_{-1}^{2}-b_{1}W_{-1}+c_{0}}{3W_{-1}^{2}-2a_{0}W_{-1}+b_{0}}.\end{eqnarray}$$

The approximate roots of (4.7) are then

(4.12a,b ) $$\begin{eqnarray}w_{F}^{\pm }=w_{f}^{\pm }+\unicode[STIX]{x1D6FE}\frac{a_{13}a_{31}(w_{f}^{\pm }-1)}{w_{f}^{\pm }(1+a_{11})-2(a_{11}-a_{12})},\quad w_{b}=\unicode[STIX]{x1D6FE}\left(\text{i}ka_{33}-\frac{a_{13}a_{31}}{a_{11}-a_{12}}\right).\end{eqnarray}$$

The above procedure sheds some light on the nature of the eigenvalues of the morphodynamic problem of flow over a mobile bed composed by a uniform sediment: (i) two of the three eigenvalues can still be labelled as flow eigenvalues, since they are found to be equal to the flow eigenvalues over a flat bed plus a small $O(\unicode[STIX]{x1D6FE})$ correction, which accounts for the presence of the mobile bed; (ii) the third eigenvalue is new, is itself of $O(\unicode[STIX]{x1D6FE})$ and can be labelled as the ‘bed’ eigenvalue, since it appears only if the bed is allowed to be modified by the flow.

In particular, the fast flow eigenvalue $w_{F}^{+}$ still describes the roll-wave instability discussed in § 4.1, which now also involves a modification of the incoherent bed. After some algebra, the bed eigenvalue eventually turns out to be

(4.13) $$\begin{eqnarray}w_{b}=\unicode[STIX]{x1D6FE}\left[(\unicode[STIX]{x1D6FC}_{m}+\unicode[STIX]{x1D6FD}_{m})\frac{T_{B}}{T_{t0}}-\unicode[STIX]{x1D6FD}_{m}\frac{\text{i}k\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D703}_{c0}}\right],\end{eqnarray}$$

where

(4.14) $$\begin{eqnarray}T_{B}=-2T_{t0}\frac{a_{13}}{a_{11}-a_{12}}\end{eqnarray}$$

is the part of the perturbation of the bed shear stress $T_{t1}$ proportional to $B_{1}$ .

The celerity and growth rate of the bed wave are then equal to

(4.15a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{b}=\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D6FC}_{m}+\unicode[STIX]{x1D6FD}_{m})\frac{T_{B}^{r}}{T_{t0}},\quad \unicode[STIX]{x1D6FA}_{b}=\unicode[STIX]{x1D6FE}\left[(\unicode[STIX]{x1D6FC}_{m}+\unicode[STIX]{x1D6FD}_{m})\frac{kT_{B}^{i}}{T_{t0}}-k^{2}\unicode[STIX]{x1D6FD}_{m}\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D703}_{c0}}\right].\end{eqnarray}$$

Moreover, ignoring for simplicity the effect of the normal stress $T_{n1}$ we obtain

(4.16) $$\begin{eqnarray}\frac{T_{B}}{T_{t0}}=-C^{2}\frac{2k^{2}C^{2}(Fr^{2}-1)+6\text{i}kFr^{2}}{k^{2}C^{4}(Fr^{2}-1)^{2}+9Fr^{4}}.\end{eqnarray}$$

By combining (4.15) and (4.16), some well-known results are recovered: (i) the celerity of the bed wave is of $O(\unicode[STIX]{x1D6FE})$ (i.e. the bed wave is much slower than a flow wave); (ii) the bed wave propagates downstream for $Fr<1$ and upstream for $Fr>1$ ; (iii) the base uniform flow is always stable with respect to perturbations of the bed elevation, since the imaginary part of the bed shear stress is consistently negative (i.e. the stress perturbation lags the bed); (iv) gravity (i.e. the streamwise bed slope) acts as a stabilizing effect, the damping increasing with the wavenumber squared.

Moreover, in the short-wave ( $k\rightarrow \infty$ ) limit

(4.17a,b ) $$\begin{eqnarray}\frac{T_{B}^{r}}{T_{t0}}\rightarrow \frac{2}{(1-Fr^{2})},\quad \frac{kT_{B}^{i}}{T_{t0}}\rightarrow -\frac{6Fr^{2}}{C^{2}(Fr^{2}-1)^{2}}\end{eqnarray}$$

so that resonance occurs when $Fr\rightarrow 1$ , whereby

(4.18a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{b}\rightarrow \pm \infty \quad \text{as}~Fr\rightarrow 1^{\mp },\;\unicode[STIX]{x1D6FA}_{b}\rightarrow -\infty .\end{eqnarray}$$

We remark that this resonance appears only when approximate roots of (4.7) are sought, because the expansion (4.9) is not uniformly valid in the transcritical $Fr\simeq 1$ region as $k\rightarrow \infty$ . Indeed, in the complete analysis of Lyn & Altinakar (2002), the celerities associated with the three eigenvalues are always two positive and one negative and there is no sign of this spurious resonance, the physical meaning of which will become clear in the next subsection.

Note also that if the normal stress $T_{n1}$ is included in the analysis, its damping, stabilizing effect modifies the behaviour of the solution in the short-wave range, as in the roll-wave example cited before

(4.19a,b ) $$\begin{eqnarray}\frac{T_{B}^{r}}{T_{t0}}\rightarrow \frac{2}{N},\quad \frac{kT_{B}^{i}}{T_{t0}}\rightarrow -\frac{2C^{2}}{N^{2}Fr^{2}}(Fr^{2}-1+N).\end{eqnarray}$$

Hence, resonance disappears and the growth rate remains consistently negative, the coefficient $N$ being typically larger than 1.

By adopting a shallow-water flow model, the linear morphodynamic analysis tends to be limited in application to processes occurring over distances larger than the flow depth. As a result, the dynamics of dunes and antidunes, the wavelength of which scales with the flow depth, cannot be handled satisfactorily (Lanzoni et al. 2006). In particular, the only bed waves that are found to be linearly unstable with uniform sediment are those associated with the formation of roll waves (Balmforth & Vakil 2012), whereby the process driving the instability has to be sought in the interactions between the flow and the free surface more than in the interactions of the flow with the erodible bed.

4.3 The quasi-steady case

One of the most popular simplifications adopted in the study of morphodynamic problems is the quasi-steady approximation, whereby times derivatives are dropped in the flow equations. The rationale behind this hypothesis stands in the smallness of the parameter $\unicode[STIX]{x1D6FE}$ , which implies that the bed dynamics evolves on a much slower time scale with respect to flow. As a consequence, the flow is assumed to adapt instantaneously to the modifications of the bed and the analysis can be split in two parts. Firstly, the linear response of the flow to a steady perturbation of the bed is determined using the flow equations and, secondly, the (only) eigenvalue of the problem is found by substituting the latter into the Exner equation.

If we enforce the quasi-steady approximation on the morphodynamic system by dropping time derivatives in the flow equations we obtain the so called ‘decoupled’ system, which reads

(4.20) $$\begin{eqnarray}\left(\begin{array}{@{}ccc@{}}a_{11} & a_{12} & a_{13}\\ 1 & 1 & 0\\ a_{31} & 0 & \text{i}ka_{33}-W\end{array}\right)\cdot \left(\begin{array}{@{}c@{}}U_{1}\\ D_{1}\\ B_{1}\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\\ 0\end{array}\right).\end{eqnarray}$$

Using the first two rows of (4.20), $U_{1}$ and $D_{1}$ can be expressed in terms of $B_{1}$

(4.21a,b ) $$\begin{eqnarray}U_{1}=-D_{1}=U_{B}B_{1},\quad U_{B}=-\frac{a_{13}}{a_{11}-a_{12}}=\frac{T_{B}}{2T_{t0}}.\end{eqnarray}$$

Substituting (4.21) into the third row of (4.20) we obtain the decoupled bed eigenvalue

(4.22) $$\begin{eqnarray}w_{bd}=\unicode[STIX]{x1D6FE}W=\unicode[STIX]{x1D6FE}(\text{i}ka_{33}+a_{31}U_{B})=w_{b},\end{eqnarray}$$

which is identical to the bed eigenvalue obtained in the previous subsection. As far as the bed eigenvalue is concerned, the two procedures are then formally equivalent. However, the decoupling procedure rules out any mode of instability associated with the flow: only the bed eigenvalue is present and flow instability cannot any more drive the formation of bed forms, as in the roll-wave example previously outlined. Indeed, it must be pointed out that the use of the quasi-steady approximation does not cancel out the flow–bed feedback mechanism that drives the process of instability: a growth (decay) of the bed perturbation is observed whenever the deformation of the bed induces a stress field able to increment (decrement) the amplitude of the perturbation itself.

The spurious resonance in the transcritical region is still present, but is amenable in this context of a more physical explanation: the system being decoupled, the slow bed disturbance is felt by the flow as an external forcing which excites a ‘natural’ frequency of the flow itself, thus leading the system to resonate. We remark once more that the resonance disappears once the complete coupled morphodynamic problem is solved (Lyn & Altinakar 2002).

Moving back to the case of heterogeneous sediment and making use of the quasi-steady hypothesis, $U_{1}$ and $D_{1}$ can be expressed as linear combinations of $B_{1}$ and $R_{1}$

(4.23a,b ) $$\begin{eqnarray}U_{1}=-D_{1}=U_{B}B_{1}+U_{R}R_{1},\quad U_{R}=-\frac{a_{14}}{a_{11}-a_{12}}=\frac{T_{R}}{2T_{t0}}.\end{eqnarray}$$

Moreover, the system (3.26) can then be rewritten as

(4.24) $$\begin{eqnarray}\left(\begin{array}{@{}cc@{}}\text{i}ka_{33}+a_{31}U_{B}-W & a_{34}+a_{31}U_{R}\\ \text{i}ka_{43}+a_{41}U_{B} & a_{44}+a_{41}U_{R}-\unicode[STIX]{x1D6FF}W\end{array}\right)\cdot \left(\begin{array}{@{}c@{}}B_{1}\\ R_{1}\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\end{array}\right)\end{eqnarray}$$

revealing itself as a new eigenvalue problem, the solution of which provides two eigenvalues, one to be associated with a bed wave, the other with a sorting wave. It must be pointed out that each wave involves a perturbation of both the bed elevation and the roughness. Which of the two defines the character of the wave will eventually become clear in the next subsection.

4.4 The sorting eigenvalue

The characteristic polynomial of (4.24) takes the form

(4.25) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}W^{2}-(T_{0}+\unicode[STIX]{x1D6FF}T_{1})W+D_{0}=0,\end{eqnarray}$$

where

(4.26a,b ) $$\begin{eqnarray}\displaystyle T_{0}=a_{44}+a_{41}U_{R},\quad T_{1}=\text{i}ka_{33}+a_{31}U_{B}, & & \displaystyle\end{eqnarray}$$

(4.27) $$\begin{eqnarray}\displaystyle D_{0}=T_{0}T_{1}-(\text{i}ka_{43}+a_{41}U_{B})(a_{34}+a_{31}U_{R}). & & \displaystyle\end{eqnarray}$$

The presence of the small parameter $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D6E4}$ suggests that it is useful to expand the roots of (4.25) as

(4.28) $$\begin{eqnarray}W=\unicode[STIX]{x1D6FF}^{-1}W_{-1}+W_{0}+O(\unicode[STIX]{x1D6FF}),\end{eqnarray}$$

where, as before, the first term of the expansion is needed because (4.25) is easily recognized as a singular perturbation problem as $\unicode[STIX]{x1D6FF}\rightarrow 0$ .

Collecting terms at various orders we obtain

(4.29a,b ) $$\begin{eqnarray}W_{-1}=0,T_{0},\quad W_{0}=\frac{D_{0}}{T_{0}},T_{1}-\frac{D_{0}}{T_{0}}\end{eqnarray}$$

and

(4.30a,b ) $$\begin{eqnarray}w_{B}=w_{b}-w_{sb},\quad w_{S}=w_{s}+w_{sb},\end{eqnarray}$$

where

(4.31a,b ) $$\begin{eqnarray}w_{s}=\unicode[STIX]{x1D6E4}(a_{44}+a_{41}U_{R}),\quad w_{sb}=\unicode[STIX]{x1D6FE}\displaystyle \frac{(\text{i}ka_{43}+a_{41}U_{B})(a_{34}+a_{31}U_{R})}{a_{44}+a_{41}U_{R}}.\end{eqnarray}$$

The approximate solution (4.30) provides some insight into the nature of the eigenvalues of the quasi-steady morphodynamic problem for heterogeneous sediment: (i) one of the two eigenvalues can still be labelled as the bed eigenvalue, since it is found equal to the one obtained for the case of uniform sediment $w_{b}$ plus a $O(\unicode[STIX]{x1D6FE})$ correction that represents the reciprocal interactions between sorting and bed waves; (ii) the same correction, with opposite sign, is present in the second eigenvalue, which is new and can be labelled as the ‘sorting’ eigenvalue, since it appears only if the bed is composed by a sediment mixture; (iii) the leading-order term of the sorting eigenvalue is of $O(\unicode[STIX]{x1D6E4})$ . Sorting waves are then found to propagate much faster than bed waves (Ribberink 1980), the order of magnitude of the ratio of the two celerities (the parameter $\unicode[STIX]{x1D6FF}$ above) being related to the non-dimensional thickness of the active layer $L_{a}$ and so, in turn, on the grain-to-depth ratio $d_{50}$ (Stecca et al. 2014).

The above considerations are used to estimate of the limits of applicability of the quasi-steady hypothesis, since not only $\unicode[STIX]{x1D6FE}$ but also $\unicode[STIX]{x1D6E4}$ must be small for this approximation to be valid. Using (2.14) and (3.28) we obtain that both conditions hold if

(4.32) $$\begin{eqnarray}d_{50}\ll 1\ll C\sqrt{\unicode[STIX]{x1D703}_{0}}.\end{eqnarray}$$

This inequality also ensures that sorting and bed perturbations do not travel faster than flow perturbations. Moreover, equation (4.32) suggests that the analysis is likely to fail for very coarse mixtures or very shallow flows (large $d_{50}$ and small $C$ ) and for values of the Shields stress close to its critical threshold. It must be pointed out that this result is totally independent of the choice of the length scale for $L_{a}^{\ast }$ . Whether the grain size, as in the present case, or the bed form amplitude, which itself scales with the flow depth, is used to express the active layer thickness, this only alters the relative celerity of sorting and bed waves, which both became of $O(\unicode[STIX]{x1D6FE})$ in the latter case. Moreover, the idea that a minimum active layer thickness must be set in order to limit the speed of the sediment waves (Sieben 1997) seems unfounded, the only limit being that the active layer cannot obviously be thinner than the median diameter of the mixture.

Another important information on sorting waves can be extracted from the approximated eigenvectors associated with (4.30). We have

(4.33) $$\begin{eqnarray}\displaystyle & \displaystyle B_{s}=\frac{1}{R_{s}}=\frac{w_{s}-\unicode[STIX]{x1D6E4}(a_{44}+a_{41}U_{R})}{\unicode[STIX]{x1D6E4}(\text{i}ka_{43}+a_{41}U_{B})}=\unicode[STIX]{x1D6FF}\frac{(a_{34}+a_{31}U_{R})}{(a_{44}+a_{41}U_{R})}, & \displaystyle\end{eqnarray}$$

(4.34) $$\begin{eqnarray}\displaystyle & \displaystyle B_{b}=\frac{1}{R_{b}}=\frac{\unicode[STIX]{x1D6FE}(a_{34}+a_{31}U_{R})}{w_{b}-\unicode[STIX]{x1D6FE}(\text{i}ka_{33}+a_{31}U_{B})}=-\frac{(a_{44}+a_{41}U_{R})}{(\text{i}ka_{43}+a_{41}U_{B})}, & \displaystyle\end{eqnarray}$$

where $B_{s}$ and $B_{b}$ are the amplitudes of the bed perturbation corresponding to a unitary perturbation of the roughness height for the sorting and the bed wave, respectively. Similarly, $R_{s}$ and $R_{b}$ are the amplitudes of the roughness height perturbation for a unitary bed perturbation.

The ratio of (4.33) and (4.34) turns out to be of $O(\unicode[STIX]{x1D6FF})$ , revealing that the sorting wave is associated with a negligible amplitude of the bed perturbation, whereas the bed wave is associated with a negligible amplitude of the roughness height perturbation and, ultimately, of the perturbation of the areal concentration of the sediment fraction. Hence, in the present one-dimensional context, the two eigenvalues well describe dune and bedload-sheet instabilities, whereby bed elevation dominates over sorting in the former and sorting rules over topography in the latter.

The stability picture emerging from the analysis of the multi-size case is indeed much richer than for the case of uniform sediment: (i) sorting and bed waves can both be stable and no instability will develop; (ii) either of them can grow, with the other decaying; (iii) they can both be linearly unstable, selecting different wavelengths and with different growth rates, the competition among them being controlled by nonlinear interactions; (iv) although perturbations of the bed elevation and of the size fraction develop side by side in both instabilities, sorting and bed waves produce quite different patterns in terms of the relative amplitudes of bed and roughness perturbations.

As already discussed at the end of § 4.2, by adopting a shallow-water flow model we have chosen to privilege the benefit of an algebraic eigenproblem, which allows for a deep insight into the characteristics of bed and sorting waves, to the detriment of an accurate representation of the phase lag between shear stress and the bed, which controls dune instability.

4.5 The weak-sorting case

Let us finally introduce the weak-sorting case, whereby the two fractions are assumed to be infinitely close one to the other. The idea is to observe how the eigenvalue problem modifies as the sediment mixture composing the bed becomes less and less heterogeneous and formally corresponds to taking the limit of the eigenvalues as $\unicode[STIX]{x1D70E}\rightarrow 0$ . To this end, every quantity evaluated at $\unicode[STIX]{x1D719}^{a,b}$ is then expanded in Taylor series, as

(4.35) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}f(\unicode[STIX]{x1D719}^{a})=f_{m}+\unicode[STIX]{x1D70E}f_{m}^{\prime }+\displaystyle \frac{\unicode[STIX]{x1D70E}^{2}}{2}f_{m}^{\prime \prime }+\displaystyle \frac{\unicode[STIX]{x1D70E}^{3}}{6}f_{m}^{\prime \prime \prime }+O(\unicode[STIX]{x1D70E}^{4}),\\ f(\unicode[STIX]{x1D719}^{b})=f_{m}-\unicode[STIX]{x1D70E}f_{m}^{\prime }+\displaystyle \frac{\unicode[STIX]{x1D70E}^{2}}{2}f_{m}^{\prime \prime }-\displaystyle \frac{\unicode[STIX]{x1D70E}^{3}}{6}f_{m}^{\prime \prime \prime }+O(\unicode[STIX]{x1D70E}^{4}),\end{array}\right\}\end{eqnarray}$$

where primes denote derivatives with respect to $\unicode[STIX]{x1D719}$ and the subscript $m$ stands for ‘evaluated at $\unicode[STIX]{x1D719}_{m0}$ ’, as in (4.6).

In particular, the functions $\unicode[STIX]{x1D6FC}^{\pm }$ and $\unicode[STIX]{x1D6FD}^{\pm }$ appearing in (3.27) become (see also appendix B)

(4.36a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}^{+}=\unicode[STIX]{x1D6FC}_{m}+\frac{\unicode[STIX]{x1D70E}^{2}}{2}\unicode[STIX]{x1D6FC}_{m}^{\prime \prime }+O(\unicode[STIX]{x1D70E}^{4}),\quad \unicode[STIX]{x1D6FC}^{-}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FC}_{m}^{\prime }+\frac{\unicode[STIX]{x1D70E}^{3}}{6}\unicode[STIX]{x1D6FC}_{m}^{\prime \prime \prime }+O(\unicode[STIX]{x1D70E}^{5}),\end{eqnarray}$$

(4.37a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}^{+}=\unicode[STIX]{x1D6FD}_{m}+\frac{\unicode[STIX]{x1D70E}^{2}}{2}\unicode[STIX]{x1D6FD}_{m}^{\prime \prime }+O(\unicode[STIX]{x1D70E}^{4}),\quad \unicode[STIX]{x1D6FD}^{-}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FD}_{m}^{\prime }+O(\unicode[STIX]{x1D70E}^{3}).\end{eqnarray}$$

With reference to (4.30), it can be easily shown that

(4.38a,b ) $$\begin{eqnarray}w_{b}=w_{b0}+\unicode[STIX]{x1D70E}^{2}w_{b2}+O(\unicode[STIX]{x1D70E}^{4}),\quad w_{s}=w_{s0}+\unicode[STIX]{x1D70E}^{2}w_{s2}+O(\unicode[STIX]{x1D70E}^{4}),\end{eqnarray}$$

(4.39) $$\begin{eqnarray}w_{sb}=\unicode[STIX]{x1D70E}^{2}w_{sb2}+O(\unicode[STIX]{x1D70E}^{4}).\end{eqnarray}$$

In the limit of vanishing $\unicode[STIX]{x1D70E}$ , we can then write

(4.40a,b ) $$\begin{eqnarray}w_{b}=w_{b0}=\unicode[STIX]{x1D6FE}\left[(\unicode[STIX]{x1D6FC}_{m}+\unicode[STIX]{x1D6FD}_{m})\frac{T_{B}}{T_{t0}}-\unicode[STIX]{x1D6FD}_{m}\frac{\text{i}k\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D703}_{c0}}\right],\quad w_{s}=w_{s0}=\unicode[STIX]{x1D6E4}\frac{2}{3}\frac{\unicode[STIX]{x1D6FC}_{m}}{\ln (2)}=\unicode[STIX]{x1D6FE}\frac{q_{u0}}{L_{a}}.\end{eqnarray}$$

Hence, as expected, for a well sorted two-grain mixture the bed eigenvalue tends to the one already discussed for the case of uniform sediment, the correction due to the effect of sorting appearing at $O(\unicode[STIX]{x1D70E}^{2})$ . Therefore, the growth rate is negative and the uniform flow is stable with respect to bed waves.

Less obviously, as soon as the mixture becomes even slightly heterogeneous a second, $O(\unicode[STIX]{x1D6E4})$ sorting eigenvalue appears, which is real and positive. Hence, at leading order the associated instability is marginal and perturbations propagate downstream without amplification nor decay. Moreover, the celerity increases with the sediment discharge, thus implying a faster downstream propagation of sorting perturbations as the Shields and the Froude numbers increase. Note that the estimate of the sorting wave celerity provided by (4.40) coincides with that of Ribberink (1987), who first studied the dynamics of small-amplitude perturbations of bed level and composition in river with non-uniform sediment.

The sorting instability appears at order $O(\unicode[STIX]{x1D70E}^{2})$ and is related to the imaginary part of the roughness-driven portion of the bed shear stress. Hence, at the leading order, the growth rate of the sorting eigenvalue is

(4.41) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{S}=-k\unicode[STIX]{x1D70E}^{2}\left[\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D6FC}_{m}^{\prime }+\unicode[STIX]{x1D6FD}_{m}^{\prime })\frac{T_{R}^{i}}{T_{t0}}-w_{sb2}^{i}\right],\end{eqnarray}$$

where $T_{R}$ is the part of the perturbation of the bed shear stress $T_{t1}$ proportional to the perturbation of the roughness height $R_{1}$ . Similarly, the celerity of the sorting wave reads

(4.42) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{S}=\unicode[STIX]{x1D6E4}\frac{2}{3}\frac{\unicode[STIX]{x1D6FC}_{m}}{\ln (2)}-\unicode[STIX]{x1D70E}^{2}\left[\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D6FC}_{m}^{\prime }+\unicode[STIX]{x1D6FD}_{m}^{\prime })\frac{T_{R}^{r}}{T_{t0}}-w_{sb2}^{r}\right],\end{eqnarray}$$

which provides the $O(\unicode[STIX]{x1D70E}^{2})$ correction of the estimate of the sorting wave celerity given by Ribberink (1987).

Moreover, neglecting as in (4.16) the normal stress $T_{n1}$ , we have

(4.43) $$\begin{eqnarray}\frac{T_{R}}{T_{t0}}=-R_{0}C^{2}\frac{2k^{2}C^{2}(Fr^{2}-1)+6\text{i}kFr^{2}}{k^{2}C^{4}(Fr^{2}-1)^{2}+9Fr^{4}}=R_{0}\frac{T_{B}}{T_{t0}}.\end{eqnarray}$$

Ignoring for the moment the $O(\unicode[STIX]{x1D6FE})$ contribution associated with $w_{sb2}$ in (4.41) and (4.42), we note that the behaviour is here reversed with respect to the bed eigenvalue: the base uniform flow is always unstable with respect to sorting perturbations, since the growth rate is proportional to the opposite of the imaginary part of the bed shear stress, which is negative (i.e. the stress perturbation leads the sorting wave).

The behaviour of the growth rate of the sorting eigenvalue in the large wavenumber range controls both the well posedness of the problem and the presence of a cutoff wavelength. Taking the limit of (4.43) as $k\rightarrow \infty$ we obtain

(4.44a,b ) $$\begin{eqnarray}\frac{T_{R}^{r}}{T_{t0}}\rightarrow \frac{2R_{0}}{(1-Fr^{2})},\quad \frac{kT_{R}^{i}}{T_{t0}}\rightarrow -\frac{6Fr^{2}R_{0}}{C^{2}(Fr^{2}-1)^{2}},\end{eqnarray}$$

so that the growth rate tends to a finite positive value and the solution is regular, except for the transcritical $Fr\rightarrow 1$ region, where, similarly to (4.17), resonance takes place

(4.45a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{S}\rightarrow \pm \infty \quad \text{as}~Fr\rightarrow 1^{\mp }\quad \unicode[STIX]{x1D6FA}_{S}\rightarrow -\infty .\end{eqnarray}$$

The latter disappears if the normal stress is included in the analysis, yielding

(4.46a,b ) $$\begin{eqnarray}\frac{T_{R}^{r}}{T_{t0}}\rightarrow \frac{2}{3N}\quad \frac{kT_{R}^{i}}{T_{t0}}\rightarrow -\frac{2C^{2}}{3}\frac{Fr^{2}-1+3R_{0}N}{N^{2}Fr^{2}}.\end{eqnarray}$$

Note that the instability is now limited to the range

(4.47a,b ) $$\begin{eqnarray}Fr>\sqrt{1-3R_{0}N}\simeq 1,\quad 3R_{0}N=d_{50}N/4\ll 1\end{eqnarray}$$

and that the solution remains regular even though a cutoff is still missing: as soon as the Froude number exceeds unity, the shortest modes are all equally unstable. In order to introduce a cutoff, the term $w_{sb2}$ , which includes the stabilizing effect of gravity, has to be accounted for in (4.41). However, by doing so, a second unstable region appears in the subcritical region.

As a matter of fact, sorting instability is related to an intricate balance between the two terms appearing in (4.41), which neither an order of magnitude analysis nor perturbation theory appear capable to disentangle. Nonetheless, the analysis of the weak-sorting case provides some useful information: (i) as the sediment mixture becomes more and more well sorted, the bed eigenvalue tends to the one obtained in the case of uniform sediment; (ii) sorting waves are found to be unstable and to propagate downstream with a faster wave speed with respect to bed waves; (iii) an $O(\unicode[STIX]{x1D70E}^{2})$ correction of the celerity of the sorting wave estimated by Ribberink (1987) is found; (iv) the growth rate is proportional to the variance of the size distribution; (v) sorting instability is regular, since the growth rate tends to a finite limit in the short-wave range.

5 Discussion of results

Figure 3. Stability plot of the sorting eigenvalue with varying $C$ and constant $\unicode[STIX]{x1D70E}=0.01$ ; $C=16$  (a) and 18 (b).

Stability plots represent the main output of the linear analysis. By plotting the isolines of the growth rate in the $(k,Fr)$ space, the regions of instability (i.e. positive growth rate) can be easily identified, bounded by the marginal curves (i.e. of vanishing growth rate). In figure 3 stability plots of the sorting eigenvalue are shown for a value of the standard deviation $\unicode[STIX]{x1D70E}$ equal to 0.01 and for two values of the non-dimensional Chézy coefficient $C$ . The growth rate is obtained from an iterative solution of the fully coupled eigenvalue problem (3.30). In each plot, the no transport area at the bottom corresponds to a base shear stress lower than the critical threshold for the coarser fraction, so that above the upper boundary of this region the sediment transport rate of both fractions is non-zero. Partial transport, whereby only the finer fraction is mobile, cannot be considered in the framework of a perturbative analysis. The white areas above the no transport line are associated with stable uniform flow (i.e. negative growth rate) with active sediment transport, whereas in the unstable regions the (positive) growth rate is displayed in shades of grey, the darkest the largest. White isolines are drawn for constant values of the growth rate (multiplied by $C^{2}/\unicode[STIX]{x1D70E}^{2}$ ) which increase on a logarithmic scale. Solid black lines identify marginal curves, whereas dashed lines mark the maximum amplification wavenumber.

Figure 4. Comparison between the stability plots of the ‘fast’ flow eigenvalue obtained with (a) and without (b) the inclusion of the depth-averaged normal Reynolds stress; $C=18$ , $\unicode[STIX]{x1D70E}=0.01$ .

Two distinct regions of instability can be observed in each plot, one in the subcritical and one in the supercritical regime, whereas in the transcritical region the uniform flow is stable with respect to sorting waves (and to bed waves as well). Both the unstable regions are bounded in the short-wave and long-wave ranges. As the conductance Chézy coefficient $C$ is raised, which corresponds to deeper flows and finer sediments, the two regions of instability move away from the critical ( $Fr=1$ ) line. Meanwhile, the subcritical region widens whereas the supercritical one slightly shrinks. As a whole, it seems relatively easier to observe unstable sorting waves in the supercritical regime, and with coarser sediments, consistently with laboratory observations. The wavenumber of maximum amplification, the one selected by the linear stability analysis, varies greatly with the Froude number, spanning approximately two orders of magnitude. Most frequent values are of order 10–100, showing that sorting waves are consistently shorter than bed waves, which are typically characterized by $O(1)$ wavenumbers.

Among the four roots of the complete eigensystem (3.30) there is only one more unstable eigenvalue, to be associated with the free-surface instability which leads, for values of the Froude number larger than 2, to the formation of roll waves. The corresponding stability plot is shown in figure 4(a), while figure 4(b) shows the effect of switching off the term associated with the depth-averaged normal Reynolds stress term. The crucial role played by this term in damping short-wave perturbations, thus introducing a cutoff and a wavelength selection mechanism (Needham & Merkin 1984; Balmforth & Vakil 2012), is more than evident. It must be remarked that each of the four eigenvalues, whether associated with a flow or with a sorting or bed instability, brings with itself, if unstable, a set of growing perturbations of velocity, depth, bed elevation and roughness all characterized by the same wavelength and wave speed. The associated eigenvectors define the relative amplitude and phase of each perturbations with respect to the others. Indeed, for Froude numbers larger than 2, both roll waves and sorting waves are found to be simultaneously unstable at a linear level.

Figure 5. Comparison between the stability plots obtained with the fully coupled (a) and with the quasi-steady (b) solutions for $C=18$ .

Figure 6. Stability plot of the quasi-steady sorting eigenvalue with varying $\unicode[STIX]{x1D70E}$ value and constant $C=18$ ; (a)  $\unicode[STIX]{x1D70E}=1$ ; (b)  $\unicode[STIX]{x1D70E}=3$ ; (c)  $\unicode[STIX]{x1D70E}=5$ ; (d) the weak-sorting solution for the same value of $C$ .

In figure 5 a comparison between the growth rate of the sorting eigenvalue obtained by means of the fully coupled problem (3.30) and of the decoupled quasi-steady problem (4.24) is presented. The quasi-steady solution presents substantially the same instability pattern observed in the coupled case; however, while the subcritical region remains mostly unchanged, the supercritical region is slightly enlarged. Note that the free-surface roll-wave instability associated with the fastest flow eigenvalue, is completely wiped out by the quasi-steady approximation, whereby the flow is assumed to adjust instantaneously to perturbations of bed elevation and roughness without being allowed to develop instabilities of its own.

Figure 7. Amplitude (a,c) and phase (b,d) of the roughness height perturbation $R_{0}R_{s}$ for a unitary amplitude perturbation of the bed elevation; (a,b)  $C=16$ , (c,d)  $C=18$ .

The effect of increasing $\unicode[STIX]{x1D70E}$ and thus the relative distance in $\unicode[STIX]{x1D719}$ -units between the two fractions is shown in figure 6. As in all the other stability plots presented above, the growth rate is multiplied by $C^{2}/\unicode[STIX]{x1D70E}^{2}$ , so that the small variations observed are related to effects of order $\unicode[STIX]{x1D70E}^{4}$ or greater. Indeed, the effect of the mixture heterogeneity is rather small and becomes visible only for values of the standard deviation of $O(1)$ . The no-transport threshold slightly increases with $\unicode[STIX]{x1D70E}$ , because a higher Shields stress (and consequently Froude number) is required to set the coarsest fraction in motion. The plot in figure 6(d) presents the solution obtained by means of the weak-sorting approximation, which can be considered as the limit of the quasi-steady solution for vanishing $\unicode[STIX]{x1D70E}$ . Note that this stability plot is almost indistinguishable from that of figure 6(a) (quasi-steady solution, $\unicode[STIX]{x1D70E}=1$ ), showing that the weak-sorting solution holds for relatively large values of the standard deviation. We recall that a value of $\unicode[STIX]{x1D70E}$ equal to 5 corresponds to a finer (coarser) fraction thirty-two times smaller (larger) than the median diameter of the mixture $d_{50}$ . If the analysis is pushed toward an even larger heterogeneity, the solution eventually breaks down, because the averaged dynamics of the mixture is not well represented by $d_{5}0$ when the two fractions are too far apart in the $\unicode[STIX]{x1D719}$ -scale. Nonetheless, according to the verbal classification scale of sorting introduced by Folk & Ward (1957), a standard deviation larger than 4 corresponds to an ‘extremely poorly sorted’ mixture, and a standard deviation larger than 2 to a ‘poorly sorted’ mixture. Hence, we can safely state that the model holds for very well sorted up to very poorly sorted mixtures and that the weak-sorting approximation can be used up to poorly sorted sediment.

Figure 7 shows the behaviour of the eigenvector associated with the sorting eigenvalue. The amplitude and phase of the roughness height perturbation $R_{0}R_{s}$ corresponding to a unitary amplitude perturbation of the bed elevation are plotted inside the instability regions for two different values of $C$ . The amplitude of the roughness height perturbation (a,c) is typically smaller in the subcritical than in the supercritical region and globally increases with $C$ . Moreover, in the supercritical range the amplitude tends to decrease with the wavenumber, whereas it is almost unaffected in the subcritical regime. The plots in figure 7(b,d) show the relative phase of the roughness with respect to the perturbation of the bed amplitude. Irrespective of the Chèzy coefficient, in the supercritical region the phase is in the range $0{-}\unicode[STIX]{x03C0}$ , whereas in the subcritical region it is in the range $\unicode[STIX]{x03C0}{-}3\unicode[STIX]{x03C0}/2$ , with values decreasing (increasing) for larger wavenumber in the supercritical (subcritical) regime. We recall that a phase lag of $\unicode[STIX]{x03C0}$ corresponds to out-of-phase perturbations, which means that the roughness height is maximum where bed elevation is minimum and vice versa. Moreover, the maximum of the roughness height corresponds to the minimum of the mean grain size (in the $\unicode[STIX]{x1D719}$ -scale) $\unicode[STIX]{x1D719}_{m}$ and to the maximum areal concentration of the coarser fraction $F^{b}$ . Hence, the present results are in good agreement with the experimental measurements of Bennett & Bridge (1995), who described bedload sheets in the subcritical regime as having ‘relatively fine-grained crests and coarse-grained troughs’ (see also Kuhnle & Southard (1988) and Kuhnle et al. (2006)).

Table 2. Main hydraulics conditions and characteristics of observed bedload sheets. In the last column, an estimate of the celerity using (4.40). BB, Bennett & Bridge (1995); C, Carbonari (2019); II, Iseya & Ikeda (1987); Ka, Kuhnle et al. (2006); KS, Kuhnle & Southard (1988); Wa, Whiting et al. (1988) (DC: Duck Creek, MC: Muddy Creek); W, Wilcock (1992).

Table 2 summarizes the hydraulic and geometrical characteristics of bedload sheets observed in other studies, made non-dimensional using the fluid mass density and the measured values of the flow depth and area velocity. In the last column, a prediction of the celerity of the sorting waves according to (4.40) is provided.

It must be noted that a direct comparison with the experimental results is made extremely difficult by the peculiarity of bedload sheets, which are ephemeral sorting patterns characterized by small heights and short wavelengths, easily disrupted, or at least heavily altered, by the interactions with bed forms and with the sorting patterns associated with them. This aspect is particularly evident in the detailed description of Run 6 contained in the seminal paper of Iseya & Ikeda (1987), which deserves some further discussion in light of the present results.

Indeed, we speculate that the rhythmic alternations of ‘congested’ and ‘smooth’ states described by Iseya & Ikeda (1987) can be interpreted as a topographic pattern, as shown by the longitudinal profile of the bed displayed in figure 6 of the paper, where a periodic change of slope (and thus of bed elevation) is quite evident. Since the bed is heterogeneous, there is necessarily a sorting pattern associated with this bed form, whereby gravel (sand) accumulates on the crests (troughs) with a spacing of approximately 2 m. Superimposed on it, there is a second sorting pattern, associated with the so-called ‘transitional’ state, which is described as ‘gravel jams one or two grain diameters thick’, migrating over a smooth sand bed with a spacing of about 8 cm. The sorting patterns observed by Iseya & Ikeda (1987) are, in our view, to be considered the result of the two different instabilities discussed in the present contribution: the bed eigenvalue has to be associated with the topographic long wave, whereas the sorting eigenvalue is representative of the sorting short wave of the bedload-sheet kind.

As a conclusive comment, in light of the experimental observations presented in table 2, we can state that the present theory correctly predicts that: (i) bedload sheets occur in both the sub-critical and super-critical regimes; (ii) bedload sheet wavenumbers are of $O(1)$ or larger; (iii) the celerity of the sorting wave increases with the Shields number and with the bedload transport rate (4.40). Moreover, there is a good, although qualitative, agreement between the theoretical predictions and the measured values of the wave speed.

On the downside, the theory seems to underestimate the wavelength of maximum amplification and it remains limited to relatively small grain-to-depth ratios and to moderate Shields numbers. Moreover, the sediment transport model is based on a simplified version of the Hirano (1971) model, whereby interactions between the active layer and the substrate are neglected and thus the analysis is unsuitable to study cases with substantial bed aggradation or degradation. Finally, the shallow-water flow model adopted does not allow for a full description of the phenomenon because the bed eigenvalue is consistently stable and no bed waves form. On the other hand, the algebraic eigenvalue problem obtained making use of the shallow-water approximation provides a deep insight into the stability of the sorting wave, which can be of great use for a future analysis based upon a more complete rotational flow model.

6 Conclusions

A linear stability analysis of a uniform flow over a bed composed of an even mixture of two grain sizes is presented. The complete eigensystem consists in a fourth degree characteristic polynomial, the roots of which can be associated with flow and bed instabilities. In particular, two eigenvalues present regions of positive growth rate, one associated with the amplification of a free-surface perturbation of the roll-wave kind, the other with the formation of sorting waves. Instability is shown to be well posed from a mathematical point of view and a cutoff in the short wavelength range is present, provided that the effect of gravity on the sediment and the diffusive term in the flow equation, which formally emerges by the integration along the depth of the normal Reynolds stress, are included in the analysis. Moreover, for the sorting eigenvalue, it has been shown that the analysis is likely to fail for relatively coarse sediment (or for relatively shallow flows) and for small values of the base Shields stress.

A simplified quasi-steady analysis is developed, whereby flow derivatives are dropped in the flow equation and the flow is assumed to adapt instantaneously to changes in bed elevation and roughness. In this case, the stability analysis reduces to a double eigenvalue problem, one of which can be associated with the formation of bedload sheets, rhythmic alternations of coarse and fine stripes with negligible variations of the bed elevation, the other with sorting over bed forms. The distinction of the two can be built on the basis of the related eigenvectors. In the present shallow-water framework, only the former is found to be unstable, predicting the formation of bedload sheets in two separate regions of instability, one for subcritical and the other for supercritical base flows. In both regimes, the perturbations travel downstream at a faster rate with respect to bed waves but at a slower speed with respect to the depth-average flow velocity. On the contrary, the bed eigenvalue is found to be invariably stable, due to the inability of the shallow-water flow model to correctly predict the phase lag of the bed shear stress perturbation driven by changes in bed elevation.

The weak-sorting case is also investigated, whereby the two fractions are assumed to collapse on the median grain size as the standard deviation of the mixture vanishes. As the sediment becomes more and more well sorted, the bed eigenvalue naturally devolves into the morphodynamic eigenvalue of the homogeneous case, the instability of which is controlled by a balance between the component of the bed shear stress associated with the perturbation of the bed elevation and the stabilizing effect of gravity. Instability of the sorting wave appears at $O(\unicode[STIX]{x1D70E}^{2})$ and is mainly related to the component of the bed shear stress associated with a perturbation of the bed roughness. The weak-sorting solution has also been used to show that the analysis is limited to relatively small values of the median grain-to-depth ratio and to relatively large values of the Shields stress. When this inequality is satisfied, the weak-sorting solution behaves surprisingly well, providing an excellent representation of the solution from very well sorted up to poorly sorted mixtures.

Finally, a comparison of the present results with some laboratory and field observations of bedload sheets has been attempted in terms of wavelength and celerity of the sorting waves, showing an acceptable, though merely qualitative, agreement.