3.1. Low settling velocity $(v_s/\kappa \ll 1)$

Figure 3.

Rescaled perturbations to the concentration field in the regime $v_s/\kappa \ll 1$ . Panels (a–c) report data for $ \textit{EQ}_1$ . (a,b) Vertical profiles of $\hat {c} = c - c_0$ at $x = 0$ , $z =$ (a) $0$ and (b) $L_z/2$ , rescaled by $v_s/\kappa$ . (c) Vertical profiles of the absolute maximum of $\hat {c}$ over the horizontal coordinates $x$ and $z$ , rescaled by $v_s/\kappa$ . Panels (d–f) report the corresponding data for the upper branch state, $ \textit{EQ}_2$ .

When $v_s = 0$ , $c_0 \equiv 1$ and (3.1) admits the trivial homogeneous perturbation $\hat {c} \equiv 0$ , regardless of the flow state. Likewise, as shown in figure 2(a), when $v_s/\kappa$ is small but non-zero, the concentration field is a small correction to this homogeneous state. We expand the concentration field as $c = c_0(y) + \epsilon c_1(x,y,z) + \ldots$ , where $\epsilon$ is a positive order parameter assumed to be small, relative to the fields $c_0$ and $c_1$ . After substituting into (2.1c ), neglecting $O(\epsilon ^2)$ terms and using (2.5c ), we obtain

(3.2) \begin{align} \epsilon (y+\hat {u} - a)\frac {\partial c_1}{\partial x} + \hat {v} \left ( -\frac {v_s}{\kappa }c_0 + \epsilon \frac {\partial c_1}{\partial y} \right ) -\epsilon v_s \frac {\partial c_1}{\partial y} + \epsilon \hat {w}\frac {\partial c_1}{\partial z} = \epsilon \kappa {\nabla} ^2 c_1. \end{align}

Provided $\textit{Sc} \sim 1$ , the settling term $\epsilon v_s\partial c_1/\partial y$ can be considered negligible, on the grounds that $v_s \ll \kappa \ll 1$ , so $\epsilon v_s \ll \epsilon$ . The size of the remaining $O(\epsilon )$ terms in (3.2) must match the vertical advection of the base concentration profile, so $\epsilon = v_s/\kappa$ . This leads to

(3.3) \begin{align} (y + \hat {u} - a) \frac {\partial c_1}{\partial x} + \hat {v}\frac {\partial c_1}{\partial y} + \hat {w}\frac {\partial c_1}{\partial z} = \hat {v} + \kappa {\nabla} ^2 c_1, \end{align}

where we have used the fact that $c_0 = 1 + O(\epsilon )$ when $v_s/\kappa \ll 1$ . We demonstrate the inferred scaling regime in figure 3 for both the solutions $ \textit{EQ}_1$ and $ \textit{EQ}_2$ . For example, in figure 3(a,d), plotted profiles of $\hat {c}$ at $(x,z) = (0,0)$ collapse under rescaling by $v_s/\kappa$ for $v_s/\kappa = 10^{-1}$ , $10^{-2}$ and $10^{-3}$ . For reference, we also show the corresponding profiles for $v_s/\kappa =1$ (dotted grey), which do not collapse onto the limiting curves, but nevertheless lie close by. Profiles taken at other points in the domain exhibit similar collapse – figure 3(b,e) demonstrates the case $(x,z)=(0,L_z/2)$ . To summarise this, in figure 3(c,f) we plot $(v_s/\kappa )^{-1}\max _{x,z}|\hat {c}|$ for the same $v_s/\kappa$ values, which attain a limiting profile for $v_s/\kappa \lesssim 10^{-1}$ .

Further insight into the asymptotic structure may be obtained by appealing to the SSP/VWI scalings for the velocity fields. Noting that these are uncoupled to $c$ , we substitute (2.6a –c ) into (3.3). Away from the asymptotically shrinking critical layer, it may be deduced that $\partial c_1 / \partial x \to 0$ as $ \textit{Re} \to \infty$ . The resulting equation may be averaged to obtain the balance that sets $\bar {c}_1$ , i.e.

(3.4) \begin{align} V\frac {\partial \bar {c}_1}{\partial y} + W\frac {\partial \bar {c}_1}{\partial z} = V + \frac {1}{\textit{Sc}} \left ( \frac {\partial ^2 \bar {c}_1}{\partial y^2} + \frac {\partial ^2 \bar {c}_1}{\partial z^2} \right )\!, \end{align}

where $ \textit{Re} ^{-1}(V, W) = (\bar {v}, \bar {w})$ and subject to the boundary conditions $\partial \bar {c}_1/\partial y=0$ at $y = \pm 1$ . We solve this numerically, using a pseudospectral representation of $\bar {c}_1(y,z)$ with Fourier modes in $z$ and Chebyshev polynomials in $y$ . The velocity fields $\bar {v}$ and $\bar {w}$ are provided by the known solutions for either $ \textit{EQ}_1$ or $ \textit{EQ}_2$ . This leads to asymptotic predictions for the concentration perturbations of each state at high $ \textit{Re}$ and low $v_s/\kappa$ , which are compared with the corresponding $c - c_0$ fields for the full states at $v_s/\kappa = 10^{-2}$ in figure 4. The $ \textit{Re} = 400$ cases are given in panels (a,d) and exhibit reasonable agreement. Sediment concentration is elevated with respect to the base profile towards the upper wall and diminished towards the lower wall. By referring to figure 1, we note that these trends are accentuated according to the locations of the rolls, which redistribute particles to form dilute and concentrated patches around the points $(y,z) = (-1,0) \equiv (-1,L_z)$ and $(1,L_z/2)$ , respectively. As $ \textit{Re}$ is increased to $1200$ in figure 4(b,e) and then to $2000$ in figure 4(c,f), the asymptotic solutions and full ECS states converge to one another.

Figure 4.

Comparison between the ECS concentration perturbations and asymptotic solutions to (3.4) in the regime $v_s/\kappa \ll 1$ for increasing $ \textit{Re}$ . (a–c) Filled contours from blue (low) to red (high) show $\bar {c} - c_0$ for $ \textit{EQ}_1$ at $v_s/\kappa = 10^{-2}$ , equally spaced between $\pm 0.5 v_s/\kappa$ and centred at zero (white). Isolines of $(v_s/\kappa )\bar {c}_1$ are overlaid at the same contour values with solid and dashed lines representing positive and negative values, respectively. The Reynolds numbers are $ \textit{Re} =$ (a) $400$ , (b) $1200$ and (c) $2000$ . Panels (d–f) present the corresponding data for $ \textit{EQ}_2$ in the same style, except the contours are taken at $\pm 0.9 v_s/\kappa$ .

3.2. High settling velocity $(v_s/\kappa \gg 1)$

When the settling velocity is relatively large, a layer of sediment forms at the bottom wall (see figure 2 c). Increasing $v_s$ concentrates the boundary layer, whose characteristic width $\varepsilon$ (to be determined below) appears to shrink to zero in the limit $v_s\to \infty$ . Therefore, we analyse this region by rescaling the vertical coordinate at the wall, defining $Y = (y + 1) / \varepsilon$ . The concentration equation then becomes

(3.5) \begin{align} \frac {\partial ^2 \hat {c}}{\partial Y^2} + \varepsilon \frac {v_s}{\kappa }\frac {\partial \hat {c}}{\partial Y}& = \frac {\varepsilon ^2}{\kappa }(-1 - a + \hat {u} + \varepsilon Y) \frac {\partial \hat {c}}{\partial x} + \frac {\varepsilon }{\kappa } \hat {v}\left ( \frac {\partial c_0}{\partial Y}+ \frac {\partial \hat {c}}{\partial Y} \right ) +\frac {\varepsilon ^2}{\kappa }\hat {w}\frac {\partial \hat {c}}{\partial z} \nonumber \\ &\quad -\varepsilon ^2\left ( \frac {\partial ^2 \hat {c}}{\partial x^2}+ \frac {\partial ^2 \hat {c}}{\partial z^2} \right )\!. \end{align}

In order for non-trivial solutions to exist as $v_s/\kappa$ becomes high, $\varepsilon = \kappa / v_s$ . (Although $\partial c_0/\partial Y$ also depends on $v_s/\kappa$ , it will become clear shortly that this term is subdominant.) The only terms that survive to leading order as $\varepsilon \to 0$ are vertical settling and vertical diffusion:

(3.6) \begin{align} \frac {\partial ^2 \hat {c}}{\partial Y^2} + \frac {\partial \hat {c}}{\partial Y} = 0. \end{align}

This equation has the general solution $\hat {c}(x,Y,z) = A(x,z) \mathrm{e}^{-Y} + B(x,z)$ . The no-flux boundary at $Y=0$ implies that $B \equiv 0$ . The amplitude $A(x,z)$ is dictated by the right-hand side terms of (3.5), which are present at higher order. The magnitudes of these terms are dependent on the velocity field. For sufficiently small $\varepsilon$ , the VWI critical layer (whose location is fixed) must lie outside the boundary layer. Consequently, the velocity structure is streamwise independent within this region and is given by the streak and roll components of (2.6a –c). These may be Taylor expanded at the wall to infer that

(3.7a) \begin{align} \hat {u}(Y,z) = \varepsilon u_1(z) Y + \ldots , \\[-30pt] \nonumber \end{align}

(3.7b) \begin{align} \hat {v}(Y,z) = \textit{Re} ^{-1}\varepsilon ^2 v_2(z) Y^2 + \ldots , \\[-30pt] \nonumber \end{align}

(3.7c) \begin{align} \hat {w}(Y,z) = \textit{Re} ^{-1}\varepsilon w_1(z) Y+\ldots , \\[0pt] \nonumber \end{align}

where $u_1$ , $v_2$ and $w_1$ are unknown functions that are $O(1)$ with respect to both $\varepsilon$ and $ \textit{Re}$ . The quadratic dependence of the wall-normal velocity field follows from the incompressibility of the flow. Substitution of these expressions into (3.5), along with the asymptotic formula for the base concentration field $c_0(Y)= 2\varepsilon ^{-1}\mathrm{e}^{-Y} + \ldots$ , implies that the right-hand side terms in (3.1) enter at $O(\varepsilon ^2)$ . Therefore, the concentration perturbation has the asymptotic form

(3.8) \begin{align} \hat {c}(x,Y,z) = A(x,z) \mathrm{e}^{-Y} + \varepsilon ^2 c_2(x,Y,z) + \ldots , \end{align}

where $c_2$ is an arbitrary $O(1)$ field. Although streamwise dependence has been retained in this expansion, in practice the concentration is streamwise invariant to leading order, because its structure is controlled by the $x$ -invariant velocity fields. Consequently, to effectively characterise $\hat {c}$ at leading order, it is enough to identify $\bar {A}(z)$ , which may be obtained by substituting (3.8) into (3.5) (along with the expressions for the velocities and $c_0$ ) and streamwise averaging. At $O(\varepsilon ^2)$ , this leads to

(3.9) \begin{align} \frac {\partial ^2 \bar c_2}{\partial Y^2}+ \frac {\partial \bar c_2}{\partial Y} =-2\textit{Sc}\, v_2(z) Y^2 \mathrm{e}^{-Y} - \bar A''(z)\mathrm{e}^{-Y}. \end{align}

Integrating this equation in the wall-normal direction between $0$ and $Y$ gives

(3.10) \begin{align} \frac {\partial \bar c_2}{\partial Y^2}+ \bar c_2 =2\textit{Sc}\, v_2(z) \mathrm{e}^{-Y}(Y^2 + 2Y + 2) + \bar {A}''(z)(\mathrm{e}^{-Y} - 1) -4 \textit{Sc} v_2(z). \end{align}

Outside the boundary layer, the concentration decays to zero. Therefore, in order for the concentration fields to match as $Y\to \infty$ , it must be the case that $4\textit{Sc} v_2 + \bar {A}'' \equiv 0$ . Note that, due to incompressibility, $v_2 = -({1}/{2}) w_1'$ . Consequently,

(3.11) \begin{align} \bar {A}(z) = \mathcal{I}(z) + \alpha z - \beta , \quad \mathrm{where}\quad \mathcal{I}(z) = 2 \textit{Sc} \int _0^z w_1(\tilde z) \,\mathrm{d}\tilde z, \end{align}

with $\alpha$ and $\beta$ constants of integration. The first of these is necessarily zero, since the boundary conditions require $\bar {A}(z)$ to be an $L_z$ -periodic function. The second constant, $\beta$ , is set by conservation of concentration, which mandates that the integral of $\hat {c}$ over the whole domain vanishes. Therefore, $\beta = \langle \mathcal{I} \rangle$ and we are left with the following asymptotic solution for the streamwise-averaged concentration field, which we write in full, up to $O(1)$ :

(3.12) \begin{align} \bar {c}(Y,z) = 2\mathrm{e}^{-Y}\left [ \frac {v_s}{\kappa } + \textit{Sc}\left ( \int ^z_0 w_1(\tilde z) \, \mathrm{d}\tilde z -\frac {1}{L_z}\int _0^{L_z}\!\! \int _0^z w_1(\tilde z)\,\mathrm{d}\tilde z \,\mathrm{d}z \right ) \right ]\!. \end{align}

Since the sign of $\bar {A}'(z)$ is everywhere identical to the sign of $w_1(z)$ , the maxima and minima of $\bar {A}(z)$ lie at cross-stream stagnation points $z = z^*$ . In particular, maxima must occur when $w_1'(z^*) \lt 0$ , meaning that sediment preferentially accumulates at interfaces between rolls where spanwise velocity is oriented towards $z^*$ from both sides, just as we have seen from the illustrative plots of $ \textit{EQ}_1$ in figure 2.

The sediment distribution within the boundary layer can be visualised more clearly by plotting contours of the streamwise-averaged perturbation to the base concentration. In figure 5(a,b) we plot these data for $ \textit{EQ}_1$ at $v_s/\kappa =14$ and $40$ , respectively, obtaining excellent agreement with the corresponding asymptotic solution of (3.12).

Figure 5.

Concentration fields of $ \textit{EQ}_1$ in the boundary layer at the lower wall for high $v_s/\kappa =$ (a) $14$ and (b) $40$ . Filled contours show the perturbation after taking a streamwise average, $\bar {c}(Y,z) - c_0(Y)$ , with values spaced at intervals of $2$ and centred around zero (white). Black contour lines at the same levels show the $O(1)$ component of the asymptotic solution from (3.12). Dashed and solid lines indicate negative and positive values, respectively.

3.3. Sediment transport

So far we have seen that the velocity structures of ECS maintain a distribution of sediment that differs from the laminar base profile in ways that are analytically tractable in the regimes of low and high settling velocity. It is natural to ask how this affects the transport of sediment in the channel. This question depends on the reference frame of an observer monitoring the flow. In applications, it is typical for the bottom boundary to be stationary from an observer’s perspective, rather than the centreline $y=0$ , as is the case for our configuration. The total streamwise flux $Q_x$ , measured in a frame with a stationary bottom boundary, is

(3.13) \begin{align} Q_x = \frac {1}{2}\int _{-1}^1 \langle \overline {(u + 1) c} \rangle \, \mathrm{d}y, \end{align}

which gives a bulk measure for the rate at which sediment is transported by the flow. The portion, $\hat {Q}_x$ , of this flux that is attributable to the presence of the ECS is

(3.14) \begin{align} \hat {Q}_x = Q_x - \frac {1}{2}\int ^{1}_{-1}(y+1)c_0\,\mathrm{d}y. \end{align}

Another important value is the vertical flux $Q_y$ , given by

(3.15) \begin{align} Q_y = \frac {1}{2}\int _{-1}^1 \langle \overline {vc} \rangle \,\mathrm{d}y, \end{align}

since this quantifies the amount of sediment uplifted by the ECS against the flux from particle settling and thereby dictates the distribution of sediment available for streamwise transport.

Figure 6.

Sediment fluxes for $ \textit{EQ}_1$ at different $v_s/\kappa$ values (coloured points). The grey dashed lines are $Q_x = 1 - 0.176 v_s/\kappa$ and $Q_x = 1.429 \kappa / v_s$ , computed (and given to 3 significant figures) from the asymptotic approximations in (3.17) and (3.19) for streamwise fluxes in the limiting settling regimes. Their counterparts in dotted grey are the same lines shifted down by $1-v_s/(3\kappa )$ (see (3.18)) and $\kappa / v_s$ , respectively, thereby giving the flux $\hat {Q}_x$ attributable to the ECS perturbation field. The grey dash–dotted lines are $Q_y = 8.72\times 10^{-4} v_s/\kappa$ and $Q_y = 0.0841(\kappa /v_s)^3$ – the asymptotic predictions for the vertical fluxes, given in (3.16) and (3.20), respectively.

In figure 6 we show the fluxes for $ \textit{EQ}_1$ at different values of settling velocity. Total streamwise sediment flux (see (3.13)) is plotted with blue crosses. This is a decreasing function of $v_s/\kappa$ , because higher settling velocities decrease the amount of sediment in the upper portions of the channel, leading to lower rates of transport. Maximum flux occurs in the limit of vanishing settling velocity, where the particle concentration becomes homogeneous. With orange circles, we plot $\hat {Q}_x$ – the corresponding values of the streamwise sediment flux that may be attributed to the ECS (see (3.14)). This component decays to zero in both extremal regimes of settling velocity, attaining its maximum at $v_s/\kappa \approx 2$ , as does the vertical sediment flux (see (3.15)), plotted with green squares. The non-monotonicity of the fluxes $\hat {Q}_x$ and $Q_y$ is consistent with the following physical interpretation. At low $v_s/\kappa$ , the ECS can do little to redistribute the homogeneous base concentration profile and, consequently, has a negligible effect on the sediment transport. Increasing $v_s/\kappa$ towards intermediate values enhances the impact of the ECS, whose vertical fluxes maintain sediment concentrations in the channel that would otherwise settle towards lower values. However, at high settling velocities, accumulation of sediment at the lower wall dominates and the trend reverses. These insights generalise to other solutions and can be made precise by using our asymptotic analyses of the passive concentration fields in §§ 3.1 and 3.2 to derive theoretical estimates of the limiting behaviour in each regime. These are the grey lines in figure 6, which are in excellent agreement with the data in both regimes. Their formulae are given below.

In the low settling velocity regime, the concentration field has the form $c = 1 + (v_s/\kappa ) \bar {c}_1 + \ldots$ , with $\bar {c}_1$ set by (3.4). Therefore, to leading order in $v_s/\kappa$ , the vertical flux is

(3.16) \begin{align} Q_y = \frac {v_s}{\kappa } \boldsymbol{\cdot }\frac {1}{2} \!\int _{-1}^1 \langle \bar {v} \bar {c}_1 \rangle \,\mathrm{d}y. \end{align}

The integrand in this expression is independent of $v_s/\kappa$ . Its value may be obtained for any ECS by numerically solving (3.4) for $\bar {c}_1$ and using an appropriate quadrature to evaluate the integrals. This was computed for $ \textit{EQ}_1$ in figure 6, using the solution obtained in the previous subsection (plotted in figure 4 a). Similarly, the streamwise flux may be written as

(3.17) \begin{align} Q_x = 1 + \frac {v_s}{\kappa } \boldsymbol{\cdot }\frac {1}{2} \int _{-1}^{1} \langle -y\bar {u}+ \bar {u}\bar {c}_1 \rangle \,\mathrm{d}y. \end{align}

The first term in this expression combined with the $\langle -y\bar {u}\rangle$ part of the integral comprises the transport of the base concentration field, which is positive. The final term, corresponding to the transport of the concentration perturbation, is expected to be negative, since the vortex structure of the ECS, which advects low-speed streaks together with high concentrations (and vice versa; see figures 1 and 2), causes $\bar {u}$ and $\bar {c}_1$ to be inversely correlated. However, the total transport due to the ECS, which is

(3.18) \begin{align} \hat {Q}_x = Q_x - 1 + \frac {1}{3}\frac {v_s}{\kappa } = \frac {v_s}{\kappa }\boldsymbol{\cdot }\frac {1}{2}\int _{-1}^1 -y\langle \bar {\hat {u}} + \bar {u}\bar {c}_1 \rangle \,\mathrm{d}y, \end{align}

to leading order, is positive in figure 6 because the advection of the base concentration by the streamwise velocity perturbation (first part of the integrand in (3.18)) is positive and outweighs the contribution from the transport of the concentration perturbation.

In the limit of high settling velocity, the sediment localises towards a boundary layer of thickness $\varepsilon = \kappa / v_s \ll 1$ at the lower wall. We determine the sediment fluxes to leading order in $\varepsilon$ by substituting the expanded velocity fields of (3.7a –c ), together with the asymptotic solution of (3.12), into the transport integrals and simplifying. Firstly, the streamwise flux (see (3.13)) becomes

(3.19) \begin{align} Q_x = \frac {\kappa }{v_s} \left ( 1 + \langle u_1 \rangle \right ). \end{align}

The first term of this expression is the flux for purely laminar flow and the second gives the particular contribution of the ECS (i.e. $\hat {Q}_x = Q_x - \kappa / v_s$ ). For $ \textit{EQ}_1$ , this latter term comprises 30 % of the total flux, while for $ \textit{EQ}_2$ , the proportion is roughly 60 %. Using the same methods for the vertical flux, we obtain

(3.20) \begin{align} Q_y = \frac {2\textit{Sc}}{\textit{Re} }\left (\frac {\kappa }{v_s}\right )^{\!3} \frac {1}{L_z} \int _0^{L_z}v_2\left (\mathcal{I} - \langle \mathcal{I} \rangle \right ) \mathrm{d}z, \end{align}

where the $z$ -dependent function $\mathcal{I}$ was defined in (3.11). The rather acute inverse cubic dependence on settling velocity in this case arises from the fact that within the $O(\kappa /v_s)$ boundary layer over which the sediment is distributed the vertical velocity is $O(\kappa ^2/v_s^2)$ .

4.1. Solutions with $v_s/\kappa = 1$

Our numerical continuations begin away from either of the limiting regimes of low and high settling velocity, by studying solutions with $v_s/\kappa = 1$ , which we take as an initial illustrative parameter choice. In figure 7 we plot some solution families connected to $ \textit{EQ}_1$ and $ \textit{EQ}_2$ at $v_s/\kappa = 1$ in terms of $ \textit{Ri}_b$ versus the mean wall stress of the perturbation to laminar flow $\hat {\tau }_y$ . Both $ \textit{EQ}_1$ (solid blue) and $ \textit{EQ}_2$ (solid green) emerge from saddle-node bifurcations in $ \textit{Ri}_b$ , at different positive $ \textit{Ri}_b$ values (upward-pointing triangles). However, in the particular projection shown in the plot, they are not directly connected. This contrasts with continuation in $ \textit{Re}$ in homogeneous plane Couette flow ( $ \textit{Ri}_b = 0$ ), where they are the lower and upper branches of a saddle node in $ \textit{Re}$ (at $ \textit{Re} \approx 218.5$ ; see Gibson et al. Reference Gibson, Halcrow and Cvitanović2009). As $ \textit{Ri}_b$ increases, the stress of both states decreases from the values of $ \textit{EQ}_1$ and $ \textit{EQ}_2$ , respectively. They reach their respective saddle nodes at $ \textit{Ri}_b = 8.01\times 10^{-4}$ and $0.0127$ , before turning back towards lower $ \textit{Ri}_b$ . It is noteworthy that these maximum $ \textit{Ri}_b$ values are relatively small, indicating that relatively weak buoyancy effects are sufficient to disrupt unstable equilibria in this case. Both solution families have a positive wave speed $a\gt 0$ when they pass back through $ \textit{Ri}_b=0$ , giving rise to two unstratified travelling waves, which we call $\textit{TW}_1^+$ and $\textit{TW}_2^+$ . These two states were originally discovered by Viswanath (Reference Viswanath2008) and Gibson et al. (Reference Gibson, Halcrow and Cvitanović2009), respectively, though not through direct connections to $ \textit{EQ}_1/\textit{EQ}_2$ . Our labelling convention for these waves diverges from Gibson et al. (Reference Gibson, Halcrow and Cvitanović2009) (who called them $\textit{TW}_2$ and $\textit{TW}_3$ ) in order to emphasise their connection with the equilibria and distinguish them from symmetry-related solutions with negative wave velocities ( $\textit{TW}_1^-$ and $\textit{TW}_2^-$ , discussed below), since this becomes necessary when the states are continued in $ \textit{Ri}_b$ .

Figure 7.

Parametric continuation of $ \textit{EQ}_1$ (blue) and $ \textit{EQ}_2$ (green) in $ \textit{Ri}_b$ as a function of mean excess wall-normal stress $\hat \tau _y$ for $v_s/\kappa = 1$ . The red dashed curves are two-dimensional roll solutions that bifurcate off from linear instabilities of the base flow at $ \textit{Ri}_b = -9.08\times 10^{-4}$ and $-5.53\times 10^{-3}$ . Inlaid figures display data for states at the corresponding points indicated on the solid curves. Filled contours show $\partial \bar c / \partial y$ using $10$ equispaced intervals from $-2.2$ (dark) to $0$ (light). These are overlaid with positive (solid) and negative (dashed) isolines of the mean streamwise vorticity $\partial \bar {w}/\partial y - \partial \bar {v}/\partial z$ , at the values $\pm 0.08 n$ for $n = 1,\ldots ,5$ . The velocity fields of the states with $ \textit{Ri}_b \geqslant 0$ are also plotted later in figure 12.

Inlaid plots in figure 7 show the mean streamwise vorticity (using solid and dashed contours for positive and negative values, respectively) and mean wall-normal concentration gradient (filled contours) of the states at various points along the solution branches. In each case, concentration gradients are strongest at the lower wall, where sediment preferentially accumulates, due to settling. The two branches of each saddle node indicate two ways in which states adapt to the presence of stratification. The trends are clearest along the $ \textit{EQ}_2/\textit{TW}_2^+$ branches, but are evident on both. From $ \textit{EQ}_2$ to the corresponding saddle node, the strength of the primary central vortex pair decreases and the structure recedes away from highly stratified regions. Conversely, in the case of $\textit{TW}_2^+$ , which is localised towards the upper channel, its vorticity increases towards the saddle node and the concentration becomes more homogenised, especially in the mid-channel where the primary vortices reside.

This picture can be complemented with a look at the energetics. The total energy density $E_{\textit {tot}}$ of each solution was given previously in (2.9). This can be decomposed into parts associated with the streak, roll and wave components of the SSP by substituting $\boldsymbol{u} = \boldsymbol{u}_0 + (\bar {u}-y,0,0)^T + (0,\bar {v},\bar {w})^T + \boldsymbol{u}_w$ with $\boldsymbol{u}_w = \boldsymbol{u} - \bar {\boldsymbol{u}}$ into (2.9) and separating out the terms corresponding to these pieces. This gives

(4.1) \begin{align} E_{\textit {tot}} = E_0 +\frac {1}{2}\int _{-1}^{1}y\langle \overline {\hat {u}}\rangle \mathrm{d}y + K_s + K_r + K_w + \hat {P}, \end{align}

where $E_0 = 1/6 + \textit{Ri}_b(\kappa /v_s - \coth (v_s/\kappa ))$ is the total energy density of the base profile, $K_s$ , $K_r$ , $K_w$ are the streak, roll and wave kinetic energy densities, and $\hat {P}$ is the potential energy density of the concentration perturbation. These constituents are defined by

(4.2a,b) \begin{align} K_s = \frac {1}{2}\int _{-1}^{1}\frac {1}{2} \big \langle \overline {\hat {u}}^2 \big \rangle \,\mathrm{d}y, \quad K_r = \frac {1}{2}\int _{-1}^{1}\frac {1}{2}\langle \bar {v}^2 + \bar {w}^2 \rangle \,\mathrm{d}y, \\[-28pt] \nonumber \end{align}

(4.2c,d) \begin{align} K_w = \frac {1}{2}\int _{-1}^{1} \frac {1}{2}\langle \overline {\boldsymbol{u}_w\boldsymbol{\cdot }\boldsymbol{u}_w}\rangle \,\mathrm{d}y, \quad \hat {P} = \frac {\textit{Ri}_b}{2}\int _{-1}^{1} \big \langle \overline {y\hat {c}}\big \rangle \,\mathrm{d}y . \\[0pt] \nonumber \end{align}

The remaining term in (4.1) (second on the right-hand side) accounts for the net loss in kinetic energy due to the reorganisation of the base shear profile by the ECS and is consequently negative. The perturbation components are plotted in figure 8. We see that the states are largely unaffected by the presence of density gradients as $ \textit{Ri}_b$ increases, until the potential energy of the perturbation reaches the size of the roll kinetic energy. Though this represents a negligible portion of the total flow energy, only one of the three SSP components needs to be fully disrupted to prohibit the existence of these turbulence-supporting solutions. In this case, the rolls feel the energetically penalising effects of stratification through their vertical component, causing them to adapt in the ways already outlined above. This interaction ultimately sets the scale of the maximum $ \textit{Ri}_b$ and mirrors observations in linearly stratified shear flow (Eaves & Caulfield Reference Eaves and Caulfield2015; Olvera & Kerswell Reference Olvera and Kerswell2017). A comparison of figures 8(a) and 8(b) indicates that the $ \textit{EQ}_2$ branch withstands higher $ \textit{Ri}_b$ values than $ \textit{EQ}_1$ because its rolls are an order-of-magnitude stronger. Following some additional continuations at $v_s/\kappa = 1$ , where we explore connections with previous plane Couette flow studies and converge some new states, we shall return to address the question of the maximum bulk Richardson number for other settling velocities in §§ 4.2 and 4.3.

Figure 8.

Energy density components associated with the perturbation fields of the (a) $ \textit{EQ}_1/\textit{TW}_1^+$ and (b) $ \textit{EQ}_2/\textit{TW}_2^+$ solution branches. The labels ‘UB’ and ‘LB’ respectively indicate which parts of each curve correspond to the upper and lower branches in figure 7.

Although $ \textit{Ri}_b$ and $v_s$ are constrained to have the same signs for the physical scenarios of this study, it is nevertheless informative to continue the branches into the $ \textit{Ri}_b \lt 0$ half-plane. We find that $ \textit{EQ}_1$ is formed via a bifurcation at $ \textit{Ri}_b= -2.02\times 10^{-3}$ from a branch of streamwise-invariant roll solutions (dashed red curves in figure 7). This branch emerges from the first eigenmode of the laminar base flow to turn unstable when $ \textit{Ri}_b$ is decreased from zero – a single counter-rotating vortex pair – which occurs at $ \textit{Ri}_b = -9.08\times 10^{-4}$ . Likewise, $\textit{TW}_1^+$ bifurcates off a streamwise-invariant branch of rolls (dashed red) at $ \textit{Ri}_b = -8.91\times 10^{-3}$ , in this case two vortex pairs stacked vertically in the channel, arising from an unstable mode at $ \textit{Ri}_b = -5.53\times 10^{-3}$ . The corresponding linear stability problem is covered in Appendix B. These observations are similar to continuations performed by Olvera & Kerswell (Reference Olvera and Kerswell2017) in stratified plane Couette flow, who observed the $ \textit{EQ}_7$ state bifurcating from a set of sheared Rayleigh–Bénard convection rolls. However, in our case, broken $\mathcal{R}$ symmetry implies that the two-dimensional roll branches can carry non-zero bulk streamwise velocity and the bifurcations off these branches break streamwise invariance to create families of three-dimensional travelling waves with non-zero $a$ . In contrast to $ \textit{EQ}_1$ and $\textit{TW}_1^+$ , the $ \textit{EQ}_2/\textit{TW}_2^+$ family does not immediately appear to connect to a nearby roll branch, at least as far as $ \textit{Ri}_b = -0.055$ (whereupon we terminated our continuations).

Figure 9 shows the travelling wave solutions $\textit{TW}_1^+$ and $\textit{TW}_2^+$ . As already indicated by the corresponding inserts in figure 7, the solutions consist of two pairs of stacked vortices, with the stronger pair of the two residing in the upper channel. These vortices redistribute momentum by advecting faster positive streamwise velocity from the top of the channel downwards and likewise bringing some velocity that is slower than the base flow towards the top. Because this takes place preferentially in the upper channel, the net effect is to skew the bulk velocity of the solutions positive. This is clear from the rightmost panel in figure 9, which shows that $\langle \bar {u} \rangle$ in both cases predominantly exceeds the corresponding value for laminar flow. Furthermore, the two travelling wave velocities are positive, with $a = 0.396$ ( $\textit{TW}_1^+$ ) and $0.465$ ( $\textit{TW}_2^+$ ). However, by the symmetry of the equations at $ \textit{Ri}_b=0$ , the states defined by $\textit{TW}_1^- = \mathcal{R} \boldsymbol{\cdot }\textit{TW}_1^+$ and $\textit{TW}_2^- = \mathcal{R} \boldsymbol{\cdot }\textit{TW}_2^+$ are also exact solutions, with corresponding velocities $a = -0.396$ and $-0.465$ . Since their primary vortices lie in the lower portion of the channel, where sediment gradients are higher, these negative travelling waves respond differently to variations in $ \textit{Ri}_b$ .

Figure 9.

Travelling wave solutions $\textit{TW}_1^+$ and $\textit{TW}_2^+$ at $ \textit{Ri}_b = 0$ . The first two panels show filled contours of the streamwise velocity perturbation $\bar {u} - y$ , coloured from brown (negative) to blue (positive), using intervals of size $0.05$ centred around zero (white). Additionally, streamlines of $(\bar {v},\bar {w})$ are overlaid and coloured according to $|(\bar {v},\bar {w})|$ , from blue (zero) to red. The red maximum value is $0.023$ for $\textit{TW}_1^+$ and $0.049$ for $\textit{TW}_2^+$ (2 significant figures). The rightmost panel plots the streamwise velocity after averaging in $x$ and $z$ for $\textit{TW}_1^+$ (solid blue) and $\textit{TW}_2^+$ (solid green). The dashed grey line is the streamwise velocity of laminar flow, i.e. $\langle \bar {u} \rangle = y$ .

Tracing out these related solution families uncovers a richer bifurcation structure. In figure 10(a) we plot the result of continuing $\textit{TW}_1^-$ and $\textit{TW}_2^-$ in $ \textit{Ri}_b$ (dotted curves), alongside the $ \textit{EQ}_1$ / $\textit{TW}_1^+$ and $ \textit{EQ}_2$ / $\textit{TW}_2^+$ branches already shown in figure 7. For reference, the corresponding streamwise wave speeds for each of the branches are plotted directly below in figure 10(b), using the same horizontal axes. This highlights that the symmetry of the positive and negative travelling wave branches breaks as soon as $ \textit{Ri}_b$ deviates from zero. For example, unlike the $\textit{TW}_1^+$ and $\textit{TW}_2^+$ branches, the $\textit{TW}_1^-$ and $\textit{TW}_2^-$ families undergo a sign change of $a$ at positive $ \textit{Ri}_b$ , passing through equilibria at $ \textit{Ri}_b = 2.3\times 10^{-3}$ and $ \textit{Ri}_b = 5.9\times 10^{-3}$ (3 significant figures). Both branches form a saddle node in $ \textit{Ri}_b$ , attaining a maximum $ \textit{Ri}_b$ and connecting back to other equilibrium solutions at $ \textit{Ri}_b = 0$ previously identified by Gibson et al. (Reference Gibson, Halcrow and Cvitanović2009). Specifically, $\textit{TW}_1^-$ connects to $ \textit{EQ}_3$ and $\textit{TW}_2^-$ connects to $ \textit{EQ}_4$ . Intriguingly, $ \textit{EQ}_3$ and $ \textit{EQ}_4$ are themselves connected, because they form a saddle-node pair when unstratified plane Couette flow is continued in $ \textit{Re}$ (Gibson et al. Reference Gibson, Halcrow and Cvitanović2009), possibly hinting at a deeper relationship between this pair and $ \textit{EQ}_1$ / $ \textit{EQ}_2$ .

Figure 10.

Parametric continuation of the $\textit{TW}_1^-$ (dotted blue) and $\textit{TW}_2^-$ (dotted green) families and related states (solid lines). States are continued in $ \textit{Ri}_b$ and plotted as a function of their (a) mean wall-normal stress perturbation $\hat \tau _y$ , and (b) streamwise wave speed $a$ . The dashed red curves in panel (a) are two-dimensional roll branches also shown in figure 7. In both plots, each labelled state is indicated with a filled circle.

When continued into negative $ \textit{Ri}_b$ , the $ \textit{EQ}_4/\textit{TW}_2^-$ branches do not obviously undergo successive bifurcations, at least within the range of our computations ( $ \textit{Ri}_b\gt -0.055$ ). The $\textit{TW}_1^-$ branch appears to asymptotically approach one of the two-dimensional roll branches as $ \textit{Ri}_b$ decreases, though the branches remain disconnected. On continuing $ \textit{EQ}_3$ , we find that it connects to a family of solutions (solid orange) featuring vortex pairs stacked vertically in the channel. The new solution branch possesses the additional half-cell shift symmetry

(4.3) \begin{align} \mathcal{T}:[u,v,w,c](x,y,z) \mapsto [u,v,w,c](x + L_x/2, y, z + L_z/2) \end{align}

that, when combined with $\mathcal{S}$ implies that states are mirror symmetric in $z$ . The connection with $ \textit{EQ}_3$ occurs when this symmetry is broken in a pitchfork bifurcation that gives rise to $ \textit{EQ}_3$ and its spanwise reflected counterpart. Since spanwise reflections do not carry implications for the stratification, we do not distinguish between such states.

The $\mathcal{T}$ -symmetric branch leads to yet further equilibrium and travelling wave solutions when continued in $ \textit{Ri}_b$ . We plot the full bifurcation diagram obtained in figure 11 and summarise the structures of all the solutions identified and the connections between them in figure 12. Though complicated, we can systematically trace these families of exact solutions and the exposition that follows is illustrated by both these figures in concert. The new branch emerges from a bifurcation off another two-dimensional roll family at $ \textit{Ri}_b = -0.0215$ (3 significant figures) (the $(m,n)=(0,2)$ secondary mode, in the nomenclature of Appendix B linear stability analysis) and contains the $ \textit{EQ}_7$ equilibrium – a solution at $ \textit{Ri}_b=0$ originally identified by Gibson et al. (Reference Gibson, Halcrow and Cvitanović2009) and Itano & Generalis (Reference Itano and Generalis2009). Following the branch further, it reaches substantially higher $ \textit{Ri}_b$ values than the solutions in figure 10 and weaves either side of $ \textit{Ri}_b=0$ multiple times. The first two axis crossings give rise to two new travelling wave solutions with $a\gt 0$ , which we call $\textit{TW}_7^+$ and $\textit{TW}_{7'}^+$ . Since these appear to be undiscovered homogeneous states for this computational cell, we report their properties in table 1. They possess particularly high wall stress values and a greater number of unstable eigenvalues than any of the solutions presented by Gibson et al. (Reference Gibson, Halcrow and Cvitanović2009) or Ahmed & Sharma (Reference Ahmed and Sharma2020), which may point to why they have not been found in these prior studies. Both the new states can be transformed, via $\mathcal{R}$ , to give their negatively directed counterparts $\textit{TW}_7^-$ and $\textit{TW}_{7'}^-$ . Continuing these states reveals that they are connected by an isola.

Table 1.

Properties of the new travelling wave states when $ \textit{Ri}_b = 0$ . The table columns are as follows: $||{\boldsymbol{\cdot }}||$ is the $L^2$ norm of the velocity perturbation $(\hat u, \hat v, \hat w)$ , $E_{\textit {tot}}$ is the energy density (2.9), $D$ is the dissipation (2.11b ), $\dim W^u$ is the dimension of the unstable manifold, $\dim W^u_H$ is the number of unstable eigenvectors that lie within the symmetry group $H$ of the solution (in this case, the group generated by $\mathcal{S}$ and $\mathcal{T}$ ), $a$ is the wave speed and the final column reports the bulk velocity.

Figure 11.

Bifurcation diagram for all solution branches identified at $v_s/\kappa = 1$ . Curves previously shown in figures 7 and 10 are plotted together with the additional $\textit{TW}_7^+$ / $\textit{TW}_{7'}^+$ (solid orange), $\textit{TW}_7^-$ / $\textit{TW}_{7'}^-$ (dotted orange), $ \textit{EQ}_8/\textit{TW}_8^+$ (solid purple) and $\textit{TW}_8^-$ (dotted purple) families. States are continued in $ \textit{Ri}_b$ and plotted as a function of their (a) mean wall-normal stress perturbation $\hat \tau _y$ , and (b) streamwise wave speed $a$ . The dashed red curves are two-dimensional roll solutions that bifurcate off from linear instabilities of the base flow. In addition to those identified in figure 7, we include a branch that emerges at $ \textit{Ri}_b = -0.0136$ (3 significant figures). In both plots, each labelled state is indicated with a filled circle.

Figure 12.

Diagram of relationships between states at $v_s/\kappa = 1$ . The first and third columns show named states at $ \textit{Ri}_b=0$ that are connected by continuation in the $ \textit{Ri}_b \gt 0$ half-plane. Connections are indicated by lines styled according to the corresponding curves plotted in figure 11. Intermediate states plotted in the middle column are located at the maximum $ \textit{Ri}_b$ attained along the relevant branch. Connections via $ \textit{Ri}_b \lt 0$ are indicated in the margins. States related by the action of $\mathcal{R}$ are joined by a dashed grey line. Each plot shows up to $15$ contours of $\bar {u} - y$ equispaced between $\pm 0.75$ and streamlines of $(\bar {v},\bar {w})$ , coloured in the same style as in figures 1 and 9.

Following the $\textit{TW}_7^+$ / $\textit{TW}_{7'}^+$ branch from $\textit{TW}_{7'}^+$ towards higher $ \textit{Ri}_b$ , it reaches a maximum at $ \textit{Ri}_b = 0.123$ (3 significant figures) before turning and approaching a structural feature suggestive of a pitchfork bifurcation at $ \textit{Ri}_b \approx 0.1$ . This is highlighted in the inset panel of figure 11. Although the continuation curve is continuous through $a = 0$ here and does not connect to other branches, we opt to restyle the solution curve dotted purple where $a\lt 0$ . As shown in figure 13, the yellow ‘central fork’ of the branch is approximately $\mathcal{R}$ symmetric as it approaches $a=0$ . Since the branch is also $\mathcal{T}$ symmetric, this implies that the solutions become close to being symmetric with respect to rotation by $\pi$ around the $z$ axis, as highlighted in the side panels of figure 13. (However, note that even at $a=0$ , where the state becomes an equilibrium, this symmetry cannot be exact, since $\mathcal{R}$ is not a symmetry of the governing equations when both $v_s$ and $ \textit{Ri}_b$ are non-zero.) As the branch continues to $a\lt 0$ , the symmetry becomes more obviously broken and this section of the branch connects through to a new negatively directed travelling wave state, which we label $\textit{TW}_8^-$ (see also figure 11). The positively directed $\textit{TW}_8^+ = \mathcal{R}\boldsymbol{\cdot }\textit{TW}_8^-$ yields a separate family that forms the other half of the disconnected pitchfork. Bulk statistics for this new solution pair are reported in table 1. As demonstrated by the contour plots in figure 13, the ‘upper’ and ‘lower’ forks are approximately related by the $\mathcal{R}$ symmetry, lending credence to the idea that this structural feature is in essence an $\mathcal{R}$ -symmetry breaking pitchfork bifurcation, fractured into disconnected pieces due to the presence of asymmetric stratification. On following the $\textit{TW}_8^+$ branch further, we note that it ultimately connects to $ \textit{EQ}_8$ , which is the upper branch counterpart of $ \textit{EQ}_7$ when these equilibria are continued in $ \textit{Re}$ (Gibson et al. Reference Gibson, Halcrow and Cvitanović2009).

Figure 13.

Close-up view of the $\textit{TW}_{7'}^+$ / $\textit{TW}_8^-$ and $\textit{TW}_8^+$ / $ \textit{EQ}_8$ connections. The left-hand plot shows a portion of the continuation curves from figure 11(b), with the $\textit{TW}_7^-$ / $\textit{TW}_{7'}^-$ branches removed for clarity. The right-hand side panels shows contours of $u(0, y, z)$ for the states at the locations indicated with the corresponding labels in the bifurcation diagram. For $\mathcal{R}$ - and $\mathcal{T}$ -symmetric states, such slices must be anti-symmetric with respect to reflection about the line $y=0$ (dashed grey). This property is approximately satisfied for the middle two panels, along the ‘central fork’ and clearly broken for the ‘upper’ and ‘lower’ forks.

The high wall stress families emerging from $\textit{TW}_7^\pm$ , $\textit{TW}_{7'}^\pm$ and $\textit{TW}_8^\pm$ have a noticeably different effect on the concentration field when compared with the states with a lower stress such as $ \textit{EQ}_1$ and $ \textit{EQ}_7$ . To demonstrate this, we compare states at $ \textit{Ri}_b=0$ in figure 14. By way of illustration, panel (a) shows slices of $\bar {c}(y,z)$ for the relatively low stress $ \textit{EQ}_8$ state versus the high stress travelling wave $\textit{TW}_8^+$ . In the latter case, the concentration field is relatively homogeneous in the channel interior. The corresponding implication for the sediment transport is indicated in figure 14(b,c), which respectively show mean streamwise and wall-normal fluxes for both states. We see that the streamwise transport of $\textit{TW}_8^+$ is higher throughout most of the channel, facilitated by its stronger vortices, but it is far from homogeneous – fluxes are elevated near the high-speed streaks where the sediment is advected more rapidly. In figure 14(d–f) we plot cross-stream averaged profiles for all of the $ \textit{Ri}_b=0$ solutions considered herein, dividing the data into low and high stress groups. To the high stress category, which includes our new travelling wave solutions, we add $ \textit{EQ}_2$ , since its excess stress is more than double that of the other equilibria (see figure 11). (Other states from the study of Gibson et al. Reference Gibson, Halcrow and Cvitanović2009 that meet this criterion would include $ \textit{EQ}_{11}$ and $ \textit{EQ}_{13}$ .) The data are clearly separated in a manner that mirrors our example contours given in figure 14(a–c), with the high stress states exhibiting a more homogeneous concentration field and greater streamwise and wall-normal fluxes. For small increments of $ \textit{Ri}_b$ above zero, these trends persist. Eventually though, many of the lower stress states cease to exist, making it difficult to compare the solutions collectively (figure 11). However, the fact that each of the high stress travelling wave branches can be continued far higher in $ \textit{Ri}_b$ than the rest is notable and is perhaps due to the homogenisation of the $c$ field, which lessens the effects of stratification in the channel interior.

Figure 14.

Transport fluxes for low and high stress states with $ \textit{Ri}_b=0$ . (a–c) Contours of (a) $\bar {c}$ , (b) $\overline {(u+1)c}$ and (c) $\overline {vc}$ for $ \textit{EQ}_8$ (left-hand panels) and $\textit{TW}_8^+$ (right-hand panels). (d–f) Wall-normal profiles of (d) $\langle \bar {c}\rangle$ , (e) $\langle \overline {(u+1)c}\rangle$ and (f) $\langle \overline {vc}\rangle$ for the low stress states (dark grey dashed) $ \textit{EQ}_1$ (squares), $ \textit{EQ}_3$ (upward-pointing triangles), $ \textit{EQ}_4$ (circles), $ \textit{EQ}_7$ (crosses), $ \textit{EQ}_8$ (downward-pointing triangles), and for the high stress states (red dashed) $ \textit{EQ}_2$ (stars), $\textit{TW}_7^\pm$ (diamonds), $\textit{TW}_{7'}^\pm$ (pluses), $\textit{TW}_8^\pm$ (pentagons).

4.2. Solutions with $v_s/\kappa \lt 1$

We now turn our attention to the effects of stratification on states away from the illustrative case of $v_s/\kappa = 1$ . To keep our numerical computations manageable and to aid their reporting, we return to focusing on the $ \textit{EQ}_1/\textit{TW}_1^\pm$ and $ \textit{EQ}_2/\textit{TW}_2^\pm$ families, with the anticipation that these states are representative of the broader class of equilibria and travelling waves underpinned by the SSP/VWI mechanism.

On decreasing the settling velocity, the concentration fields for each state are expected to become increasingly homogeneous, since they must eventually become uniform in the limit $v_s/\kappa \to 0$ , regardless of the magnitude of $ \textit{Ri}_b$ . In § 3.1 we used this to determine that in the passive regime ( $ \textit{Ri}_b=0$ ) the deviation of $c$ away from a homogeneous state is $O(v_s/\kappa )$ . As $ \textit{Ri}_b$ is increased away from this limiting case, we can ask how large it must become in order for the effect of sediment-induced stratification to be non-negligible. Since $\boldsymbol{u}$ and $c$ remain uncoupled to leading order until this regime is reached, the magnitude of the density perturbation in the passive case implies that the perturbation to the buoyancy force is $O(\textit{Ri}_b v_s/\kappa )$ , while $ \textit{Ri}_b$ remains sufficiently small. This reaches the size of the $O(\textit{Re} ^{-2})$ terms in the vertical momentum equation when $ \textit{Ri}_b = O(\kappa v_s^{-1}\textit{Re} ^{-2})$ . If $ \textit{Ri}_b$ is increased beyond this regime, solutions that obey the VWI balances must either cease to exist or adapt to the presence of stratification by modifying their structure to match the larger buoyancy forcing term. This can happen by developing inner scales that enhance the spatial gradients of the flow fields, i.e. localisation (Olvera & Kerswell Reference Olvera and Kerswell2017).

Figure 15.

Collapse of solution branches related to $ \textit{EQ}_1$ and $ \textit{EQ}_2$ , as $v_s\to 0$ . The two rows of plots show continuations in $ \textit{Ri}_b$ versus the (a,b) mean wall stress perturbation $\hat {\tau }_y$ and (c,d) streamwise wave velocity $a$ , with the horizontal axes rescaled by $v_s$ to demonstrate collapse of the data.

Within the minimal flow cell and moderate Reynolds number considered herein, this latter possibility does not occur. Instead, the states exist up to a maximum $ \textit{Ri}_b$ that is proportional to $1/v_s$ . This is evidenced in figure 15, in which continuations of our two example solution families in $ \textit{Ri}_b$ are plotted for $v_s/\kappa$ decreasing from $1$ down to $0.01$ . We see that they collapse almost exactly onto a set of asymptotic curves when rescaled under the combination $ \textit{Ri}_b v_s$ . Furthermore, the continuation clarifies the relationship between the equilibria and travelling wave states. The symmetry-related pairs $\textit{TW}_1^\pm$ and $\textit{TW}_2^\pm$ are created in subcritical pitchfork bifurcations that are asymptotically realised as $v_s\to 0$ , where the system approaches invariance under $\mathcal{R}$ . Since the equations decouple in this limit, the $ \textit{Ri}_b$ at which the bifurcation occurs is sent to infinity. The equilibrium pairs $ \textit{EQ}_1$ / $ \textit{EQ}_3$ and $ \textit{EQ}_2$ / $ \textit{EQ}_4$ are (asymptotically) connected either side of their respective pitchforks, which may explain why the computations of Gibson et al. (Reference Gibson, Halcrow and Cvitanović2009) showed that the portions of $ \textit{EQ}_3$ and $ \textit{EQ}_4$ ’s unstable manifolds that lie within the solutions’ symmetry groups are exactly one dimension larger than those corresponding to $ \textit{EQ}_1$ and $ \textit{EQ}_2$ , respectively.

Figure 16.

Sequence of continuations in $ \textit{Ri}_b$ for increasing $v_s/\kappa$ (as labelled). (a) Continuation of $ \textit{EQ}_1$ and $\textit{TW}_1^+$ (solid lines, first four panels), which become disconnected by $v_s/\kappa = 4$ to leave separate branches. In the final panel, the $ \textit{EQ}_1$ branch is shown as a solid line and ultimately connects to $ \textit{EQ}_2$ (not shown) and the $\textit{TW}_1^+$ branch is dotted and connects to $\textit{TW}_2^+$ (not shown). (b) Continuation of $ \textit{EQ}_3$ and $\textit{TW}_1^-$ (dotted lines, first three panels), which likewise become disconnected by $v_s/\kappa = 2$ , with $ \textit{EQ}_3$ connecting to $ \textit{EQ}_4$ (solid lines) and $\textit{TW}_1^-$ connecting to $\textit{TW}_2^-$ (dotted lines).

Figure 17.

Continuation of the $ \textit{EQ}_1$ / $ \textit{EQ}_2$ branch in $ \textit{Ri}_b$ for high settling velocities. The solution curves are coloured from light orange to black depending on $v_s/\kappa$ . Panel (b) demonstrates their collapse under rescaling $ \textit{Ri}_b$ by $(v_s/\kappa )\exp (-v_s/\kappa )$ .

4.3. Solutions with $v_s/\kappa \gt 1$

Increasing $v_s/\kappa$ beyond unity results in some structural changes to the connections plotted in figures 7 and 15. These are shown in figure 16. Firstly, we see in row (a) that between $v_s/\kappa = 3.7$ and $v_s/\kappa = 4$ , the link between $ \textit{EQ}_1$ and $\textit{TW}_1^+$ becomes severed. At higher settling velocities, $ \textit{EQ}_1$ instead connects to $ \textit{EQ}_2$ (as it does under $ \textit{Re}$ continuation) and $\textit{TW}_1^+$ connects to $\textit{TW}_2^+$ . Likewise, figure 16(b) shows that between $v_s/\kappa = 1.75$ and $2$ , the $ \textit{EQ}_3/\textit{TW}_1^-$ and $ \textit{EQ}_4/\textit{TW}_2^-$ connections ‘swap partners’, leaving $ \textit{EQ}_3$ / $ \textit{EQ}_4$ and $\textit{TW}_1^-$ / $\textit{TW}_2^-$ as saddle-node pairs in $ \textit{Ri}_b$ at higher $v_s/\kappa$ .

After these reorganisations at intermediate settling velocities, further increases in $v_s/\kappa$ lead the solution branches towards a high settling regime, in which states persist up to much higher $ \textit{Ri}_b$ values. Focusing on the $ \textit{EQ}_1/\textit{EQ}_2$ branch, further continuations are shown figure 17, up to $v_s/\kappa = 14$ . Panel (a) of this figure demonstrates the trend, while panel (b) shows collapse of the curves under an exponential rescaling that shall be rationalised later. For example, denoting the largest bulk Richardson number attained along the relevant solution branch by $ \textit{Ri}_{b,\textit {max}}$ , there is an order-of-magnitude increase in this quantity between the $v_s/\kappa = 6$ and $8$ cases, from $ \textit{Ri}_{b,\textit {max}} =$ $0.028$ to $0.14$ , respectively. Our prior analysis of the passive scalar regime points to why this might occur. As observed in figure 2, increasing $v_s/\kappa$ from unity when $ \textit{Ri}_b = 0$ causes the sediment field to concentrate at the lower wall and consequently dilutes the remainder of the channel. This suggests that, for fixed $ \textit{Ri}_b\in (0,\textit{Ri}_{b,\textit {max}})$ and increasing $v_s/\kappa$ , the coupled interaction between the concentration field and the velocity structure must ebb away, as the buoyancy forces in the bulk of the channel become weaker due to reduced sediment load. This is demonstrated in figures 18 and 19, in which lower branch states at $ \textit{Ri}_b=0.002$ are plotted for increasing $v_s/\kappa$ . Relative to the passive case, which is also depicted in the two figures, the kinetic energy of the velocity perturbation (figure 18) is preferentially located away from the highly concentrated boundary layers (figure 19). As $v_s/\kappa$ increases, the width of these layers diminishes, reducing the effect of the sediment on the flow fields.

Figure 18.

Effect of stratification as $v_s/\kappa$ is increased along the lower branch. In each panel, seven filled contours of the streamwise-averaged perturbation kinetic energy $( {1}/{2})(\overline {\hat {u}^2+\hat {v}^2+\hat {w}^2})$ are plotted, equispaced between $0$ and $0.105$ . The leftmost panel shows $ \textit{EQ}_1$ at $ \textit{Ri}_b = 0$ and the trio of right-hand panels show the corresponding states obtained by continuing to $ \textit{Ri}_b = 0.002$ , with $v_s/\kappa = 4$ , $8$ and $12$ , as labelled. Grey dashed lines indicate the distance of the first contour from the lower wall in each panel.

Figure 19.

Streamwise-averaged concentration fields corresponding to the lower branch states in figure 18. Each panel shows seven filled contours of $\min (\overline {c}, 15)$ . The applied thresholding removes very large values from the rightmost two plots below $y = -0.95$ (which are not expected to affect the velocity fields significantly) so that the same contour scale can be used for each panel. The grey dashed lines lie at the same locations as those in figure 18.

Figure 20.

Plots of the  $ \textit{EQ}_1$ / $ \textit{EQ}_2$ branch states at the saddle node, where $ \textit{Ri}_b = \textit{Ri}_{b,\textit {max}}$ . Each panel shows streamwise-averaged perturbation kinetic energy contours, coloured in the same style as figure 18, with the leftmost plot showing $ \textit{EQ}_1$ at $ \textit{Ri}_b=0$ and the subsequent plots showing $v_s/\kappa = 4$ , $8$ and $12$ . Grey dashed lines indicate the distance of the first contour from the lower wall in each panel.

Conversely, for fixed $v_s/\kappa$ and increasing $ \textit{Ri}_b$ , the vertical location for which an equivalent buoyancy effect is felt increases, since this is proportional to $ \textit{Ri}_b$ . This enlarges the effective width of the stratified lower region and ultimately drives states away from the bottom wall. At $ \textit{Ri}_{b,\textit {max}}$ , the velocity fields approach a limiting solution as $v_s/\kappa$ becomes large. This is summarised in figure 20, where we plot the kinetic energy of the velocity perturbation at $ \textit{Ri}_{b,\textit {max}}$ for $v_s/\kappa = 4,8,12$ , which is seen to be concentrated within the portion of the channel above $y\approx -0.6$ for all three solutions. Presumably, beyond this point the vortices are no longer able to adapt to the presence of stratification by localising towards the upper part of the channel and the solution branch must turn back on itself, as shown in figure 17.

Compared with the low settling case (§ 4.2), the question of what the asymptotic structures should be in the $v_s/\kappa \gg 1$ regime is less clear cut. On increasing $ \textit{Ri}_b$ from zero with $v_s/\kappa \gg 1$ fixed, the buoyancy force becomes non-negligible when $ \textit{Ri}_b\hat {c}$ is comparable to the dominant balances in the vertical momentum equation. From our analysis of the passive case in § 3.2, this occurs first near the bottom wall, where $\hat {c}$ is $O(1)$ and the leading physics is controlled by vertical diffusion of the $O(\textit{Re} ^{-1})$ velocity field. Therefore, the flow momentum should begin to feel the presence of the sediment when $ \textit{Ri}_b = O(\textit{Re} ^{-2})$ . However, we have seen that states cannot be fully disrupted by near-wall stratification alone, so this cannot be the scale that determines $ \textit{Ri}_{b,\textit {max}}$ . To identify the point beyond this scale, at which the velocity structure in the bulk of the channel becomes disrupted by stratification, we adopt an empirical approach. Figure 21(a) shows vertical profiles of the size of the buoyancy perturbation (measured by $ \textit{Ri}_b\max _{x,z}|\hat {c}|$ ) halfway along the lower branch for the $v_s/\kappa = 8,10,12,14$ curves, alongside the equivalent data for the contribution of the base concentration, $ \textit{Ri}_b c_0$ . We see that, at least by $v_s/\kappa = 10$ , the buoyancy perturbation curves have a similar magnitude and dependence in the lower half of the channel and become considerably weaker in the upper half, as $v_s/\kappa$ increases. The picture is qualitatively similar for states at $ \textit{Ri}_{b,\textit {max}}$ , whose profiles are plotted in figure 21(b) and possess higher values overall, due to increased $ \textit{Ri}_b$ . Based on these curves, we speculate that the asymptotic velocity fields drive a concentration perturbation that collapses under an appropriate rescaling in the lower half of the channel, but tends towards vanishing concentration in the upper half as $v_s/\kappa \to \infty$ . We note that, since the corresponding $c_0(y)$ curves all approximately intersect at a single point $y=y^*$ , the rescaling factor can be measured by leveraging the fact that $ \textit{Ri}_b c_0(y^*) \approx \mathrm{const}$ at this point. This gives $ \textit{Ri}_b \sim 1/c_0(y^*)\sim ( {\kappa }/{v_s})\mathrm{e}^{{( {v_s}/{\kappa })(y^*+1)}}$ as the scale that controls the structure of the high settling velocity states. We estimate $y^*\approx 0$ from the figure. As demonstrated in figure 17(b), this is sufficient to collapse the available data reasonably well.

Figure 21.

Decay of the buoyancy term for high $v_s/\kappa$ states (a) halfway along the $ \textit{EQ}_1$ lower branch, and (b) at the saddle node ( $ \textit{Ri}_b=\textit{Ri}_{b,\textit {max}}$ ). Dashed curves plot $ \textit{Ri}_b c_0(y)$ for the corresponding cases.

4.4. Summary

Taken together, our investigations of the $ \textit{EQ}_1/\textit{EQ}_2$ family in both extremal regimes of settling velocity imply that $ \textit{Ri}_{b,\textit {max}}$ is non-monotonic with respect to $v_s/\kappa$ , because $ \textit{Ri}_{b,\textit {max}}$ becomes unboundedly large as either $v_s/\kappa \to 0$ or $v_s/\kappa \to \infty$ . This echoes our findings in the passive scalar regime, where it was seen that either limit causes the concentration field to homogenise, leading to non-monotonic sediment transport curves. We summarise our continuation data in the stratified regime, by plotting $ \textit{Ri}_{b,\textit {max}}$ for each computed $v_s/\kappa$ in figure 22. Both limits are seen to agree well with the asymptotic scalings deduced in §§ 4.2 and 4.3 (grey curves). Other states that are underpinned by the SSP/VWI mechanism in the passive scalar regime can reasonably expect to scale similarly, at least in the low settling regime. In the high settling regime, where geometric considerations appear to play a role in constraining the interaction between the sediment and velocity fields, it may be the case that solutions with different vortex morphologies, such as $ \textit{EQ}_7$ , behave quantitatively differently. Nevertheless, the scaling regimes demonstrated carry implications for the range of bulk Richardson numbers that support turbulent flow and may be compared with prior DNS explorations of the laminar–turbulent boundary, which are quantitatively similar to figure 22 (Cantero et al. Reference Cantero, Shringarpure and Balachandar2012b ; Langham & Hogg Reference Langham and Hogg2024). This is discussed in greater detail in the following section. We note additionally that the minimum of the $ \textit{Ri}_{b,\textit {max}}$ curve occurs at $v_s/\kappa \approx 2.5$ , which is close to the point at which vertical and streamwise transport of the concentration perturbation is maximised for the passive case (figure 6).

Figure 22.

Maximum bulk Richardson number attained by the $ \textit{EQ}_1$ / $ \textit{EQ}_2$ family of states. The fitted grey curves are $B_1 (v_s/\kappa )^{-1}$ and $B_2 (v_s/\kappa )^{-1}\exp (v_s / \kappa )$ with $B_1 = 0.0125$ and $B_2 = 3.4\times 10^{-4}$ .