Structured input–output analysis of stably stratified plane Couette flow

Abstract

We employ a recently introduced structured input–output analysis (SIOA) approach to analyse streamwise and spanwise wavelengths of flow structures in stably stratified plane Couette flow. In the low-Reynolds-number ($Re$) low-bulk Richardson number ($Ri_b$) spatially intermittent regime, we demonstrate that SIOA predicts high amplification associated with wavelengths corresponding to the characteristic oblique turbulent bands in this regime. SIOA also identifies quasi-horizontal flow structures resembling the turbulent–laminar layers commonly observed in the high-$Re$ high-$Ri_b$ intermittent regime. An SIOA across a range of $Ri_b$ and $Re$ values suggests that the classical Miles–Howard stability criterion ($Ri_b\leq 1/4$) is associated with a change in the most amplified flow structures when the Prandtl number is close to one ($Pr\approx 1$). However, for $Pr\ll 1$, the most amplified flow structures are determined by the product $PrRi_b$. For $Pr\gg 1$, SIOA identifies another quasi-horizontal flow structure that we show is principally associated with density perturbations. We further demonstrate the dominance of this density-associated flow structure in the high $Pr$ limit by constructing analytical scaling arguments for the amplification in terms of $Re$ and $Pr$ under the assumptions of unstratified flow (with $Ri_b=0$) and streamwise invariance.

1. Introduction

Statically stable density stratification in wall-bounded shear flows plays an important role in many industrial and environmental applications, e.g. in cooling equipment (Zonta & Soldati 2018), and the turbulent boundary layers governing atmospheric and oceanic flows (Pedlosky 2013; Vallis 2017). In the atmospheric boundary layer, stable stratification arising from strong ground cooling effects is of particular importance at night (Nieuwstadt 1984; Mahrt 1999, 2014) and near the polar region (Grachev et al. 2005). At the ocean floor, stable density stratification is known to influence the boundary layer thickness (Weatherly & Martin 1978; Lien & Sanford 2004).

(Stably) stratified plane Couette flow (PCF) is a canonical model for stratified wall-bounded shear flow. When the density as well as the velocity is maintained at different values at the two horizontal boundary planes, with gravity acting vertically, stratified PCF has the added benefit (as defined more precisely below) that a natural bulk Richardson number $Ri_b$ can be defined, thereby enabling the relative significance of the imposed stratification and shear to be captured. Furthermore, unstratified PCF has no linear instability for any Reynolds number ($Re$, again defined more precisely below) (Romanov 1973), and yet is observed to transition at Reynolds numbers as low as $Re=360\pm 10$ (Tillmark & Alfredsson 1992). Stratified PCF is a convenient model flow for investigating the effect of stable stratification on the transition dynamics (Deusebio, Caulfield & Taylor 2015).

Stable stratification provides a restoring buoyancy force inhibiting vertical motion (Turner 1979; Davidson 2013). Thus, transition to turbulence in stably stratified PCF typically occurs at a higher Reynolds number than unstratified PCF; see e.g. Deusebio et al. (2015); Eaves & Caulfield (2015); Deguchi (2017) and Olvera & Kerswell (2017). In the transitional regime, both stratified PCF and unstratified PCF exhibit spatial intermittency; i.e. the coexistence of laminar and turbulent regions. In the relatively low-$Re$ low-$Ri_b$ intermittent regime, the spatial intermittency in stratified PCF is characterized by oblique turbulent bands (Deusebio et al. 2015; Taylor et al. 2016) at least qualitatively similar to those seen in unstratified PCF (Prigent et al. 2003; Duguet, Schlatter & Henningson 2010) with a very large channel size (${\sim }O(100)$ times the channel half-height). In the high-$Re$ high-$Ri_b$ intermittent regime, flow structures are instead characterized by turbulent and laminar layers over the vertical direction due to the strong effect of buoyancy (Deusebio et al. 2015). This spatial intermittency directly imposes challenges for the computation of averaged measurements of flow behaviour (such as the efficiency of mixing or the dissipation rate), and thus an understanding of the underlying mechanisms is important for the parameterization of turbulence properties, in particular the irreversible mixing in stratified flows (Caulfield 2020, 2021).

The existence of a unique critical Richardson number that separates flow into laminar and turbulent regimes is questionable, to put it mildly (Andreas 2002; Galperin, Sukoriansky & Anderson 2007). A threshold value close to $1/4$ is supported by some field measurements (Kundu & Beardsley 1991) and experiments (Rohr et al. 1988), although other field measurements reported sustained turbulence in flows with Richardson numbers ${\simeq }1$ (Lyons, Panofsky & Wollaston 1964). More recently, increasing evidence has been found that vertically sheared stably stratified flow appears to self-organize to maintain an appropriately defined Richardson number near $1/4$, both in field observations (Smyth & Moum 2013; Smyth, Nash & Moum 2019) and in simulations (Salehipour, Peltier & Caulfield 2018). This threshold value of $1/4$ also appears in the classical ‘Miles–Howard’ theorem (Howard 1961; Miles 1961), which provides a necessary condition for linear instability in inviscid, non-diffusive steady parallel flow. In particular, it states that instability requires that the local or gradient Richardson number be less than $1/4$ somewhere. Therefore, it is of interest to consider stratified PCF as a well-controlled sheared stratified flow to investigate whether some kind of self-organized criticality and/or marginal stability naturally emerges in a viscous and diffusive flow.

The Prandtl number ($Pr=\nu /\kappa$, where $\nu$ is the kinematic viscosity and $\kappa$ is the diffusivity of the density field) plays a perhaps unsurprisingly important role in determining flow structures. For example, for sufficiently small $Pr$, flows with the same value of the product $Pr Ri_b$ develop the same averaged vertical density profile (Langham, Eaves & Kerswell 2020). This observation that the product $PrRi_b$ determines flow behaviour in the low Prandtl number limit is widely observed in stratified shear flows; see e.g. Lignieres (1999); Garaud, Gallet & Bischoff (2015); Garaud, Gagnier & Verhoeven (2017) and Garaud (2021). Conversely, in the high Prandtl number regime, studies of exact coherent structures in stratified PCF (Langham et al. 2020) show that a nearly uniform density region forms near the channel centre, and the influence of bulk Richardson number on the averaged properties of the flow are significantly reduced. Moreover, Taylor & Zhou (2017) proposed a multi-parameter criterion for the formation of a ‘staircase’ in the density distribution (i.e. a distribution with relatively deep ‘layers’ of nearly uniform density separated by relatively thin interfaces of enhanced density gradient) which suggests that this staircase formation is actually favoured for a large Prandtl number (Taylor & Zhou 2017). The sharpness of the density interfaces also appears to increase as the Prandtl number increases (Zhou et al. 2017b). In addition, increasing the Prandtl number has a larger influence on the mean density profiles than on the mean velocity profiles (Zhou, Taylor & Caulfield 2017a).

The oblique turbulent bands observed in the intermittent regime of stratified PCF (Deusebio et al. 2015; Taylor et al. 2016) require a very large channel size to accommodate them fully, which poses challenges for both simulations and experiments. The three different flow parameters of interest, $Re$, $Ri_b$ and $Pr$ also lead to computational challenges in exploring the full range of flow regimes. To overcome these challenges to direct numerical simulation (DNS), we use an input–output (resolvent) analysis based approach. Such methods, built upon the spatio-temporal frequency response, have been widely employed in unstratified wall-bounded shear flows (Farrell & Ioannou 1993a; Bamieh & Dahleh 2001; Jovanović & Bamieh 2005; McKeon & Sharma 2010; McKeon 2017). This analysis framework has the advantage of computational tractability and is not subject to finite channel effects. Related analysis has shown promise in studying stratified flows, including inviscid stratified shear flow with constant shear (Farrell & Ioannou 1993b), stratified PCF (Jose et al. 2015, 2018) and stratified turbulent channel flow (Ahmed et al. 2021).

In this work, we extend the structured input–output analysis (SIOA) originally developed for unstratified PCF (Liu & Gayme 2021) to stratified PCF. Prior application of the SIOA approach to unstratified transitional wall-bounded shear flows (Liu & Gayme 2021) demonstrated that including the componentwise structure of the nonlinearity uncovers a wider range of known key flow features identified through nonlinear analysis, experiments and DNS, but not captured through traditional (unstructured) input–output approaches. Here, SIOA for stratified PCF includes the effect of nonlinearity in the momentum and density equations (under the Boussinesq approximation) within a computationally tractable linear framework through a feedback interconnection between the linearized dynamics and a structured forcing that is explicitly constrained to preserve the componentwise structure of the nonlinearity. The structured singular value (Doyle 1982; Safonov 1982) of the spatio-temporal frequency response associated with this feedback interconnection can then be calculated at each streamwise and spanwise length scale. This structured singular value can be interpreted as the flow structure that shows the largest input–output gain (amplification) given the structured feedback interconnection.

Here, we apply the SIOA to characterize highly amplified flow structures in the intermittent regime of stratified PCF and investigate the behaviour of the flow across a range of $Re$, $Ri_b$ and $Pr$. Our aims are twofold. First, we wish to investigate whether the structures predicted by the SIOA can be quantitatively identified with fully nonlinear structures that have been observed in previously reported DNS of stratified PCF with specific values of the control parameters $Re$, $Ri_b$ and $Pr$. Second, we wish to explore the dependence on the control parameters of predictions from the SIOA in parameter regimes which are not (as yet) accessible to DNS. More specifically, to address our first aim, we examine how $Re$ and $Ri_b$ affect flow structures with Prandtl number set at $Pr=0.7$, i.e. the value for air. We demonstrate that SIOA does indeed predict the characteristic wavelengths and angle of the oblique turbulent bands observed in very large channel size DNS of the low-$Re$ low-$Ri_b$ intermittent regime of stratified PCF at the same $Pr$ (Deusebio et al. 2015; Taylor et al. 2016). We further show that, in the high-$Re$ high-$Ri_b$ intermittent regime, the SIOA identifies quasi-horizontal flow structures resembling turbulent–laminar layers (Deusebio et al. 2015).

Having achieved our first aim, and demonstrated the usefulness of the SIOA for identifying realistic nonlinear flow structures, we then turn our attention to our second aim. We demonstrate that increasing the bulk Richardson number reduces the amplification of streamwise-varying flow structures. These results show that the classical Miles–Howard stability criterion ($Ri_b\leq 1/4$) appears (perhaps fortuitously) to be associated with a change in the most amplified flow structures, which is robust for a wide range of $Re$ and valid at $Pr\approx 1$.

We then examine flow behaviour at different $Ri_b$ and $Pr$. For flows with $Pr\ll 1$, a larger bulk Richardson number is required to reduce the amplification of streamwise-varying flow structures to the same level as streamwise-independent ones compared with $Pr\approx 1$. The largest amplification also is predicted to occur at the same value of the product $Pr Ri_b$ consistent with the observation of the averaged density profile only varying with the product $Pr Ri_b$ in the $Pr\ll 1$ regime (Langham et al. 2020). For flows with $Pr\gg 1$, the SIOA identifies another quasi-horizontal flow structure independent of $Ri_b$. By decomposing input–output pathways into separate components associated with velocity and density fluctuations, we show that these quasi-horizontal flow structures at $Pr\gg 1$ are primarily associated with fluctuations in the density field. We further highlight the importance of this density-associated flow structure at $Pr\gg 1$ by constructing an analytical scaling argument for the input–output amplification in terms of $Re$ and $Pr$ under the assumptions of unstratified flow (with $Ri_b=0$) and streamwise invariance. The above observations using SIOA distinguish two types of quasi-horizontal flow structures, one associated with the high-$Re$ high-$Ri_b$ regime and the other one associated with density perturbations that emerges in the high $Pr$ regime.

To achieve our twin aims, and to demonstrate the above summarized results, the remainder of this paper is organized as follows. Section 2 describes the flow configuration of stratified PCF and then develops the SIOA for this flow. Section 3 analyses the results obtained from SIOA focusing on the wall-parallel length scale of flow structures. In § 4, we develop analytical scaling arguments with respect to $Re$ and $Pr$ to investigate behaviour for flows in the high $Pr$ limit. Finally we draw conclusions and suggest some avenues of future work in § 5.

2. Structured input–output response of stratified flow

2.1. Governing equations

We consider stably stratified PCF between two infinite parallel plates and employ $x$, $y$ and $z$ to denote the streamwise, wall-normal (or vertical) and spanwise directions. The corresponding (assumed incompressible) velocity components are denoted as $u$, $v$ and $w$. The coordinate frames and configurations for this stratified PCF are shown in figure 1. We express the velocity field as a vector $\boldsymbol {u}_{{tot}}=[u_{{tot}} \enspace v_{{tot}} \enspace w_{{tot}}]^\text {T}$ with ${}^\text {T}$ indicating the transpose. We then decompose the velocity field into the sum of a laminar linearly varying base flow $U(y)=y$ and fluctuations about the base flow $\boldsymbol {u}$; i.e. $\boldsymbol {u}_{{tot}}=U(y)\boldsymbol {e}_x+\boldsymbol {u}$ with $\boldsymbol {e}_x$ denoting the $x$-direction (streamwise) unit vector. The pressure field is similarly decomposed as $p_{{tot}}=P+p$. We decompose the density $\rho _{{tot}}$ as the sum of a reference density $\rho _r$, a base, again linearly varying, density $\bar {\rho }=-y$ and a density fluctuation $\rho$; i.e. $\rho _{{tot}}=\rho _r+\bar {\rho }+\rho$. We use $\rho _0$ to denote half of the density difference between the top and bottom walls, which is assumed to be much smaller than the reference density $\rho _0\ll \rho _r$ so that the Boussinesq approximation can be used.

Figure 1. (a) Schematic of stably stratified PCF with laminar base flow $\boldsymbol {U}(y)=y\boldsymbol {e}_x$ and background density $\bar {\rho }=-y$. The gravity vector $\boldsymbol {g}=-g\boldsymbol {e}_y$ is orthogonal to the channel walls. The blue and red colours represent high-density and low-density fluids, respectively. (b) Block diagram of the feedback interconnection between the linearized dynamics and structured forcing (outlined by the blue dashed line) in (2.4) modelling the nonlinearity.

The dynamics of the fluctuations $\boldsymbol {u}$, $p$ and $\rho$ is hence governed by the Navier–Stokes equations for an incompressible velocity field under the Boussinesq approximation and an advection–diffusion equation for the density

(2.1a)\begin{gather} \partial_{t} \boldsymbol{u} + U\partial_x \boldsymbol{u} + v\frac{{\rm d} U }{{\rm d}y}\boldsymbol{e}_x+Ri_b\rho\boldsymbol{e}_y +\boldsymbol{\nabla} p -\frac{1}{Re}{\nabla}^2 \boldsymbol{u} ={-} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}, \end{gather}

(2.1b)\begin{gather}\partial_t \rho+U\partial_x \rho+v\frac{{\rm d}\bar{\rho}}{{\rm d}y}-\frac{1}{RePr}\nabla^2 \rho ={-}\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\rho, \end{gather}

(2.1c)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0. \end{gather}

Here, the spatial variables are normalized by the channel half-height $h$, and the velocity is normalized by half of the velocity difference between the top and bottom walls $U_w$, where $\pm U_w$ is the velocity at the channel walls. Time and pressure are normalized by $h/U_w$ and $\rho _r U_w^2$, respectively. The base density field $\bar {\rho }(y)$ and the density fluctuations $\rho$ are normalized by $\rho _0$. Under this normalization, the base density profile $\bar {\rho }=-y$ is balanced by a hydrostatic pressure $P=Ri_b y^2/2$.

The non-dimensional control parameters are the Reynolds number $Re$, the Prandtl number $Pr$ and the bulk Richardson number $Ri_b$, naturally defined as

(2.2a–c)\begin{equation} Re:=\frac{U_wh}{\nu},\quad Pr:=\frac{\nu}{\kappa},\quad Ri_b:=\frac{g \rho_0 h}{\rho_r U_w^2}, \end{equation}

where $\nu$ is the kinematic viscosity, $\kappa$ is the molecular diffusivity of the density scalar and $g$ is the magnitude of gravity. The gravity vector is in the direction orthogonal to the wall $\boldsymbol {g}=-g\boldsymbol {e}_y$ with $\boldsymbol {e}_y$ denoting the $y$-direction (wall-normal, or vertical) unit vector. In equation set (2.1), $\boldsymbol {\nabla }:=[\partial _x \enspace \partial _y \enspace \partial _z]^\text {T}$ represents the gradient operator, and $\nabla ^2:=\partial _{x}^2+\partial _{y}^2+\partial _{z}^2$ represents the Laplacian operator. We impose no-slip boundary conditions at the wall $\boldsymbol {u}(y=\pm 1)=\boldsymbol {0}$ and Dirichlet boundary conditions for density fluctuations $\rho (y=\pm 1)=0$ that can be maintained by e.g. constant temperatures at the wall with a linear equation of state (with the hotter plate at the top).

A large body of linear analysis techniques views the nonlinear terms as a forcing, which enables these terms to be represented as an unmodelled effect (which can be thought of as some type of ‘uncertainty’ in the equations). There are a wide range of such models, but a common approach is a delta-correlated or coloured stochastic forcing that captures a wide range of the unmodelled effects, see e.g. the discussion in Farrell & Ioannou (1993a), Bamieh & Dahleh (2001), Jovanović & Bamieh (2005), McKeon & Sharma (2010), McKeon (2017) and Zare, Jovanović & Georgiou (2017). Here, we similarly write the nonlinear terms as the forcing

(2.3a)\begin{gather} \boldsymbol{f}_{\boldsymbol{u}}:={-}\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}=\begin{bmatrix} -\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}u & -\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}v & -\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}w\end{bmatrix}^\text{T}=:\begin{bmatrix} f_x & f_y & f_z\end{bmatrix}^\text{T}. \end{gather}

(2.3b)\begin{gather}f_{\rho}:={-}\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\rho, \end{gather}

which turns (2.1) into a set of linear evolution equations subject to the forcing terms $\boldsymbol {f}_{\boldsymbol {u}}$ and $f_{\rho }$.

We now construct a model of the nonlinearity, where the velocity field $-\boldsymbol {u}$ in (2.3) associated with the forcing components can be viewed as the gain operator of an input–output system in which the velocity and density gradients $\boldsymbol {\nabla } u$, $\boldsymbol {\nabla } v$, $\boldsymbol {\nabla } w$ and $\boldsymbol {\nabla }\rho$ act as the respective inputs and the forcing components $f_x$, $f_y$, $f_z$ and $f_{\rho }$ act as the respective outputs. This input–output models of the nonlinear components in the momentum and density equations (2.3), are respectively given by

(2.4a)\begin{gather} \boldsymbol{f}_{\boldsymbol{u},\xi}:={-}\boldsymbol{u}_{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}=\begin{bmatrix} -\boldsymbol{u}_{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}u & -\boldsymbol{u}_{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}v & -\boldsymbol{u}_{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}w \end{bmatrix}^\text{T}=:\begin{bmatrix} f_{x,\xi} & f_{y,\xi} & f_{z,\xi}\end{bmatrix}^\text{T}, \end{gather}

(2.4b)\begin{gather}f_{\rho,\xi}:={-}\boldsymbol{u}_{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}\rho. \end{gather}

Here, $-\boldsymbol {u}_{\xi }$ in (2.4) maps the corresponding velocity and density gradients into each component of the modelled forcing driving the linearized dynamics. This forcing in (2.4) is referred to as structured forcing because it preserves the componentwise structure of the nonlinear terms in (2.3). Figure 1(b) shows a block diagram of the feedback interconnection between the linearized dynamics and this forcing whose block-diagonal structure mirrors the nonlinear interactions in the Navier–Stokes equations, i.e. the forcing does not include terms such as $-\boldsymbol {u}\boldsymbol {\cdot}\boldsymbol {\nabla } v$ in the forcing $f_{x,\xi }$ since this term does not appear in the Navier–Stokes equations.

Although the nonlinearity in (2.3) can be written in many different ways, the current formulation leads to a straightforward and unified formulation for structured forcing in each momentum and density equation in (2.4). We next exploit this form of the equations to construct an input–output map using the structured singular value formalism (Packard & Doyle 1993; Zhou, Doyle & Glover 1996). This map will enable us to analyse the fluctuations which are prominent in the intermittent regime.

2.2. Structured input–output response

We need to define the spatio-temporal frequency response $\mathcal {H}_\nabla ^S(y;k_x,k_z,\omega )$ of stratified PCF that will form the basis of the structured input–output response. We use the superscript $^S$ to distinguish this operator from its counterpart for unstratified wall-bounded shear flow (Liu & Gayme 2021). We employ the standard transformation to express the velocity field dynamics in (2.1) in terms of the wall-normal velocity $v$ and the wall-normal vorticity $\omega _y:=\partial _z u-\partial _x w$ (Schmid & Henningson 2012). This transformation enforces the incompressibility constraint in (2.1c) and eliminates the pressure by construction. We exploit shift invariance in the $(x,z)$ spatial directions and assume shift invariance in time $t$, which allows us to perform the following triple Fourier transform:

(2.5)\begin{equation} \hat{\psi}(y;k_x,k_z,\omega):=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\psi(x,y,z,t)\exp({-\text{i}(k_x x + k_z z+\omega t)})\,{{\rm d}\kern0.7pt x}\,{\rm d} z\,{\rm d} t,\end{equation}

where $\text {i}=\sqrt {-1}$, $\omega$ is the temporal frequency and $k_x = 2{\rm \pi} /\lambda _x$ and $k_z = 2{\rm \pi} /\lambda _z$ are the $x$ and $z$ wavenumbers, respectively. The sign of the temporal frequency $\omega$ in (2.5) is chosen for ease of employing control-oriented toolboxes in our computations.

The resulting equations describing the transformed linearized equations subject to the forcing $\left [\begin {smallmatrix}\boldsymbol {f}_{\boldsymbol {u}, \xi }\\ f_{\rho,\xi }\end {smallmatrix}\right ]$ are given by

(2.6a)\begin{gather} \text{i}\omega \begin{bmatrix} \hat{v}\\ \hat{\omega}_y\\ \hat{\rho} \end{bmatrix}=\hat{\mathcal{A}}^S\begin{bmatrix} \hat{v}\\ \hat{\omega}_y\\ \hat{\rho} \end{bmatrix}+\hat{\mathcal{B}}^S\begin{bmatrix}\hat{f}_{x,\xi}\\ \hat{f}_{y,\xi}\\ \hat{f}_{z,\xi}\\ \hat{f}_{\rho,\xi}\end{bmatrix}, \end{gather}

(2.6b)\begin{gather}\begin{bmatrix}\hat{u}\\ \hat{v}\\ \hat{w}\\ \hat{\rho}\end{bmatrix}=\hat{\mathcal{C}}^S\begin{bmatrix} \hat{v}\\ \hat{\omega}_y\\ \hat{\rho} \end{bmatrix}. \end{gather}

The operators in equation set (2.6) are defined as

(2.7a)\begin{gather} \hat{\mathcal{A}}^S(k_x,k_z):= \hat{\mathcal{M}}^{{-}1} \begin{bmatrix} -\text{i}k_xU{\widehat{\nabla}}^2+\text{i}k_xU''+\dfrac{\widehat{{\nabla}}^4}{Re} & 0 & Ri_b(k_x^2+k_z^2)\\ -\text{i}k_zU' & -\text{i}k_x U+\dfrac{\widehat{{\nabla}}^2}{Re} & 0\\ -\bar{\rho}' & 0 & -\text{i}k_x U+\dfrac{\widehat{\nabla}^2}{Re Pr} \end{bmatrix}, \end{gather}

(2.7b)\begin{gather}\mathcal{\hat{B}}^S(k_x,k_z):=\hat{\mathcal{M}}^{{-}1} \begin{bmatrix} -\text{i}k_x\partial_y & -(k_x^2+k_z^2) & -\text{i}k_z \partial_y & 0\\ \text{i}k_z & 0 & -\text{i}k_x & 0\\ 0 & 0 & 0 & \mathcal{I} \end{bmatrix},\quad \hat{\mathcal{M}}:=\begin{bmatrix} \widehat{{\nabla}}^2 & 0 & 0\\ 0 & \mathcal{I} & 0\\ 0 & 0 & \mathcal{I} \end{bmatrix}, \end{gather}

(2.7c)\begin{gather}\mathcal{\hat{C}}^S(k_x,k_z):= \frac{1}{k_x^2+k_z^2}\begin{bmatrix} \text{i}k_x\partial_y & -\text{i}k_z & 0 \\ k_x^2+k_z^2 & 0 & 0\\ \text{i}k_z \partial_y & \text{i}k_x & 0 \\ 0 & 0 & k_x^2+k_z^2 \end{bmatrix}, \end{gather}

where $U':={\rm d} U(y)/{{\rm d}y}$, $U'':={\rm d}^2U(y)/{{\rm d}y}^2$, $\bar {\rho }':={\rm d}\bar {\rho }(y)/{{\rm d}y}$, $\widehat {{\nabla }}^2:=\partial _{yy}-k_x^2-k_z^2$, $\widehat {\nabla }^4:=\partial ^{(4)}_{y}-2(k_x^2+k_z^2)\partial _{yy}+(k_x^2+k_z^2)^2$ and $\mathcal {I}$ is the identity operator. The equation associated with the $\hat {\mathcal {A}}^S$ operator in (2.7a) can also be obtained by generalizing the classical Taylor–Goldstein equation (Goldstein 1931; Taylor 1931; Smyth & Carpenter 2019) to include viscosity, density diffusivity and coupling with wall-normal vorticity $\hat {\omega }_y$ with $k_z\neq 0$. The boundary conditions associated with (2.6) are $\hat {v}(y=\pm 1)= ({\partial \hat {v}}/{\partial y})(y=\pm 1)=\hat {\omega }_y(y=\pm 1)=\hat {\rho }(y=\pm 1)=0$.

The spatio-temporal frequency response $\mathcal {H}^S$ of the system in (2.6), which maps the input forcing to the velocity and density fields at the same spatial-temporal wavenumber–frequency triplet; i.e. $\left [\begin {smallmatrix}\boldsymbol {\hat {u}}(y;k_x,k_z,\omega )\\ {\hat {\rho }}(y;k_x,k_z,\omega ) \end {smallmatrix}\right ]=\mathcal {H}^S(y;k_x,k_z,\omega )\left [\begin {smallmatrix} \boldsymbol {\hat {f}}_{\boldsymbol {u}, \xi }(y;k_x,k_z,\omega )\\ \hat {f}_{\rho,\xi }(y;k_x,k_z,\omega )\end {smallmatrix}\right ]$, is given by

(2.8)\begin{equation} \mathcal{H}^S(y;k_x,k_z,\omega):=\hat{\mathcal{C}}^S \left(\text{i}\omega\mathcal{I}_{3\times 3}-\hat{\mathcal{A}}^S\right)^{{-}1}\hat{\mathcal{B}}^S. \end{equation}

Here, $\mathcal {I}_{3\times 3}:=\text {diag}(\mathcal {I},\mathcal {I},\mathcal {I})$, where $\text {diag}(\,\cdot\,)$ indicates a block-diagonal operation.

The linear form of (2.4a)–(2.4b) allows us to perform the same spatio-temporal Fourier transform on the model of the nonlinearity, which can be decomposed as

(2.9)\begin{equation} \begin{bmatrix} \hat{f}_{x,\xi}\\ \hat{f}_{y,\xi}\\ \hat{f}_{z,\xi}\\ \hat{f}_{\rho,\xi} \end{bmatrix}=\text{diag} \left(-\hat{\boldsymbol{u}}_{\xi}^{\text{T}},-\hat{\boldsymbol{u}}_{\xi}^{\text{T}}, -\hat{\boldsymbol{u}}_{\xi}^{\text{T}}, -\hat{\boldsymbol{u}}_{\xi}^{\text{T}} \right) \text{diag}\left(\widehat{\boldsymbol{\nabla}},\widehat{\boldsymbol{\nabla}}, \widehat{\boldsymbol{\nabla}},\widehat{\boldsymbol{\nabla}}\right) \begin{bmatrix} \hat{u}\\ \hat{v}\\ \hat{w}\\ \hat{\rho} \end{bmatrix}. \end{equation}

A block diagram illustrating this decomposition of the modelled nonlinearity is shown inside the blue dashed line ($\text {- -}$, blue) in figure 2(a). This block-diagonal structure constrains the modelled nonlinear interactions, i.e. provides structured forcing.

Figure 2. Block diagram showing the feedback interconnection of the SIOA framework applied to stratified PCF. (a) Redraws figure 1(b), where the blocks inside of ($\text {- -}$, blue) lines represent the modelled forcing in (2.9) corresponding to the bottom block in figure 1(b) also inside of ($\text {- -}$, blue). (b) Redraws (a) after discretization with the top block corresponding to the combination of the two top blocks in (a) and the bottom block corresponding to the bottom block of (a).

In order to isolate the gain operator $-\boldsymbol {u}_{\xi }$, we combine the linear gradient operator with the spatio-temporal frequency response of the linearized system (2.8). The resulting modified frequency response operator with outputs that are the vectorized gradients of the velocity and density components is defined as

(2.10)\begin{equation} \mathcal{H}^S_{\boldsymbol{\nabla}}(y;k_x,k_z,\omega):=\text{diag} \left(\widehat{\boldsymbol{\nabla}},\widehat{\boldsymbol{\nabla}}, \widehat{\boldsymbol{\nabla}},\widehat{\boldsymbol{\nabla}}\right) \mathcal{H}^S(y;k_x,k_z, \omega). \end{equation}

The resulting system model can be redrawn as a feedback interconnection between this linear operator and the structured uncertainty

(2.11)\begin{equation} \hat{\boldsymbol{u}}^S_{\varXi}:=\text{diag}\left(-\hat{\boldsymbol{u}}_{\xi}^{\text{T}}, -\hat{\boldsymbol{u}}_{\xi}^{\text{T}} ,-\hat{\boldsymbol{u}}_{\xi}^{\text{T}}, -\hat{\boldsymbol{u}}_{\xi}^{\text{T}} \right). \end{equation}

Here, the structure is introduced in terms of the diagonal form of $\hat {\boldsymbol {u}}_{\varXi }^S$ that enforces the componentwise structure of the nonlinearity in the forcing model defined in (2.4). Figure 2(b) illustrates this feedback interconnection between the modified spatio-temporal frequency response and the structured uncertainty, where $\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla }}$ and $\hat {\boldsymbol{\mathsf{u}}}^S_{\varXi }$, respectively, represent the spatial discretizations (numerical approximations) of $\mathcal {H}^S_{\boldsymbol {\nabla }}$ in (2.10) and $\boldsymbol {\hat {u}}^S_{\varXi }$ in (2.11).

We are interested in characterizing the horizontal length scales of the most amplified flow structures under this structured forcing. This amplification can be quantified in terms of the structured singular value of the modified frequency response operator $\mathcal {H}^S_{\boldsymbol {\nabla }}$; see e.g. Packard & Doyle (1993, definition 3.1) and Zhou et al. (1996, definition 11.1), which is defined as follows.

Definition 2.1 Given wavenumber and frequency pair $(k_x,k_z,\omega )$, the structured singular value $\mu _{\hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }}[\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla }}(k_x,k_z,\omega )]$ is defined as

(2.12)\begin{equation} \mu_{\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi}} \left[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega)\right]:= \frac{1}{\min\{\bar{\sigma}[\hat{\boldsymbol{\mathsf{u}}}^S_{\varXi}]: \hat{\boldsymbol{\mathsf{u}}}^S_{\varXi}\in \hat{\boldsymbol{\mathsf{U}}}^S_{\varXi}, \det[\boldsymbol{\mathsf{I}}-\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega) \hat{\boldsymbol{\mathsf{u}}}^S_{\varXi}]=0\}}.\end{equation}

If no $\hat {\boldsymbol{\mathsf{u}}}^S_{\varXi }\in \hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }$ makes $\boldsymbol{\mathsf{I}}-\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla }}\hat {\boldsymbol{\mathsf{u}}}^S_{\varXi }$ singular, then $\mu _{\hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }}[\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla }}]:=0$.

Here, $\bar {\sigma }[\,\cdot\,]$ is the largest singular value, $\det [\,\cdot\,]$ is the determinant of the argument, and $\boldsymbol{\mathsf{I}}$ is the identity matrix. The subscript of $\mu$ in (2.12) is a set $\hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }$ containing all uncertainties having the same block-diagonal structure as $\hat {\boldsymbol{\mathsf{u}}}^S_{\varXi }$; i.e.

(2.13)\begin{equation} \hat{\boldsymbol{\mathsf{U}}}^S_{\varXi}:=\left\{{\rm {diag}} \left(-\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}},-\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}}, -\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}},-\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}}\right): -\hat{\boldsymbol{\mathsf{u}}}^{\text{T}}_{\xi}\in \mathbb{C}^{N_y\times 3N_y}\right\}, \end{equation}

where $N_y$ denotes the number of grid points in $y$.

Definition 2.1 suggests that the inverse of the structured singular value $1/\mu$ is the minimal norm of the perturbation $\hat {\boldsymbol{\mathsf{u}}}^S_{\varXi }$ that destabilizes the feedback interconnection in figure 2(b) in the input–output sense defined by the small gain theorem, see Liu & Gayme (2021, proposition 2.2) and Zhou et al. (1996, theorem 11.8). This interpretation suggests that the flow field is more sensitive to perturbations with the flow structures associated with a larger amplification measured by the structured singular value $\mu$. A similar notion of destabilizing perturbation was also employed to interpret the largest (unstructured) singular value (Trefethen et al. 1993, p. 581).

Here, the form of structured uncertainty in (2.13) allows the full degrees of freedom for the complex matrix $-\hat {\boldsymbol{\mathsf{u}}}^{\text {T}}_{\xi }\in \mathbb {C}^{N_y\times 3N_y}$ for ease of computation. While $\boldsymbol {u}_\xi$ is not constrained to be incompressible, the incompressibility of $\boldsymbol {u}$ and the role of pressure are accounted for within the current $v$-$\omega _y$ formulation. Further refinement to better represent the physics and uncover the form of $\boldsymbol {u}_\xi$ requires an extension of both the analysis method and computational tools. These extensions are beyond the scope of the current work.

We then define the structured response following Liu & Gayme (2021) as

(2.14)\begin{equation} \|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\mu}(k_x,k_z):=\underset{\omega \in \mathbb{R}}{\text{sup}} \mu_{\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi}}\left[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega)\right], \end{equation}

where $\text {sup}$ represents a supremum (least upper bound) operation. Here, we abuse the notation and terminology by writing $\|\,\cdot\,\|_{\mu }$ (Packard & Doyle 1993), although this quantity is not a proper norm. We employ this notation in analogy with the corresponding unstructured response of the feedback interconnection, which is given by

(2.15)\begin{equation} \|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\infty}(k_x,k_z):=\underset{\omega \in \mathbb{R}}{\text{sup}} \bar{\sigma}\left[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega)\right].\end{equation}

This quantity is the unstructured counterpart of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$, which is obtained by replacing the structured uncertainty set $\hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }$ with the set of full complex matrices $\mathbb {C}^{4N_y\times 12N_y}$. In both cases, a larger value indicates that the corresponding flow structures (associated with a particular $k_x$ and $k_z$ pair) have larger amplification under either structured or unstructured feedback forcing. For example, a larger value of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }(k_x,k_z)$ indicates that the corresponding flow structures (associated with a particular $k_x$ and $k_z$ pair) have larger amplification under structured feedback in figure 2(b).

2.3. Numerical method

We employ the Chebyshev differential matrix (Trefethen 2000; Weideman & Reddy 2000) to discretize the operators in equation set (2.7). Our code is validated through comparison with the unstratified PCF and Poiseuille flow results in Jovanović (2004), Jovanović & Bamieh (2005) and Schmid (2007). The implementation of stratification is validated by reproducing the maximum growth rate of the linear normal mode in a layered stratified PCF determined by Eaves & Caulfield (2017, figures 3 and 6a), as well as the linear stability predictions for the unstable stratification configuration in Olvera & Kerswell (2017, figure 1 and appendix B). We use $N_y=60$ collocation points not including the boundary points over the wall-normal extent, as well as $48$ and $36$ logarithmically spaced streamwise and spanwise wavenumbers in the respective spectral ranges $k_x \in [10^{-4},10^{0.48}]$ and $k_z \in [10^{-2},10^{1.2}]$, unless otherwise mentioned. To verify that this resolution is sufficient to achieve grid convergence we recomputed selected results with 1.5 times the number of collocation points in the wall-normal direction and verified that the results did not change. The quantity $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ in (2.14) for each wavenumber pair $(k_x,k_z)$ is computed using the mussv command in the Robust Control Toolbox (Balas et al. 2005) of MATLAB. The arguments of mussv employed here include the state-space model of $\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla }}$ that samples the frequency domain adaptively. The BlockStructure argument comprises four full $N_y\times 3N_y$ complex matrices, and we use the ‘Uf’ algorithm option.

3. Structured spatio-temporal frequency response of stratified flow

In this section, we use the SIOA approach described in § 2.2 to characterize the flow structures that are most amplified in stably stratified PCF.

3.1. Low-$Re$ low-$Ri_b$ vs high-$Re$ high-$Ri_b$ intermittency

In this subsection, we analyse flow structures that are prominent in either the low-$Re$ low-$Ri_b$ or the high-$Re$ high-$Ri_b$ intermittent regime (Deusebio et al. 2015). Here, we keep the Prandtl number fixed at $Pr=0.7$. This value corresponds to thermally stratified air and is the same value studied by Deusebio et al. (2015). We first consider a flow with $Re=865$, $Ri_b=0.02$ and $Pr=0.7$, where oblique turbulent bands have been observed (Deusebio et al. 2015; Taylor et al. 2016). In order to evaluate the relative effect of the feedback interconnection and the imposed structure, we also compute $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }(k_x,k_z)$ defined in (2.15) and

(3.1)\begin{equation} \|\mathcal{H}^S\|_{\infty}(k_x,k_z):= \underset{\omega \in \mathbb{R}}{\text{sup}} \bar{\sigma}\left[\boldsymbol{\mathsf{H}}^S(k_x,k_z,\omega)\right].\end{equation}

Here, $\boldsymbol{\mathsf{H}}^S$ is the discretization of spatio-temporal frequency response operator $\mathcal {H}^S$ in (2.8), i.e. the spatio-temporal frequency operator governing the linearized dynamics without the feedback interconnection. The values of $\|\mathcal {H}\|_{\infty }$ for unstratified plane Couette and plane Poiseuille flows were previously analysed in Jovanović (2004), Schmid (2007) and Illingworth (2020). The quantity in (3.1) describes the maximum singular value of the frequency response operator $\mathcal {H}^S$, which represents the maximal gain of $\mathcal {H}^S$ over all temporal frequencies; i.e. the worst-case amplification over harmonic inputs.

Figure 3 shows $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ in (a) alongside (b) $\|\mathcal {H}^S\|_{\infty }$ and (c) $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }$. We indicate the characteristic wavelength pair $\lambda _x=32{\rm \pi}$, $\lambda _z=16{\rm \pi}$ corresponding to the oblique turbulent bands observed in DNS under the same flow regime (Deusebio et al. 2015; Taylor et al. 2016, figure 2b) in these panels using the symbol ($\triangle$, black). These structures are observed to have a characteristic inclination angle (measured from the streamwise direction in $x$–$z$ plane) of $\theta :=\tan ^{-1}(\lambda _z/\lambda _x)\approx 27^\circ$, which is indicated in all panels by the black solid line (–) that plots $\lambda _z=\lambda _x\tan (27^\circ )$. While there is some footprint of these structures and this angle in all three panels, the correspondence with the peak amplitude is most prominent in (a). In fact, the peak value of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ in (a) occurs at streamwise and spanwise wavenumbers associated with the characteristic wavelengths and angle of the oblique turbulent bands reported in DNS (Deusebio et al. 2015; Taylor et al. 2016), and the line representing the angle of the oblique turbulent bands crosses through the centre of the narrow roughly elliptical peak region whose principal axis coincides with this angle. The results in figure 3(a) suggest that the SIOA captures both the wavelengths and angle of the oblique turbulent bands in the low-$Re$ low-$Ri_b$ intermittent regime of stratified PCF. This analysis suggests that these oblique turbulent bands arise in the intermittent regime of stratified PCF due to their large amplification, or equivalently their sensitivity to disturbances.

Figure 3. Plots of (a) $\log _{10}[\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }(k_x,k_z)]$, (b) $\log _{10}[\|\mathcal {H}^S\|_{\infty }(k_x,k_z)]$ and (c) $\log _{10}[\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }(k_x,k_z)]$ for stratified PCF at $Re=865$, $Ri_b=0.02$ and $Pr=0.7$. Here, the symbols (${\triangle }$) are characteristic wavelengths ($\lambda _x=32{\rm \pi}$, $\lambda _z=16{\rm \pi}$) corresponding to the oblique turbulent band observed in DNS in the same flow regime (Deusebio et al. 2015; Taylor et al. 2016). The lines (–) are $\lambda _z=\lambda _x\tan (27^\circ )$ indicating the $27^\circ$ angle of the oblique turbulent bands.

The traditional input–output analysis results, $\|\mathcal {H}^S\|_\infty$ in (b), provide a noticeable improvement compared with growth rate analysis (as presented in more detail in Appendix A) and are also able to identify the preferred wavenumber pair in this intermittent regime. However, this analysis suggests larger amplification of the streamwise elongated modes. Moreover, the inclusion of an unstructured feedback loop quantified through $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }$ in (c) correctly orders the relative amplification between the oblique turbulent bands and streamwise elongated structures ($k_x\approx 0$). The differences between $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }$ and $\|\mathcal {H}^S\|_{\infty }$ are likely associated with the additional $\widehat {\boldsymbol {\nabla }}$ operator in defining $\mathcal {H}_{\boldsymbol {\nabla }}$ in (2.10), which emphasizes flow structures with a larger horizontal wavenumber. The difference between $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ and $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }$ is associated with the structured feedback interconnection that constrains the permissible feedback pathway, which weakens the amplification associated with the lift-up mechanism; see similar discussion on unstratified PCF (Liu & Gayme 2021, § 3.3). A comparison of the results in figure 3 suggests that it is the imposition of the componentwise structure from the nonlinear terms in (2.3) that further improves the prediction of the oblique turbulent bands.

We now consider the high-$Re$ high-$Ri_b$ intermittent regime, which was shown to be qualitatively different in behaviour from the low-$Re$ low-$Ri_b$ intermittent regime (Deusebio et al. 2015). We first isolate the effect of increasing either $Re$ or $Ri_b$. Figure 4(a) presents $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ for a flow with $Ri_b=0.02$ and $Re=4250$. The larger colour bar range vs figure 3 highlights the expected higher magnitudes vs those for a flow with a lower Reynolds number ($Re=865$). We can see that the wavenumber pair of the peak region extends towards smaller values (larger wavelengths) than those associated with the oblique turbulent bands that were in the peak region in figure 3(a). Figure 4(b) presents $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ for a higher bulk Richardson number $Ri_b=0.2541$ and the same $Re$ and $Pr$ values as figure 3(a). Here, the amplification associated with the streamwise-varying flow structures such as the oblique turbulent bands observed in figure 3(a) is reduced and quasi-horizontal flow structures $(k_x\approx 0, k_z\approx 0)$ show a similar level of amplification (see the bottom left corner in figure 4b). Note that this flow structure associated with $k_x\approx 0$, $k_z\approx 0$ is referred to as quasi-horizontal to distinguish it from a horizontally uniform mode ($k_x=0$, $k_z=0$).

Figure 4. Plots of $\log _{10}[\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }(k_x,k_z)]$ at $Pr=0.7$ and: (a) $Re=4250$ and $Ri_b=0.02$; (b) $Re=865$ and $Ri_b=0.2541$ and (c) $Re=52\,630$ and $Ri_b=0.15$.

Armed with these insights, we consider the combined high-$Re$ high-$Ri_b$ intermittent regime ($Re=52\,630$ and $Ri_b=0.15$); these values correspond to the results shown in figure 7 of Deusebio et al. (2015). Figure 4(c) presents $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ for these parameter values with an increased wall-normal grid with $N_y=90$. Here, the amplification of the oblique turbulent band is of a similar order to that of flow structures with a wide range of wavenumber pairs ranging from $k_x\lesssim 10^{-2}$ and $k_z\lesssim 1$ down to $k_x\approx 0$ and $k_z\approx 0$. These latter flow structures resemble the quasi-horizontal flow structures that have a horizontal length scale much larger than their vertical length scales, which are limited by the channel height and therefore restricted to non-dimensional scales of the order of $2$. The response in this regime, therefore, shows a large qualitative difference from that in the low-$Re$ low-$Ri_b$ ($Re=865$ and $Ri_b=0.02$) intermittent regime shown in figure 3(a). This qualitative difference mirrors the different features in intermittent regimes described by Deusebio et al. (2015), where oblique turbulent bands are prevalent in the low-$Re$ low-$Ri_b$ intermittent regime, but the high-$Re$ high-$Ri_b$ intermittent regime is characterized by turbulent–laminar layers indicating a large horizontal length scale.

In figure 4(c), we also observe that the quasi-horizontal flow structures have a streamwise wavelength much larger than their spanwise wavelength ($\lambda _x\gg \lambda _z$), which is also consistent with the observation in Deusebio et al. (2015) that the turbulent and laminar layers in the high-$Re$ high-$Ri_b$ intermittent regime are homogeneous in the streamwise direction. This behaviour can be understood through an order of magnitude estimation of the terms in the continuity equation. We assume highly anisotropic flow with $v\approx 0$ under strong stratification, which simplifies the continuity equation to ${\partial u}/{\partial x}+{\partial w}/{\partial z}=0$. We further assume that the restoring buoyancy force due to stratification does not have a preference between the streamwise or spanwise direction and, therefore, we also assume ${\partial u}/{\partial x}$ and ${\partial w}/{\partial z}$ are the same order of magnitude. However, in the current stratified PCF configuration, streamwise velocity is associated with a characteristic velocity much larger than its spanwise counterpart due to the base flow velocity. As a result, streamwise variation is reduced much faster than spanwise variation (i.e. $k_x\ll k_z$) to keep ${\partial u}/{\partial x}$ and ${\partial w}/{\partial z}$ the same order of magnitude.

3.2. Case $Ri_b>1/4$: a change in the most amplified flow structures

The Miles–Howard theorem (Howard 1961; Miles 1961) implies that the laminar base flow would be linearly stable in the limits where $\nu$ and $\kappa$ both are zero if $Ri_b > 1/4$. Although the theorem is not applicable to unsteady flows with finite $Re$ and $Pr$, it has also been observed that a ‘marginal’ or ‘critical’ Richardson number near this threshold value appears to emerge naturally in simulations (Salehipour et al. 2018) and field measurement (Smyth et al. 2019). In the previous section, we noted that increasing $Ri_b=0.02$ to $Ri_b=0.2541$, for a fixed Reynolds number of $Re=865$, reduces the overall response and changes the types of flow structures, $(k_x, k_z)$ wavenumber pairs, that exhibit the largest response, see figures 3(a) and 4(b). In this subsection, we further investigate whether this apparently marginal threshold $Ri_b \simeq 1/4$ is associated with a change in flow structures and whether this behaviour is independent of $Re$. As in the previous subsection, the Prandtl number is fixed at $Pr=0.7$.

Here, we aggregate results varying over a range of $(k_x,k_z)$ wavenumber pairs in terms of the maximum value

(3.2)\begin{equation} \|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\mu}^M:=\underset{k_z,k_x}{\text{max}} \|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\mu}(k_x,k_z),\end{equation}

over the wavenumber domain $k_x\in [10^{-6},10^{0.48}]$ and $k_z\in [10^{-6}, 10^{0.48}]$. Lowering the minimum value of the wavenumber ranges vs those considered in the previous subsection is motivated by the observation in figure 4 that both the $k_x$ and $k_z$ values corresponding to the peak region decrease with increasing Reynolds and Richardson numbers. To separate streamwise-varying and streamwise-independent flow structures, we similarly evaluate

(3.3)\begin{equation} \|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\mu}^{sc}:= \underset{k_z,k_x=10^{{-}6}}{\text{max}}\|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\mu}(k_x,k_z). \end{equation}

This quantity restricts the streamwise wavenumber to $k_x=10^{-6}$ to approximate the streamwise constant modes and computes the maximum value over $k_z\in [10^{-6}, 10^{0.48}]$, where we use an upper bound of $10^{0.48}$ rather than the larger value of $10^{1.2}$ to save computation time. This change in the upper bound was not found to affect the results since the $k_z$ associated with the maximum value was consistently found to be below this upper bound. The restriction in streamwise wavelengths to $k_x=10^{-6}$ naturally includes the quasi-horizontal flow structures prevalent in the high-$Re$ high-$Ri_b$ regime ($k_x\approx 0$, $k_z\approx 0$) as an extreme case, but excludes streamwise-varying flow structures such as the oblique turbulent bands observed in the low-$Re$ low-$Ri_b$ regime discussed in § 3.1.

Figure 5 shows the variation of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^M$ (solid lines) and $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}$ (dashed lines) with bulk Richardson number $Ri_b\in [0,6]$ and Reynolds number $Re\in [865,15\,000]$. The quantities $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{M}$ including streamwise-varying flow structures are very close to $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}$ when $Ri_b \gtrsim 1/4$ for the full range of Reynolds numbers $Re\in [865,15\,000]$ in figure 5(a). This phenomenon is also reflected in figure 5(b), where for flows with $Ri_b=0.24$ and $Ri_b=0.75$, the curves for $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{M}$ and $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}$ largely overlap. These trends suggest that the inviscid marginal stability value $Ri_b=1/4$ predicted by the Miles–Howard theorem (Howard 1961; Miles 1961) for laminar flow is apparently associated with a change in the most amplified flow structure in stratified PCF at finite $Re$ and $Pr \sim 1$.

Figure 5. The dependence on $Ri_b$ and $Re$ of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ (–, black) and $\|\mathcal {H}^{S}_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ ($\text {- -}$, red) at $Pr=0.7$. Each black marker on the lines of $\|\mathcal {H}^{S}_{\boldsymbol {\nabla }}\|_{\mu }^{M}$ indicates the associated value of (a) $Re$ and (b) $Ri_b$ defined in the legend. Each red marker on the lines of $\|\mathcal {H}^{S}_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ indicates the same parameter value for (a) $Re$ and (b) $Ri_b$ as the corresponding black one.

The plots in figure 5 show that the largest amplification of the streamwise-invariant modes represented by $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}$ (dashed lines) do not appear to be influenced by $Ri_b$, as shown by the horizontal dashed lines in figure 5(a) and the overlapping dashed lines in figure 5(b). We further explore this $Ri_b$ independence for streamwise constant flow structures by considering the limit of horizontal invariance $\partial _x(\,\cdot\,)=0$ and $\partial _z(\,\cdot\,)=0$ ($k_x=0$ and $k_z=0$), which directly results in $\partial _y v=0$ due to the continuity equation. The boundary condition $v(y=\pm 1)=0$ then directly results in $v=0$. Using these assumptions, the advection terms vanish (i.e. $U\partial _x(\,\cdot\,)=0$, and $\boldsymbol {u}\boldsymbol {\cdot} \boldsymbol {\nabla }(\,\cdot\,)=0$) in each momentum and density equation. The terms associated with the background shear and density gradient also vanish due to zero wall-normal velocity; i.e. $vU'=v\bar {\rho }'=0$.

These observations lead to a simplification of the momentum and density equations in (2.1) to

(3.4a,b)\begin{gather} \partial_t u=\frac{1}{Re}\partial_{yy} u,\quad \partial_y p={-}Ri_b \rho, \end{gather}

(3.4c,d)\begin{gather}\partial_t w=\frac{1}{Re}\partial_{yy} w,\quad \partial_t \rho=\frac{1}{RePr}\partial_{yy} \rho. \end{gather}

Here, we can see that horizontal momentum and density field equations are all reduced to the diffusion equation, and the wall-normal momentum equation is reduced to a balance between the buoyancy force and the vertical pressure gradient i.e. to a hydrostatic balance. This balance suggests that the only dependence on $Ri_b$ can be absorbed into the pressure by rescaling pressure and thus does not influence the quasi-horizontal flow structures. The results presented in figures 3(a) and 4 suggest that flow structures with $k_x=10^{-4}$ lead to the same structured response $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_\mu$ at a wide range of spanwise wavenumbers $k_z\leq 1$, and this value is consistent with that of quasi-horizontal flow structures $(k_x\approx 0,k_z\approx 0)$ that are independent of $Ri_b$.

This analysis provides evidence that quasi-horizontal flow structures are associated with amplification which is independent of $Ri_b$. Instead, high $Ri_b$ (i.e. strong stratification) will reduce the amplification of other horizontally varying flow structures such as the oblique turbulent bands observed in the low-$Re$ low-$Ri_b$ intermittent regime of stably stratified PCF. Furthermore, it appears in figure 5(b) that quasi-horizontal flow structures are also increasingly amplified as $Re$ increases with scaling law $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}\sim Re$. This behaviour indicates that quasi-horizontal flow structures may prefer a high-$Re$ high-$Ri_b$ regime. In § 4, we further explore this scaling law $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}\sim Re$ by developing analytical scaling arguments for $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ in unstratified and streamwise-invariant flow.

3.3. Effects of low and high $Pr$

The Prandtl number is known to play an important role in the types of flow structures characterizing stratified PCF (Taylor & Zhou 2017; Zhou et al. 2017a,b; Langham et al. 2020). The Prandtl number also varies over a wide range in different applications. For example, $Pr\ll 1$ is relevant for flow in the stellar interior; see e.g. Garaud (2021), while a Prandtl number of $Pr=7$ corresponds to thermally stratified water. The Schmidt number (the analogous parameter to the Prandtl number for compositionally induced density variations) for salt-stratified water is significantly larger $Sc \simeq 700$. Moreover, the Prandtl number is obviously not well defined under the inviscid and non-diffusive assumptions of the Miles–Howard theorem. In this subsection, we explore the effect of low or high Prandtl number on flow structures. Here, we keep the Reynolds number fixed at $Re=4250$ following Zhou et al. (2017a,b). In order to resolve fully the additional scales introduced at high $Pr$, we increase the number of wall-normal grid points to $N_y=120$ at $Pr=70$, which is chosen as a more computationally accessible ‘large’ value, as previously considered by Zhou et al. (2017a,b).

We first investigate the effect of low $Pr$. The cross-channel density profiles of exact coherent structures in stratified PCF were shown to match in flows with the same $PrRi_b$ at $Pr\in [10^{-4},10^{-2}]$ (Langham et al. 2020, figure 3). This combined measure $PrRi_b$ has been proposed as the natural control parameter for stably stratified shear flows at the low Prandtl number limit $Pr\ll 1$ (Lignieres 1999; Garaud et al. 2015). In order to further explore this dependence, we plot $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^M$ as a function of $PrRi_b$ for Prandtl numbers in the range $Pr\in [10^{-4}, 7]$ in figure 6. Here, the results $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^M$ for $Pr\in [10^{-4},10^{-1}]$ (magenta dashed lines) again show a natural matching dependence on $PrRi_b$. This behaviour breaks down for flows with $Pr\geq 1$, as shown in figure 6. A similar end to the region of matched dependence on $PrRi_b$ alone is observed in studies using exact coherent structures, where the density profile at $Pr= 0.1$ deviates from the matching profiles for flows with $Pr\in [10^{-4}, 10^{-2}]$ yet fixed $Pr Ri_b$ (Langham et al. 2020, figure 3).

Figure 6. The dependence on $PrRi_b$ of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ at $Re=4250$.

In figure 7, we plot $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ and $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ as a function of $Ri_b$ and $Pr$ over the respective ranges $Ri_b\in [0,6]$ and $Pr\in [10^{-4},70]$. Figure 7(a) shows that the marginal stability value $Ri_b=1/4$ is not associated with any significant changes in the types of flow features undergoing the largest amplification for flows with $Pr=0.01$ ($\triangle$). Figure 7(b) further suggests that flows with smaller $Pr$ require a larger $Ri_b$ to reduce the amplification of streamwise-varying flow structures to the same level as streamwise constant structures. This behaviour is consistent with the observation that the exact coherent structures in the low $Pr$ limit require a larger $Ri_b$ to be affected by stratification in PCF (Langham et al. 2020). Figure 7 further shows that, for flows with high $Pr$, the quantities $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ and $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ are the same over a wide range of $Ri_b\in [0,6]$. In particular, the horizontal lines plotted in figure 7(a) for flows with different $Ri_b$ collapse to one line in the high $Pr\gg 1$ regime shown in figure 7(b). This observation is also consistent with Langham et al. (2020), who noted that, in the high Prandtl number limit $Pr\gg 1$, the effect of increasing $Ri_b$ is mitigated.

Figure 7. The dependence on $Ri_b$ and $Pr$ of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ (–, black) and $\|\mathcal {H}^{S}_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ ($\text {- -}$, red) at $Re=4250$. Each black marker on the lines of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ indicates the associated value of (a) $Pr$ and (b) $Ri_b$ indicated in the legend. Each red marker on the lines of $\|\mathcal {H}^{S}_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ indicates the same parameter as the corresponding black marker.

In order to investigate in isolation the effect of Prandtl number on the amplification of each wavenumber pair $(k_x,k_z)$, in figure 8 we plot $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }(k_x,k_z)$ for flows with (a) $Pr=10^{-4}$, (b) $Pr=7$ and (c) $Pr=70$. The peak region in figure 8(a) at $Pr=10^{-4}$ resembles the shape of the peak region for unstratified PCF (Liu & Gayme 2021, figure 4a). We have also computed the results for unstratified PCF at the same Reynolds number $Re=4250$ (not shown here), and we find almost the same results as shown in figure 8(a). This similarity suggests that, for the same bulk Richardson number, a lower Prandtl number will result in a weakening of the stabilizing effect of stratification. Comparing $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ in figure 4(a) at $Pr=0.7$ with the same quantity at $Pr=7$ and $Pr=70$, respectively shown in figures 8(b) and 8(c), we notice that the amplification associated with the wavenumber pair $k_x=10^{-4}$ and $k_z=10^{-2}$ increases with $Pr$. More specifically, the value of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ associated with $k_x=10^{-4}$ and $k_z=10^{-2}$ becomes comparable to the values associated with the $k_x\approx 10^{-2}$ and $k_z\approx 10^{-1}$ at $Pr=7$ as shown in 8(b). The wavenumber pair $k_x=10^{-4}$ and $k_z=10^{-2}$ is associated with the largest magnitude over the $(k_x,k_z)$ contour region at $Pr=70$, as shown in figure 8(c).

Figure 8. Plots of $\log _{10}[\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }(k_x,k_z)]$ at $Re=4250$, $Ri_b=0.02$ with three different Prandtl numbers at (a) $Pr=10^{-4}$, (b) $Pr=7$ and (c) $Pr=70$.

The quasi-horizontal flow structures ($k_x\approx 0$, $k_z\approx 0$) observed in flows at high $Pr$ have different features from those previously observed in the high-$Re$, high-$Ri_b$ regime (e.g. results for flows with $Re=52\,630$ and $Ri_b=0.15$ shown in figure 4c) and described in § 3.1. This indicates that a new type of quasi-horizontal flow structure appears in flows with sufficiently high $Pr$. The appearance of this flow structure at a high $Pr$ suggests that this flow structure is associated with fluctuations in the density field. This can be further explored by isolating the input–output pathway for each component of the momentum and density equations, i.e. inputs $f_x, f_y,\ f_z,\ f_{\rho }$ in (2.3) to outputs $u$, $v$, $w$ and $\rho$. These input–output pathways can be studied through the definition of operators $\mathcal {H}^S_{ij}$, where $j$ defines the forcing input component ($j=x,y,z,\rho$) and $i=u,v,w,\rho$ describes each velocity or density output component

(3.5)\begin{equation} \mathcal{H}^S_{ij}=\hat{\mathcal{C}}^S_{i}\left(\text{i}\omega \mathcal{I}_{3\times 3}-\hat{\mathcal{A}}^S\right)^{{-}1}\hat{\mathcal{B}}^S_{j},\end{equation}

with

(3.6a,b)\begin{gather} \mathcal{\hat{B}}^S_x:=\mathcal{\hat{B}}^S\begin{bmatrix} \mathcal{I} & 0 & 0 & 0 \end{bmatrix}^{\text{T}},\quad \mathcal{\hat{B}}^S_y:=\mathcal{\hat{B}}^S\begin{bmatrix} 0 & \mathcal{I} & 0 & 0 \end{bmatrix}^{\text{T}}, \end{gather}

(3.6c,d)\begin{gather}\mathcal{\hat{B}}^S_z:=\mathcal{\hat{B}}^S\begin{bmatrix} 0 & 0 & \mathcal{I} & 0 \end{bmatrix}^{\text{T}},\quad \mathcal{\hat{B}}^S_\rho:=\mathcal{\hat{B}}^S\begin{bmatrix} 0 & 0 & 0 & \mathcal{I} \end{bmatrix}^{\text{T}}, \end{gather}

(3.6e,f)\begin{gather}\mathcal{\hat{C}}^S_u:=\begin{bmatrix}\mathcal{I} & 0 & 0 & 0 \end{bmatrix} \mathcal{\hat{C}}^S ,\quad \mathcal{\hat{C}}^S_v:=\begin{bmatrix}0 & \mathcal{I} & 0 & 0 \end{bmatrix} \mathcal{\hat{C}}^S, \end{gather}

(3.6g,h)\begin{gather}\mathcal{\hat{C}}^S_w:=\begin{bmatrix}0 & 0 & \mathcal{I} & 0 \end{bmatrix} \mathcal{\hat{C}}^S,\quad \mathcal{\hat{C}}^S_\rho:=\begin{bmatrix}0 & 0 & 0 & \mathcal{I} \end{bmatrix} \mathcal{\hat{C}}^S. \end{gather}

Figures 9 and 10 present quantities $\|\mathcal {H}^S_{ij}\|_{\infty }$ $ij=ux,vy,wz$ and $\rho \rho$ for the control parameters ($Re=865, Ri_b=0.02, Pr=0.7)$ and $(Re=4250,Ri_b=0.02, Pr=70)$ respectively. The combined effect of these four panels associated with the input and output in the same component resembles the shape of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ at the same flow regime in figures 3(a) and 8(c). This correspondence is because the structured feedback interconnections in (2.4a)–(2.4b) constrain the permissible feedback pathways. In figure 10, (d) $\|\mathcal {H}_{\rho \rho }^S\|_{\infty }$ (associated with input $f_\rho$ and output $\rho$) is significantly larger than the responses shown in the other panels for a flow with $Pr=70$, suggesting the strong role of density in the amplification for this parameter range. We can further isolate each component of the frequency response operator $\mathcal {H}_{\boldsymbol {\nabla }}^S$ by defining

(3.7)\begin{equation} \mathcal{H}_{\boldsymbol{\nabla} ij}^S=\widehat{\boldsymbol{\nabla}}\mathcal{H}_{ij}^S. \end{equation}

The values $\|\mathcal {H}_{\boldsymbol {\nabla } ij}^S\|_{\infty }$, not shown here for brevity, show qualitatively similar behaviour to $\|\mathcal {H}^S_{ij}\|_{\infty }$ as plotted in figures 9 and 10. This componentwise analysis demonstrates that the quasi-horizontal flow structures appearing at high $Pr$ are associated with density fluctuations. The appearance of this type of quasi-horizontal flow structure associated with density fluctuations in a high $Pr$ regime is qualitatively consistent with previous observations that sharp density gradients or even density ‘staircases’ can be observed when $Pr$ is increased (Taylor & Zhou 2017; Zhou et al. 2017b).

Figure 9. Componentwise values of $\log _{10}[\|\mathcal {H}^S_{ij}\|_{\infty }](k_x,k_z)$ at $Re=865$, $Ri_b=0.02$, $Pr=0.7$ ($ij=ux,vy,wz,\rho \rho$). Note the range of the colour bar for each panel is modified based on the maximum value of the particular quantity: (a) $\|\mathcal {H}^S_{u x}\|_{\infty }$; (b) $\|\mathcal {H}^S_{v y}\|_{\infty }$; (c) $\|\mathcal {H}^S_{w z}\|_{\infty }$; (d) $\|\mathcal {H}^S_{\rho \rho }\|_{\infty }$.

Figure 10. Componentwise values of $\log _{10}[\|\mathcal {H}^S_{ij}\|_{\infty }](k_x,k_z)$ at $Re=4250$, $Ri_b=0.02$, $Pr=70$ ($ij=ux,vy,wz$, and $\rho \rho$). Note the range of the colour bar for each panel is modified based on the maximum value of the particular quantity: (a) $\|\mathcal {H}^S_{u x}\|_{\infty }$; (b) $\|\mathcal {H}^S_{v y}\|_{\infty }$; (c) $\|\mathcal {H}^S_{wz}\|_{\infty }$; (d) $\|\mathcal {H}^S_{\rho \rho }\|_{\infty }$.

4. Scaling for density-associated flow structures when $Pr\gg 1$

The previous subsection reveals the appearance of quasi-horizontal flow structures associated with density fluctuations in the $Pr\gg 1$ limit. In this section, we construct an analytical scaling of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ in terms of $Re$ and $Pr$ to provide further evidence that such flow structures prefer the $Pr\gg 1$ regime. The analytical scaling in terms of $Re$ and $Pr$ can also further provide insight into high $Re$ and $Pr$ flow regimes beyond the current computation range achievable through DNS.

We assume streamwise-invariance ($k_x=0$) and unstratified flow ($Ri_b=0$) to facilitate analytical derivation. The importance of streamwise-invariant flow structures is suggested by the quasi-horizontal flow structures ($k_x\approx 0$ and $k_z\approx 0$), which are nearly streamwise constant. The independence with respect to variations in $Ri_b$ of the amplification of streamwise-invariant flow structures $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|^{sc}_{\mu }$ shown in figures 5 and 7 and the analysis of Langham et al. (2020) suggest that $Ri_b=0$ (i.e. density fluctuations can be treated as a passive scalar) is a reasonable regime to consider to obtain further insight. The analytically derived $Re$ and $Pr$ scalings of each component of $\mathcal {H}_{ij}^S$ in (3.5) and $\mathcal {H}_{\boldsymbol {\nabla } ij}^S$ in (3.7) are presented in theorems 4.1(a) and (b), respectively.

Theorem 4.1 Consider streamwise-invariant ($k_x=0$) unstratified ($Ri_b=0$) plane Couette flow with a passive scalar ‘density’ field.

(a) Each component of $\|\mathcal {H}^S_{ij}\|_{\infty }$ ($i=u,v,w,\rho$ and $j=x,y,z,\rho$) scales as

(4.1)\begin{align} & \begin{bmatrix} \|\mathcal{H}^S_{ux}\|_{\infty} & \|\mathcal{H}^S_{uy}\|_{\infty} & \|\mathcal{H}^S_{uz}\|_{\infty} & \|\mathcal{H}^S_{u\rho}\|_{\infty}\\ \|\mathcal{H}^S_{vx}\|_{\infty} & \|\mathcal{H}^S_{vy}\|_{\infty} & \|\mathcal{H}^S_{vz}\|_{\infty} & \|\mathcal{H}^S_{v\rho}\|_{\infty}\\ \|\mathcal{H}^S_{wx}\|_{\infty} & \|\mathcal{H}^S_{wy}\|_{\infty} & \|\mathcal{H}^S_{wz}\|_{\infty} & \|\mathcal{H}^S_{w\rho}\|_{\infty}\\ \|\mathcal{H}^S_{\rho x}\|_{\infty} & \|\mathcal{H}^S_{\rho y}\|_{\infty} & \|\mathcal{H}^S_{\rho z}\|_{\infty} & \|\mathcal{H}^S_{\rho \rho}\|_{\infty}\\ \end{bmatrix} \nonumber\\ &\quad =\begin{bmatrix} Re\, h^S_{ux}(k_z) & Re^2\, h^S_{uy}(k_z) & Re^2\, h^S_{uz}(k_z) & 0\\ 0 & Re\,h^S_{vy}(k_z) & Re\,h^S_{vz}( k_z) & 0\\ 0 & Re\, h^S_{wy}(k_z) & Re\, h^S_{wz}(k_z) & 0\\ 0 & Re^2Pr\,h^S_{\rho y}(k_z) & Re^2Pr\,h^S_{\rho z}(k_z) & RePr\, h^S_{\rho \rho}(k_z) \end{bmatrix}, \end{align}

where functions $h^S_{ij}(k_z)$ are independent of $Re$ and $Pr$.

(b) Each component of $\|\mathcal {H}^S_{\boldsymbol {\nabla } ij}\|_{\infty }$ ($i=u,v,w,\rho$ and $j=x,y,z,\rho$) scales as

(4.2)\begin{align} & \begin{bmatrix} \|\mathcal{H}^S_{\boldsymbol{\nabla} ux}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla} uy}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla} uz}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla} u\rho}\|_{\infty}\\ \|\mathcal{H}^S_{\boldsymbol{\nabla} vx}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla} vy}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla} vz}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla} v\rho}\|_{\infty}\\ \|\mathcal{H}^S_{\boldsymbol{\nabla} wx}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla} wy}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla} wz}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla} w\rho}\|_{\infty}\\ \|\mathcal{H}^S_{\boldsymbol{\nabla} \rho x}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla}\rho y}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla}\rho z}\|_{\infty} & \|\mathcal{H}^S_{\boldsymbol{\nabla}\rho \rho}\|_{\infty}\\ \end{bmatrix} \nonumber\\ &\quad =\begin{bmatrix} Re\, h^S_{\boldsymbol{\nabla} ux}(k_z) & Re^2\, h^S_{\boldsymbol{\nabla} uy}(k_z) & Re^2\, h^S_{ \boldsymbol{\nabla} uz}(k_z) & 0\\ 0 & Re\,h^S_{\boldsymbol{\nabla} vy}(k_z) & Re\,h^S_{ \boldsymbol{\nabla} vz}( k_z) & 0\\ 0 & Re\, h^S_{\boldsymbol{\nabla} wy}(k_z) & Re \,h^S_{ \boldsymbol{\nabla} wz}(k_z) & 0\\ 0 & Re^2Pr\,h^S_{\boldsymbol{\nabla}\rho y}(k_z) & Re^2Pr\,h^S_{\boldsymbol{\nabla}\rho z}(k_z) & RePr\, h^S_{\boldsymbol{\nabla}\rho \rho}(k_z) \end{bmatrix}, \end{align}

where functions $h^S_{\boldsymbol {\nabla } ij}(k_z)$ are independent of $Re$ and $Pr$.

The first three columns and three rows presented in (4.1) are the same as those derived in Jovanović (2004, theorem 11) for unstratified wall-bounded shear flows with no passive scalar field. The details of the proof are presented in Appendix B. These results demonstrate that $Pr$ only contributes to the scaling associated with the density field (here of course assumed to be a passive scalar); i.e. the bottom rows of (4.1) and (4.2) corresponding to the density output. We also note that the rightmost columns of (4.1) and (4.2) show that the forcing in the density equation $f_{\rho }$ does not influence the output corresponding to velocity components $u$, $v$ and $w$, which is consistent with the assumption that $Ri_b=0$, in that the density perturbation behaves as a passive scalar.

The effect of imposing a componentwise structure for the nonlinearity within the feedback is analogous to the effect seen in unstratified PCF (Liu & Gayme 2021, § 3.3). The imposed correlation between each component of the modelled forcing $f_{x,\xi }$, $f_{y,\xi }$, $f_{z,\xi }$, $f_{\rho,\xi }$, and the respective velocity and density components $u$, $v$, $w$, $\rho$ constrain the influence of the forcing to its associated component of the velocity or density field. Thus, the overall scaling of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ is related to the worst-case scaling of the diagonal terms in (4.2) in theorem 4.1. The concept is formalized in theorem 4.2, and we relegate the details of the proof to Appendix B.

Theorem 4.2 Given a wavenumber pair $(k_x,k_z)$

(4.3)\begin{equation} \|\mathcal{H}_{\boldsymbol{\nabla}}^S\|_{\mu}\geq \max[\|\mathcal{H}_{\boldsymbol{\nabla} ux}^S\|_{\infty}, \|\mathcal{H}_{\boldsymbol{\nabla} vy}^S\|_{\infty}, \|\mathcal{H}_{\boldsymbol{\nabla} wz}^S\|_{\infty}, \|\mathcal{H}_{\boldsymbol{\nabla} \rho \rho}^S\|_{\infty}]. \end{equation}

We can combine results in theorems 4.1(b) and 4.2 to obtain the scaling of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ in corollary 4.3

Corollary 4.3 Consider streamwise-invariant ($k_x=0$) unstratified ($Ri_b=0$) PCF with a passive scalar ‘density’ field

(4.4)\begin{equation} \|\mathcal{H}_{\boldsymbol{\nabla}}^S\|_{\mu}(0,k_z)\geq \max[Re\, h^S_{\boldsymbol{\nabla} ux }(k_z), Re\, h^S_{\boldsymbol{\nabla} vy }(k_z),Re\, h^S_{\boldsymbol{\nabla} wz }(k_z), RePr\, h^S_{\boldsymbol{\nabla} \rho \rho}(k_z)], \end{equation}

where functions $h^S_{\boldsymbol {\nabla } ij}(k_z)$ with $ij=ux,vy,wz,\rho \rho$ are independent of $Re$ and $Pr$.

Although corollary 4.3 provides a lower bound on $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$, the numerical results suggest that $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ follows the same $Re$ and $Pr$ scaling as the right-hand side of (4.4) in corollary 4.3. For example, corollary 4.3 suggests that the lower bound of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }(0,k_z)$ will scale as ${\sim }Re$ at a fixed $Pr$, which is consistent with the red dashed lines of figure 5(b). At a fixed $Re$, corollary 4.3 also suggests that $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }(0,k_z)\sim Pr$ in the limit $Pr\gg 1$, but $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }(0,k_z)$ will become independent of $Pr$ in the limit $Pr\ll 1$. This is also consistent with the numerical results shown in the red dashed lines of figure 7(b) that suggest $\|\mathcal {H}_{\boldsymbol {\nabla }}^{S}\|^{sc}_{\mu }\sim Pr$ when $Pr\gg 1$ and independently of $Pr$ when $Pr\ll 1$. For $Pr\gg 1$, theorem 4.2 and corollary 4.3 further suggest that the component $\|\mathcal {H}^{S}_{\boldsymbol {\nabla } \rho \rho }\|_{\infty }$ associated with the density will dominate the overall behaviour of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$, which is consistent with the large amplification of quasi-horizontal flow structures associated with density fluctuations, i.e. $\|\mathcal {H}^S_{\rho \rho }\|_\infty$ shown in figure 10(p). Corollary 4.3 further supports the notion that the flow structures associated with density fluctuations prefer the flow regime with $Pr\gg 1$ under the assumptions of streamwise-invariant ($k_x=0$) and unstratified ($Ri_b=0$) flow.

5. Conclusions and future work

In this paper, we have extended the SIOA originally developed for unstratified wall-bounded shear flows (Liu & Gayme 2021) to stratified PCF. We first apply SIOA to characterize highly amplified flow structures in the intermittent regimes of stratified PCF. We examine how variations in $Re$ and $Ri_b$ affect flow structures with $Pr=0.7$. SIOA predicts the characteristic wavelengths and angle of the oblique turbulent bands observed in very large channel size DNS of the low-$Re$ low-$Ri_b$ intermittent regime of stratified PCF (Deusebio et al. 2015; Taylor et al. 2016). In the high-$Re$ high-$Ri_b$ intermittent regime, SIOA identifies quasi-horizontal flow structures resembling turbulent–laminar layers (Deusebio et al. 2015).

Having validated the ability of the SIOA approach to predict important structures in the intermittent regime, we next investigate the behaviour of the flow across a range of important control parameters $Re$, $Ri_b$ and $Pr$. Increasing $Ri_b$ is shown to reduce the amplification of streamwise-varying flow structures. The results indicate that the classical marginally stable $Ri_b= 1/4$ for the laminar base flow appears to be associated with a change in the most amplified flow structures, an observation which is robust for a wide range of $Re$ and valid at $Pr\approx 1$.

We then examine flow behaviour at different $Ri_b$ and $Pr$. For flows with $Pr\ll 1$, a larger value of $Ri_b$ is required to reduce the amplification of streamwise-varying flow structures to the same level as streamwise-invariant ones compared with flows with $Pr\approx 1$. The largest amplification also occurs at the same value of $Pr Ri_b$, consistent with the observation of matching averaged density profile for flows with the same value of $PrRi_b$ in the $Pr\ll 1$ regime (Langham et al. 2020). For flows with $Pr\gg 1$, the SIOA identifies another quasi-horizontal flow structure that is independent of $Ri_b$. By decomposing input–output pathways into each velocity and density component, we show that these quasi-horizontal flow structures for flows with $Pr\gg 1$ are associated with density fluctuations. The importance of this density-associated flow structure for flows with $Pr\gg 1$ is further highlighted through a derived analytical scaling of amplification with respect to $Re$ and $Pr$ under the assumptions that the flow is streamwise invariant ($k_x=0$) and unstratified (i.e. $Ri_b=0$ and the density behaves as a passive scalar). The above observations using SIOA distinguish two types of quasi-horizontal flow structures, one emerging in the high-$Re$ high-$Ri_b$ regime and the other one (associated with density fluctuations) emerging in the high $Pr$ regime.

The results here suggest the promise of this computationally tractable approach in identifying horizontal length scales of prominent flow structures in stratified wall-bounded shear flows and opens up many directions for future work. For example, this framework may be extended to other stratified wall-bounded shear flows such as stratified channel flow (Garcia-Villalba & del Alamo 2011), stratified open channels (Flores & Riley 2011; Brethouwer, Duguet & Schlatter 2012; Donda et al. 2015; He & Basu 2015, 2016) and the stratified Ekman layer (Deusebio et al. 2014), where intermittent regimes of flow dynamics were also observed. This framework may also be extended to configurations where the background density gradient and velocity gradient are orthogonal, e.g. spanwise stratified PCF (Facchini et al. 2018; Lucas, Caulfield & Kerswell 2019) and spanwise stratified plane Poiseuille flow (Le Gal et al. 2021).