2 Formulation

We begin by considering the following version of the monochromatic body-forced, incompressible, Boussinesq equations

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}^{\ast }}{\unicode[STIX]{x2202}t^{\ast }}+\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\ast }\boldsymbol{u}^{\ast }+\frac{1}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x1D735}^{\ast }p^{\ast }=\unicode[STIX]{x1D708}\unicode[STIX]{x0394}^{\ast }\boldsymbol{u}^{\ast }+\unicode[STIX]{x1D712}\sin (2\unicode[STIX]{x03C0}ny^{\ast }/L_{y})\hat{\boldsymbol{x}}-\frac{\unicode[STIX]{x1D70C}^{\ast }g}{\unicode[STIX]{x1D70C}_{0}}\hat{\boldsymbol{z}}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x2202}t^{\ast }}+\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D70C}^{\ast }+\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D70C}_{B}=\unicode[STIX]{x1D705}\unicode[STIX]{x0394}^{\ast }\unicode[STIX]{x1D70C}^{\ast } & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\boldsymbol{u}^{\ast }=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}^{\ast }(x,y,z,t)=u^{\ast }\hat{\boldsymbol{x}}+v^{\ast }\hat{\boldsymbol{y}}+w^{\ast }\hat{\boldsymbol{z}}$ is the three-dimensional velocity field, $n$ is the forcing wavenumber, $\unicode[STIX]{x1D712}$ the forcing amplitude, $\unicode[STIX]{x1D708}$ is the kinematic viscosity, $\unicode[STIX]{x1D705}$ is the density diffusivity, $p^{\ast }$ is the pressure, $\unicode[STIX]{x1D70C}_{0}$ is an appropriate reference density and $\unicode[STIX]{x1D70C}^{\ast }(x,y,z,t)$ the varying part of the density away from the background linear density profile $\unicode[STIX]{x1D70C}_{B}=-\unicode[STIX]{x1D6FD}z,$ i.e. $\unicode[STIX]{x1D70C}_{total}=\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x1D70C}_{B}(z)+\unicode[STIX]{x1D70C}^{\ast }(x,y,z,t)$ . Gravity acts in the negative $z$ -direction. We impose periodic boundary conditions in all directions $(x,y,z)\in [0,L_{x}]\times [0,L_{y}]\times [0,L_{z}]$ on $\boldsymbol{u}^{\ast }$ and $\unicode[STIX]{x1D70C}^{\ast }.$

Figure 1. Schematic showing the base flow and the background stratification.

For simplicity we set $L_{f}=L_{y}=L_{z}$ . The system is naturally non-dimensionalised using the characteristic length scale $L_{f}/2\unicode[STIX]{x03C0}$ , characteristic time scale $\sqrt{L_{f}/2\unicode[STIX]{x03C0}\unicode[STIX]{x1D712}}$ and density gradient scale $\unicode[STIX]{x1D6FD}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}_{B}\cdot \hat{\boldsymbol{z}}$ to give

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}p=\frac{1}{Re}\unicode[STIX]{x0394}\boldsymbol{u}+\sin (ny)\hat{\boldsymbol{x}}-B\unicode[STIX]{x1D70C}\hat{\boldsymbol{z}}, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}=w+\frac{1}{RePr}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C} & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where we define the Reynolds number $Re$ , a buoyancy parameter $B$ , the Prandtl number $Pr$ and the aspect ratio of the domain $\unicode[STIX]{x1D6FC}$ as

(2.7a-d ) $$\begin{eqnarray}Re:=\frac{\sqrt{\unicode[STIX]{x1D712}}}{\unicode[STIX]{x1D708}}\left(\frac{L_{y}}{2\unicode[STIX]{x03C0}}\right)^{3/2},\quad B:=\frac{g\unicode[STIX]{x1D6FD}L_{y}}{\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D712}2\unicode[STIX]{x03C0}}=\frac{N_{B}^{2}L_{y}}{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D712}}=\frac{1}{4\unicode[STIX]{x03C0}^{2}F_{hB}^{2}},\qquad Pr:=\frac{\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D705}},\qquad \unicode[STIX]{x1D6FC}=\frac{L_{f}}{L_{x}},\end{eqnarray}$$

where $N_{B}$ is the (dimensional) buoyancy frequency associated with the background density field, and the characteristic velocity in this particular (background) definition of the horizontal Froude number $F_{hB}=U/(L_{f}N_{B})$ has been constructed using the characteristic length scale divided by the characteristic time scale, exactly as in the definition of the Reynolds number.

It is important to note that the distinguished limit formally required of the scaling proposed by Billant & Chomaz (Reference Billant and Chomaz2001) is that the Reynolds number is large (and so it is expected that the flow is turbulent), while an appropriate horizontal Froude number $F_{h}$ is small such that $ReF_{h}^{2}$ is still large. For our flows, if we use this background definition of $F_{hB}$ , these conditions correspond to both $Re$ and (much more stringently) $Re/(4\unicode[STIX]{x03C0}^{2}B)$ being large. However, as discussed for example in Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016), there are two significant, inter-related issues which are of interest. First, in practice, it is important to know the limits of applicability of a formally asymptotic scaling, identifying the finite limits for the various parameters to be sufficiently large or small so that the predicted regime actually arises. Second, and related to this, there is of course a freedom (not least in terms of the selection of factors of $2\unicode[STIX]{x03C0}$ for example) in the specific form of the various characteristic length scales, and so comparison of specific delimiting numerical values of these non-dimensional parameters quoted in different studies must be done with care. As discussed in more detail by Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016), parameters based on the internal properties of any ensuing turbulent flow are more straightforward to compare from one flow to another (or indeed from one subregion of a flow to another), and so we will also consider alternative measures for the key scaling of the internal Froude number in terms of the internal properties of the turbulence.

The equations then are solved over the cuboid $[0,2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}]\times [0,2\unicode[STIX]{x03C0}]^{2}$ . We define the diagnostics involved in the energetic budgets as

(2.8a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{K}}={\textstyle \frac{1}{2}}\langle |\boldsymbol{u}|^{2}\rangle _{V},\qquad {\mathcal{I}}=\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{f}\rangle _{V}=\langle u\sin (ny)\rangle _{V}, & \displaystyle\end{eqnarray}$$

(2.9a-c ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{B}}=\langle \boldsymbol{u}\boldsymbol{\cdot }B\unicode[STIX]{x1D70C}\hat{\boldsymbol{z}}\rangle _{V}=\langle wB\unicode[STIX]{x1D70C}\rangle ,\qquad {\mathcal{D}}=\frac{1}{Re}\langle |\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\rangle _{V},\qquad {\mathcal{D}}_{lam}=\frac{Re}{2n^{2}}, & \displaystyle\end{eqnarray}$$

where ${\mathcal{K}}$ is the total kinetic energy density, ${\mathcal{I}}$ is the energy input by the forcing, ${\mathcal{B}}$ is the buoyancy flux, and ${\mathcal{D}}$ is the dissipation rate. ${\mathcal{D}}_{lam}$ is the dissipation rate associated with the basic state

(2.10) $$\begin{eqnarray}\boldsymbol{u}_{lam}=\frac{Re}{n^{2}}\sin (ny)\hat{\boldsymbol{x}},\end{eqnarray}$$

where the forcing and dissipation precisely balance and $\langle (\cdot )\rangle _{V}:=\unicode[STIX]{x1D6FC}\int \!\!\int \!\!\int (\cdot )\,\text{d}x\,\text{d}y\,\text{d}z/(2\unicode[STIX]{x03C0})^{3}$ denotes a volume average. We consider only $n=1$ throughout, i.e. the flow is forced with $\sin (y)\hat{\boldsymbol{x}}$ and denote a time average with an overbar, i.e. $\bar{(\cdot )}=[\int _{0}^{T}(\cdot )\text{d}t]/T$ , where $T$ is normally the full simulation time (having removed the initial 5 time units of transient spin-up from the initial condition). Vorticity $\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D735}\times \boldsymbol{u}$ is used as the prognostic variable and direct numerical simulations (DNS) are performed using the fully dealiased (two-thirds rule) pseudospectral method with mixed fourth-order Runge–Kutta and Crank–Nicolson timestepping implemented in CUDA to run on GPU cards. This code is a further extension, to stratified flow, of that used in Lucas & Kerswell (Reference Lucas and Kerswell2017). Each simulation is run on one GPU card and due to memory limitations on the current generation of NVIDIA chips we are restricted to resolutions of $N_{x}=N_{y}=N_{z}=256$ . We initialise the flow velocity field’s Fourier components with uniform amplitudes and randomised phases in the range $2.5\leqslant |\boldsymbol{k}|\leqslant 9.5$ such that the total enstrophy $\langle |\boldsymbol{\unicode[STIX]{x1D714}}|^{2}\rangle _{V}=1$ and leave the density field unperturbed initially i.e. $\unicode[STIX]{x1D70C}=0$ . Spatial convergence is checked by comparing the Kolmogorov microscale $\unicode[STIX]{x1D702}=\left(Re^{3}{\mathcal{D}}\right)^{-1/4}$ with the smallest permitted scale in the system, i.e. the maximum wavenumber in simulation $k_{max}=N_{x}/3=85$ , in other words we ensure $k_{max}\unicode[STIX]{x1D702}>1$ .

3 Direct numerical simulations

To determine the effect of stratification on this flow we begin by performing direct numerical simulations with various $B$ at fixed $Re=100$ , $Pr=1$ , $n=1$ , $\unicode[STIX]{x1D6FC}=1$ , integrating until $T=100.$ These calculations are designed to span a range of $B$ from the essentially unstratified to the point at which well-formed layers are observed. They serve primarily as motivation for the stability and exact coherent structure analysis to follow. It is important to stress that the key aim of this paper is not to consider properties of flows in the LAST regime with extreme parameter values, but rather to explore the connections between linear instabilities, exact coherent structures and layer formation in a sheared and stratified fluid.

Figure 2 shows $yz$ planes (at $x=\unicode[STIX]{x03C0}$ ) of $t=50$ snapshots of $\unicode[STIX]{x1D70C}$ and $u$ along with $\langle \bar{\unicode[STIX]{x1D70C}}\rangle _{x}$ and $\langle \bar{u}\rangle _{x},$ i.e. streamwise averages of the mean (in time over the whole simulation). Immediately we recognise the emergence of coherent structures organising the flow. Most evident in the mean, but also noticeable in the snapshots, is an inclining of the background horizontal shear into characteristic chevron, or (rotated) ‘v’ shapes. Associated with this structure is the organisation of the density field into vertical layers, the number of layers (and associated v-shapes in $u$ ) increasing with $B$ . As the stratification is increased, the mean flow becomes progressively more layered. The formation mechanism of these layers, and in particular their relationship to both linear instabilities and exact coherent structures, will be discussed in the subsequent sections.

Figure 2. Snapshots in a $yz$ midplane (at $x=\unicode[STIX]{x03C0}$ for simulations with $Re=100$ , $Pr=1$ , $n=1$ and $\unicode[STIX]{x1D6FC}=1$ ) of $\unicode[STIX]{x1D70C}$ (a, e, i, m) and $u$ (b, f, j, n) at $t=50.$ (c, g, k, o) show $\langle \bar{\unicode[STIX]{x1D70C}}\rangle _{x}$ and (d, h, l, p) show $\langle \bar{u}\rangle _{x},$ i.e. the $x$ -average of the time-averaged (across the full $t\in [0,T]$ interval with $T=100$ ) perturbation density and streamwise velocity. Rows show $B=1,5,10,50$ (simulations A1–A4 from table 1) from top to bottom. Notice the increase in vertical structure as $B$ increases.

At $B\geqslant 50$ the simulations exhibit significant bursts between quiescence and turbulence due to the combination of stratification suppressing vertical shear instability and body forcing across the domain. To obtain statistical stationarity we apply a throttling method to modulate the forcing strength to maintain turbulence and thus observe the formation of layers at increased stratification; without the throttling, bursting overturns the density field intermittently and does not permit well-defined layers to form. The results from large $B$ in the throttled cases are summarised in figures 3 and 4. The important observation is the persistence, at large stratification, of the coherent chevron structures that are responsible for layer formation. Further details for the cases shown in figures 3 and 4 are given in the appendix, as well as a detailed discussion of the bursting phenomena and our throttling protocol.

Figure 3. Snapshots in a $yz$ midplane (at $x=\unicode[STIX]{x03C0}$ ) of $\unicode[STIX]{x1D70C}$ (a, e, i, m) and $u$ (b, f, j, n) at $t=25$ . (c, g, k, o) show $\langle \bar{\unicode[STIX]{x1D70C}}\rangle _{x}$ and (d, h, l, p) show $\langle \bar{u}\rangle _{x}$ , i.e. the $x$ -average of the time-averaged perturbation density and streamwise velocity. Time averages are over $t\in [0,T]$ with $T=50$ . Rows show data from simulations D1, B3, C2 and C3 (as defined in table 1) from top to bottom. Notice the increase in the number of layers as the buoyancy parameter $B$ increases and the increasingly angled structure of the velocity field. All simulations have $Pr=1.$

Figure 4. Three-dimensional rendering of: (a) the streamwise velocity $u$ ; (b) the vertical shear $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z$ ; and (c) the total density field $\unicode[STIX]{x1D70C}_{B}(z)+\unicode[STIX]{x1D70C}$ , for the throttled simulation C2 with parameters as listed in table 1. Note that regions of increased density gradient coincide once again with quiescent regions in $u$ where the vertical shear is minimum.

3.1 Layer length scale

For the throttled simulations, it is now possible to observe smaller vertical length scales in the developing layered structures. The increased buoyancy flux due to the sustained turbulence in the inclined shear layers now spontaneously leads to sharpening interfaces separating relatively well-mixed ‘layers’. This spontaneous interface–layer formation can be most clearly seen in figure 5 (a) where profiles of the total density (including the background linear component $\unicode[STIX]{x1D70C}_{B}(z)$ ) are plotted. Figure 4 shows three-dimensional renderings of snapshots of the streamwise velocity, vertical shear $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z$ , and the total density field for the simulation C2 (see table 1) with $B=1000$ and ${\mathcal{D}}_{0}=100$ . It can be seen that gradients of density are approximately coincident with planes of small vertical shear of the underlying coherent structure, where the vertical buoyancy flux is small, or more precisely away from regions of strong vertical shear.

Figure 5. (a) Profiles of total density $\unicode[STIX]{x1D70C}_{tot}=\unicode[STIX]{x1D70C}_{B}+\unicode[STIX]{x1D70C}$ at $t=50$ , $x=y=\unicode[STIX]{x03C0}$ for cases B1, B3 and C1 from table 1. (b) Vertical length scale $l_{\bar{u}}$ as defined in (3.2) plotted against $U/N_{B}=\langle |\bar{(u)}|\rangle _{V}/\sqrt{B}$ for each DNS in table 1. We find the scaling $l_{v}\sim U/N_{B}$ is recovered, particularly for smaller $U/N_{B}$ . A least-squares linear fit is also plotted: $0.78l_{\bar{u}}-0.12.$

A common characterisation of vertical length scales in stratified turbulent flows is the scaling proposed by Billant & Chomaz (Reference Billant and Chomaz2001) leading to the ‘layered anisotropic stratified turbulence’ (LAST) regime mentioned in the introduction. The characteristic vertical scale of the layers $l_{v}\sim U/N$ is commonly observed in simulation and experiments, where $U$ is a typical velocity scale and $N$ is the appropriate buoyancy frequency. As mentioned in the introduction, this scaling is strictly valid in a distinguished limit of small horizontal Froude number, $F_{h}=U/l_{h}N\rightarrow 0,$ implying that $l_{v}\ll l_{h}$ . Due not least to the periodicity of our computational domain, the formal requirement that $l_{v}\ll l_{h}$ cannot be satisfied in our computations as the layers typically span the full horizontal extent of the box. However, under the assumption that the dissipation rate has an inertial scaling (see Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007) the LAST regime may equivalently be associated with sufficiently large values of the buoyancy Reynolds number $Re_{B}$ , defined as

(3.1) $$\begin{eqnarray}Re_{B}=\frac{\unicode[STIX]{x1D716}}{\unicode[STIX]{x1D708}N_{B}^{2}}=\frac{{\mathcal{D}}Re}{B}\end{eqnarray}$$

( $\unicode[STIX]{x1D716}$ is the dimensional dissipation rate). As discussed in more detail in Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016), $Re_{B}$ is determined from the internal local characteristics of the turbulence (captured by the dissipation rate) relative to the joint stabilising effects of the fluid’s viscosity and the background density gradient.

Although it is now becoming widely accepted (see, for example, Venayagamoorthy & Koseff (Reference Venayagamoorthy and Koseff2016)) that $Re_{B}$ alone cannot uniquely identify the properties of stratified turbulence driven by shear, the evidence is nevertheless strong that $Re_{B}\gtrsim O(10)$ is necessary for there to be any opportunity for the LAST regime to occur. As explained in detail by Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), this is due to the necessity for there to be large-scale separations between three scales: the Kolmogorov microscale, $\unicode[STIX]{x1D702}=(Re^{3}{\mathcal{D}})^{-1/4}$ ; the Ozmidov scale $l_{O}=({\mathcal{D}}/B^{3/2})^{1/2}$ and some vertical scale of the layers $l_{v}$ , which must be smaller than the vertical extent of the computational domain. The Ozmidov scale is the largest vertical scale which is largely unaffected by the background stratification, and so the ordering $\unicode[STIX]{x1D702}\ll l_{O}<l_{v}$ ensures that there is both an inertial range for the turbulent motions largely unaffected by the stratification, and a range of scales (between $l_{O}$ and $l_{v}$ ) where the LAST regime occurs. Since by definition $Re_{B}=(l_{O}/\unicode[STIX]{x1D702})^{4/3}$ , $Re_{B}$ must be sufficiently large for an inertial range to exist.

The formal accessibility of this hierarchy of length scales within our simulations is challenging. Due to the horizontal (periodic) boundary conditions, there is no natural way in which the actual horizontal extent of any layer can be identified. Specifically it is not appropriate to use the characteristic length scale of the shear forcing as the horizontal scale of the ‘layer’, as the flow has clearly adjusted so that the (density) layer extends horizontally across the entire computational domain. Furthermore, the observed vertical length scale of layers must be sufficiently small so that it can be measured within the vertical extent of the computational domain, yet still sufficiently large so that it is both resolved computationally and also larger than the Ozmidov length scale $l_{O}$ . In turn, $l_{O}$ must still be sufficiently large so that there is separation between it and the Kolmogorov microscale, $\unicode[STIX]{x1D702}$ , which, for the entire simulation to be resolved adequately, must be at worst of the same order as the smallest scale of the simulation. These various length scales are listed in table 1. It is important to appreciate that several of the simulations, particularly when $B$ is large, have values of $Re_{B}$ which are too small for the flow to be in the LAST regime. In particular such small values of $Re_{B}$ can arise when the dissipation rate ${\mathcal{D}}$ in the definition (3.1) is taken to be the time average of the actual (in general time-varying) volume-averaged dissipation rate within our simulations. That being said, the parameter choices for these DNS were purposely made to span situations from very weakly stratified cases, to cases where significant deformation of the mean flows are observed due to the stratification leading subsequently to the formation of layers.

It is therefore still of interest to investigate whether the spontaneous layers observed here are consistent with the scaling ideas at the heart of this proposed regime, even if the relevant parameters are not in the formally necessary self-consistent asymptotic scale separation. We choose to define an appropriate vertical ‘layer’ length scale $l_{\bar{u}}$ for the layers as

(3.2) $$\begin{eqnarray}l_{\bar{u}}=\frac{2\unicode[STIX]{x03C0}\int \hat{\bar{u}}\,\text{d}k_{z}}{\int k_{z}\hat{\bar{u}}\,\text{d}k_{z}},\end{eqnarray}$$

where $\hat{\bar{u}}$ is the Fourier transform of the horizontally and temporally averaged ( $t\in [0,T]$ ) streamwise velocity $\langle \bar{u}\rangle _{x,\,y}$ . For the LAST regime as discussed above, the layer length scale should scale with the ratio of an appropriate large scale or background horizontal velocity scale $U$ and the (dimensionless) buoyancy frequency $N_{B}=\sqrt{B}$ from (2.7). An appropriate choice for the velocity scale $U$ is the volume and time-averaged absolute streamwise velocity $U=\langle |\bar{u}|\rangle _{V}$ . Figure 5 shows $l_{\bar{u}}$ plotted against $U/N_{B}$ . We find that for small $U/N_{B},l_{\bar{u}}$ approaches a linear trend, and is much more scattered as $U/N_{B}$ reaches larger values, though still largely following the asymptotic prediction of Billant & Chomaz (Reference Billant and Chomaz2001), which is formally expected to occur as $U/(L_{h}N_{B})\rightarrow 0$ . Significantly, the scaling appears to be satisfied even for the flows with relatively small buoyancy Reynolds number (see table 1) associated with relatively large values of $B$ and hence $N_{B}$ , where the LAST regime scaling arguments are not strictly valid. As is shown from the least-squares fit, the vertical Froude number $F_{v}$ , defined as

(3.3) $$\begin{eqnarray}F_{v}=\frac{U}{N_{B}l_{u}}\simeq 0.78=O(1),\end{eqnarray}$$

consistently with the scaling regime of Billant & Chomaz (Reference Billant and Chomaz2001).

It should also be noted that similar estimates of layer depth could be constructed using the density field. However, due to the strong spatiotemporal intermittency at large $B$ the trend to small $U/N_{B}$ is harder to observe in the density field, since measurable layering of density requires sustained buoyancy fluxes and mixing, which are challenging to resolve. At higher resolutions, where turbulence can be maintained with larger target dissipation rates at larger $B$ , it is perfectly reasonable to expect that robust layering of the density field obeying the $U/N_{B}$ scaling will be more straightforward to observe, as is strongly suggested in figure 4. Note several other studies of stratified turbulence have observed this $U/N$ layer scaling, (e.g. Holford & Linden (Reference Holford and Linden1999a ), Waite & Bartello (Reference Waite and Bartello2004), Lindborg (Reference Lindborg2006), Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), Khani & Waite (Reference Khani and Waite2013), Oglethorpe et al. (Reference Oglethorpe, Caulfield and Woods2013), Augier et al. (Reference Augier, Billant and Chomaz2015), Thorpe (Reference Thorpe2016)). However, a key open question remains surrounding the precise mechanisms for the formation of such layers, which will be the focus of the subsequent sections. In particular, while a connection to known linear instabilities has been conjectured as the mechanism for such layers, there does not yet exist a robust connection between the output of (inherently nonlinear) direct numerical simulations and such an instability.

4 Linear stability analysis

There appears to be at least qualitative similarity between the streamwise velocity structure observed in our forced flows and the streamwise velocity eigenstructure at onset of the now classic ‘zigzag’ instability mechanism of Billant & Chomaz (Reference Billant and Chomaz2000a ,Reference Billant and Chomaz b ,Reference Billant and Chomaz c ). Since this class of instability is also often invoked as the precursor to layer formation (see e.g. Thorpe Reference Thorpe2016), it seems appropriate to investigate the linear stability properties of the flows we are considering, in particular to identify whether they are prone to instabilities which may be identified as being of ‘zigzag’ type. For the case when the forcing scale and streamwise integral scale are the same (i.e. using our notation $\unicode[STIX]{x1D6FC}=1$ and $n=1$ ), the unstratified three-dimensional Kolmogorov flow is known to undergo subcritical transition to turbulence with the base flow remaining linearly stable at all $Re$ (Marchioro Reference Marchioro1986; van Veen & Goto Reference van Veen and Goto2016). Any linear instability in this geometry is therefore due to the added physics provided by the statically stable stratification.

The inviscid linear stability properties of a background horizontal $\tanh$ shear profile with linear background stratification was considered by Deloncle et al. (Reference Deloncle, Chomaz and Billant2007). They reported that, while two-dimensional (in the horizontal plane) homogeneous perturbations are the most unstable, leading to instabilities of KH type associated with the inflection point in the background shear profile, new inherently three-dimensional stratified instabilities are present and can have comparable growth rates. The growth rates of the stratified instabilities are also observed to follow the self-similar scaling with respect to an appropriate horizontal Froude number as proposed by Billant & Chomaz (Reference Billant and Chomaz2001).

Figure 6. Variation of the maximum (across all horizontal wavenumbers, associated with $k_{x}=0.59$ ) and scaled growth rate $\unicode[STIX]{x1D70E}_{m}/Re$ with: vertical wavenumber $k_{z}$ (a); and vertical wavenumber scaled with stability horizontal Froude number (as defined in (4.2))  $F_{hS}k_{z}$ (b). Red pluses mark results for a flow with $Re=100$ , $B=10^{4}$ ( $F_{hS}=1$ ), green crosses mark results for a flow with $Re=1000$ , $B=10^{8}$ ( $F_{hS}=0.1$ ), red stars mark results for a flow with $Re=1000,$ $B=10^{6}$ ( $F_{hS}=1$ ) and magenta squares mark results for a flow with $Re=5000$ , $B=10^{10}$ ( $F_{hS}=0.05$ ).

In the conventional fashion, we consider normal mode disturbances proportional to $\exp [i(k_{x}x+k_{z}z)+\unicode[STIX]{x1D70E}t]$ , away from the base state $S_{B}=(\boldsymbol{U}_{B},\unicode[STIX]{x1D70C}_{B})$ defined by

(4.1a,b ) $$\begin{eqnarray}\boldsymbol{U}_{B}=\frac{Re}{n^{2}}\sin (ny)\hat{\boldsymbol{x}},\quad \unicode[STIX]{x1D70C}_{B}=-z,\end{eqnarray}$$

and solve the ensuing eigenvalue problem for the linearised equations using the Python NUMPY package, which is itself a front end for LAPACK (van der Walt, Colbert & Varoquaux Reference van der Walt, Colbert and Varoquaux2011). We find very good agreement with the results of Deloncle et al. (Reference Deloncle, Chomaz and Billant2007) for the horizontal sinusoidal shear flow considered here provided $Re\gg 1$ , which is a natural requirement as their analysis is inviscid. Figure 6 shows the variation with vertical wavenumber $k_{z}$ of the maximal growth rate $\unicode[STIX]{x1D70E}_{m}$ (across all streamwise wavenumbers, corresponding to $k_{x}=0.59$ ) rescaled with $Re$ (since $\unicode[STIX]{x1D70E}_{m}\sim Re$ ), plotted for a range of $Re$ and $B.$ Similarly to the results of Deloncle et al. (Reference Deloncle, Chomaz and Billant2007), we also find that two-dimensional disturbances (with $k_{z}=0$ ) are most unstable and that the most unstable horizontal wavenumber $k_{x}$ is independent of $B.$ This is unsurprising, as this two-dimensional instability is once again of KH type associated with the inflection point in the background velocity shear.

For further comparison with their results, it is necessary to define an appropriate horizontal Froude number for this stability problem. An appropriate velocity scale for comparison for these linear instabilities is the maximum magnitude of $\boldsymbol{U}_{B}$ i.e. $Re/n^{2}$ (half the total velocity jump across the shear layer), while an appropriate length scale is (within our non-dimensionalisation) $1/n$ , and so the stability horizontal Froude number $F_{hS}$ may be defined as

(4.2) $$\begin{eqnarray}F_{hS}=\frac{Re}{n\sqrt{B}}.\end{eqnarray}$$

As already noted, for all the calculations we present here, we set $n=1$ . Replotting the growth rates against $F_{hS}k_{z}$ leads to a very good collapse of the curves for a wide range of stability horizontal Froude numbers $F_{hS}=(0.05,0.1,1)$ as shown in figure 6(b), once again in pleasing agreement with Billant & Chomaz (Reference Billant and Chomaz2001) and Deloncle et al. (Reference Deloncle, Chomaz and Billant2007). In particular, see figure 4 in Deloncle et al. (Reference Deloncle, Chomaz and Billant2007), where an equivalent definition of horizontal Froude number (i.e. using the maximum magnitude of the velocity distribution and the scale of the shear layer) is used and the curves are extremely similar to those in figure 6; the shape of the shear profile serves only to shift the growth rates slightly. We stress again that this is the appropriate Froude number for comparison of linear stability results, and not as the characteristic horizontal $F_{h}$ of the layers which develop in the DNS, as they extend horizontally across the entire computational domain.

In the (nonlinear) DNS, only horizontal wavenumbers $k_{x}=m/\unicode[STIX]{x1D6FC}$ (where $m$ is an integer) are admissible. Therefore, any linear instabilities must have $k_{x}\geqslant 1/\unicode[STIX]{x1D6FC}$ to have any chance of occurring within our computational domain. In figure 7(a) for flows with $B=50$ we plot neutral curves of linear stability for various vertical wavenumbers $k_{z}$ on the $Re-k_{x}$ plane. Although the unstratified long-wave instability is clearly evident for $k_{z}=0$ , the neutral curve asymptotes near $k_{x}=1$ but crucially never crosses as $Re$ increases. Similarly, the neutral curve for $k_{z}=1$ always has $k_{x}<1,$ i.e. a long streamwise wavelength is required in this range of $Re$ , and so neither the instabilities with $k_{z}=0$ nor $k_{z}=1$ are expected to arise in the direct numerical simulations discussed above, with aspect ratio $\unicode[STIX]{x1D6FC}=1$ . However, the neutral curves for the instabilities with $k_{z}>1$ do intersect $k_{z}=1$ (and indeed cross to even higher wavenumbers), suggesting that these instabilities can actually develop within our direct numerical simulations where $\unicode[STIX]{x1D6FC}=1$ . The plot in figure 7(b) shows the neutral curves on a $B-k_{x}$ plane, for fixed $Re=15$ , illustrating that increasing $B$ and $Re$ introduces more instabilities with higher values of vertical wavenumber $k_{z}$ .

Figure 7. (a) Neutral stability curves on the $Re-k_{x}$ plane for modes with vertical wavenumbers $k_{z}=0,1,2,3,4$ in a flow with $B=50$ ; (b) neutral stability curves on the $B-k_{x}$ plane for modes with vertical wavenumbers $k_{z}=0,1,2,3,4$ in a flow with $Re=15$ . Note that the $k_{z}=0$ instability has $k_{x}<1$ for all $Re$ and $B$ .

5 Exact coherent structures

We can determine if the finite-amplitude states observed in the direct numerical simulations can be connected to the inherently stratified linear instabilities discussed in the previous section by attempting to converge invariant states using a Newton-GMRES-hookstep algorithm (Viswanath Reference Viswanath2007). The approach used here is simply to sample the time series of a certain DNS which exhibits a close approach to the coherent structure of interest and output the full state vector for post-processing. The state vector in this case consists of the full flow field and density degrees of freedom;

(5.1) $$\begin{eqnarray}\boldsymbol{X}=\left[\begin{array}{@{}c@{}}u\\ v\\ w\\ \unicode[STIX]{x1D70C}\end{array}\right].\end{eqnarray}$$

The post-processing takes the form of a high-dimensional root-finding algorithm which solves

(5.2) $$\begin{eqnarray}\boldsymbol{F}(\boldsymbol{X}_{0},T_{p}):=\boldsymbol{X}(\boldsymbol{X}_{0},T_{p})-\boldsymbol{X}_{0}=0,\end{eqnarray}$$

where $\boldsymbol{X}_{0}$ is the starting condition and $\boldsymbol{X}(\boldsymbol{X}_{0},T_{p})$ is the final state after some small, fixed period trajectory of length $t=T_{p}$ . This period $T_{p}$ is arbitrary in the cases considered here since searches are only conducted for steady states. Note that the final state vector is a highly nonlinear function of the starting state $\boldsymbol{X}_{0}.$ Each Newton iteration updates $\boldsymbol{X}_{0}^{new}=\boldsymbol{X}_{0}+\unicode[STIX]{x1D6FF}\boldsymbol{X}_{0}$ , where $\unicode[STIX]{x1D6FF}\boldsymbol{X}_{0}$ is given by

(5.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{F}}{\unicode[STIX]{x2202}\boldsymbol{X}_{0}}\unicode[STIX]{x1D6FF}\boldsymbol{X}_{0}=-\boldsymbol{F}(\boldsymbol{X}_{0},T_{p}).\end{eqnarray}$$

The Newton-GMRES-hookstep method then forms the solution to the ensuing linear system of derivatives $\unicode[STIX]{x2202}\boldsymbol{F}/\unicode[STIX]{x2202}\boldsymbol{X}_{0}$ via GMRES, a Krylov subspace method and a simple forward difference to approximate the action of the matrix. A hookstep then constrains each Newton step to a trust region inside which the Newton linearisation is expected to hold. The method used was first described by Viswanath (Reference Viswanath2007) and the implementation here is identical to that described in detail in Chandler & Kerswell (Reference Chandler and Kerswell2013), except for the inclusion of the GPU Boussinesq timestepping code and additional degrees of freedom for density. Note in general we search for translations in the solution (as explained in Chandler & Kerswell (Reference Chandler and Kerswell2013)) due to the continuous symmetries in $x$ and $z$ which allow travelling waves to be converged. However, all the states under discussion here are stationary. Successfully converging invariant states and performing arclength continuation of the solutions in a certain control parameter (for example $\unicode[STIX]{x1D6FC}$ , $Re$ , or $B$ ) enables us to build a bifurcation diagram for these invariant states. This continuation method includes the relevant parameter into the state vector $\boldsymbol{X}$ as an unknown and the solution is extrapolated along the solution branch via a local arclength parameterisation which permits the solution curves to turn corners.

In practice, the greatest difficulty in this case (and in general) is in finding a starting $\boldsymbol{X}_{0}$ sufficiently nearby in phase space to an underlying unstable steady solution for the Newton method to converge. For this reason we concentrated on unthrottled flows with low $Re$ in relatively long domains with $\unicode[STIX]{x1D6FC}<1$ when attempting to identify invariant states. In such situations, the flow is less ‘unstable’, in the specific sense that relatively simple coherent structures are observed to be approached closely. Therefore, we expect our attempts at convergence to meet with more success. A set of new simulations were performed at a range of low $Re$ , and at various values of $B$ , and state vectors were sampled and fed into the Newton-GMRES method. Our successes at $Re=8$ and 15 at $B=50$ and $\unicode[STIX]{x1D6FC}=0.9$ are reported here. Figure 8 shows a time series of the dissipation for these cases to highlight that, despite the low $Re,$ these flows are unsteady and, at least for these short times, they exhibit chaotic dynamics. It is quite possible the long-time attractor is periodic or quasi-periodic, particularly at $Re=8$ . This issue is not of immediate interest as we use these chaotic trajectories to sample for a nearby unstable steady solution. We observe the closest approach to the coherent structure when ${\mathcal{D}}/{\mathcal{D}}_{lam}$ has a minimum, and so use this as a guess for the Newton-GMRES algorithm.

Figure 8. Time series of ${\mathcal{D}}/{\mathcal{D}}_{lam}$ for the cases $Re=8$ and $Re=15$ for $B=50$ and $\unicode[STIX]{x1D6FC}=0.9$ .

Figure 9 shows the results arising from converging the unstable invariant state (labelled state (a1)) in a flow with $Re=8$ , $B=50$ and $\unicode[STIX]{x1D6FC}=0.9$ . This state clearly has one chevron-shaped structure in its streamwise velocity. Due to the assumed periodicity in the vertical direction, this state will naturally form a ‘zigzag’, as clearly observed in the experiments of Billant & Chomaz (Reference Billant and Chomaz2000a ).

Figure 9. Bifurcation diagram showing solution branches for flows with $Re=8$ and $B=50$ identified through continuation in $\unicode[STIX]{x1D6FC}$ (or equivalently streamwise wavenumber $k_{x}$ ) against scaled total dissipation $D/D_{lam}$ . The spatial structures of streamwise velocity for states (a1)–(c) are shown in the top left inset. Thick lines correspond to primary states generated by linear instabilities of the laminar velocity profile $u_{lam}$ , (equation (2.10)), thin lines correspond to secondary states generated by linear instabilities of the primary states and dashed lines correspond to tertiary states, in turn arising from linear instabilities of the secondary states. The green curve follows the two-dimensional solution state (c) (with $k_{z}=0$ ) connecting to the long-wave homogeneous instability of the laminar base flow. The top right inset shows a zoom of the bifurcation leading to state (c) and also shows the (secondary) stratified instability from this state leading to the three-dimensional inclined ‘V’-shaped coherent structure (a3) which follows the lighter blue curve as $\unicode[STIX]{x1D6FC}$ is continued, and which then through a further tertiary bifurcation (marked with a dashed red line) switches to state (a1). State (a1) follows the red dashed curve in this plane, and in turn attaches to the secondary state (a2) plotted with a thin orange line. State (a2) finally connects to the primary state (b), plotted with a thick dark blue line, which originates from the $k_{z}=1$ stratified linear instability of the laminar base flow $u_{lam}$ . Note that the labelling of the curves corresponds to the locations in parameter space for the three-dimensional visualisations plotted above.

The lower panel of figure 9 shows the bifurcation diagram mapped out by continuing this solution in $\unicode[STIX]{x1D6FC}$ , projected onto ${\mathcal{D}}/{\mathcal{D}}_{lam}$ such that when this quantity is identically 1, the solution has returned to $u_{lam}$ defined in equation (2.10). Solution (a1) follows the (dashed) red line and is observed to connect, via two bifurcations (and via state (a2)) to the state (b), which is the ensuing state arising from the $k_{z}=1$ linear instability of the laminar state. One can check the bifurcation of state (b) against the linear instability by noticing that its long-wave streamwise dependence means that as the solution curve goes to ${\mathcal{D}}/{\mathcal{D}}_{lam}\rightarrow 1$ the $\unicode[STIX]{x1D6FC}$ parameter can be associated with the $k_{x}$ in the linear stability analysis. In other words, the bifurcation of the nonlinear steady state (b) from $u_{lam}$ at $\unicode[STIX]{x1D6FC}\approx 0.717$ corresponds directly with crossing the $k_{z}=1$ neutral curve in figure 7 at $Re=8$ and $B=50$ and $k_{x}\approx 0.717$ from the linear analysis. Since state (b) arises from a linear instability of the laminar flow, it is marked with a thick line (blue for consistency with the $k_{z}=1$ neutral curves in figure 7) and we refer to it as a ‘primary’ state, while we refer to state (a2) as a ‘secondary’ state, and plot its continuation with a thin line, arising as it does from a bifurcation from a primary state. We call (a1) a ‘tertiary state’ since it bifurcates from a secondary state, and denote it with a dashed line. Note, consistently with the linear stability calculations presented in figure 7, state (b) does not exist in a flow geometry with $\unicode[STIX]{x1D6FC}=1$ , but requires $\unicode[STIX]{x1D6FC}<1$ .

Moreover, when increasing $\unicode[STIX]{x1D6FC}$ , state (a1) connects via a bifurcation to another secondary state (a3), marked with a thin blue line. State (a3) is found to connect to the primary state, labelled (c), and marked with a thick green line (again for consistency with figure 7) arising from the $k_{z}=0$ unstratified linear instability of the laminar flow. As before, the bifurcation at $\unicode[STIX]{x1D6FC}\approx 0.992$ is associated with the corresponding neutral curve for $k_{z}=0,$ at $Re=8,$ $B=50$ and $k_{x}\approx 0.992.$ In other words, the state labelled (c) undergoes its own stratified linear instability (to the state (a3) marked with the light blue line) in much the same fashion as the laminar background flow $u_{lam}$ undergoes the linear instability, ultimately leading to the saturated state (b). The different varieties of chevron-shaped states (a1)–(a3) are very similar and only on close inspection is one able to discern the subtle symmetry differences, notably between (a2) and (a3), which are switched by traversing the ‘rung’ of the tertiary state (a1). Note that all states are inherently three-dimensional apart from the states lying along the $k_{z}=0$ green curve associated with state (c). In summary, we have identified an inherently nonlinear and coherent chevron-shaped structure which is manifest in the states labelled (a). Furthermore, these states are highly reminiscent of the finite-amplitude structures observed experimentally by Billant and Chomaz. Crucially, we establish that the chevron-shaped structure that we have isolated arises from a secondary bifurcation away from the nonlinear state associated with the primary instability, and should not be thought of as the straightforward finite-amplitude manifestation of the primary instability at least for this particular instability with $k_{z}=1$ within our flow domain. This is consistent with the numerically based observations of Deloncle et al. (Reference Deloncle, Billant and Chomaz2008) that the zigzag instability does not ‘saturate’. To be specific, the nonlinear state labelled (b) has a translational invariance $(x,y,z)\rightarrow (x+\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC},y,z+\unicode[STIX]{x03C0})$ , and this symmetry is broken by the bifurcation to the secondary state (a2).

The states shown in figure 9 ultimately arise from the $k_{z}=0$ and $k_{z}=1$ primary instabilities of the laminar base flow, and so are associated with flow geometries with relatively long domains with $\unicode[STIX]{x1D6FC}$ strictly less than one (though in some cases quite close to one). However, this is not a necessary restriction. Figure 10 shows the analogous situation for states bifurcating from the $k_{z}=0$ linear instability (once again plotted with a thick green line for consistency with the neutral curves plotted in figure 7) and from the $k_{z}=2$ linear instability (plotted with a magenta line, also consistently with figure 7). Here, the nonlinear coherent structure is initially found from a guess obtained from the $Re=15,$ $B=50$ and $\unicode[STIX]{x1D6FC}=0.9$ case. The state (denoted (A)) has two chevron-shaped inclined shear layers, as shown in figure 10, and once again is reminiscent of finite-amplitude ‘zigzag’ structures.

Perhaps unsurprisingly, through continuation (plotted with a thin blue line in the figure) this state is once again found to be ‘secondary’, in that it attaches through bifurcations both to the primary state (shown in figure 10 as state (C), and plotted with a thick magenta line for consistency with figure 7) arising from the linear instability with $k_{z}=2$ , and the primary state arising from the two-dimensional ( $k_{z}=0$ ) linear instability, whose continuation with respect to aspect ratio $\unicode[STIX]{x1D6FC}$ is once again plotted with a thick green line. Note that analogously to the $k_{z}=1$ case discussed in figure 9, the secondary state (A) breaks the $(x,y,z)\rightarrow (x+\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC},y,z+\unicode[STIX]{x03C0}/2)$ translational symmetry of the primary state (C).

Using continuation, we also identify another secondary state, labelled (B) and plotted with a thin orange line. This state is, in some sense, lower amplitude, having ${\mathcal{D}}/{\mathcal{D}}_{lam}$ nearer to 1, and does not have such strongly inclined layers, although a chevron-shaped structure is still apparent in the streamwise velocity, as shown in the middle upper panel of figure 10. Importantly however, we find that state (B) still persists to $\unicode[STIX]{x1D6FC}\geqslant 1.$ Notably, the local maximum in the solution curve is found at $\unicode[STIX]{x1D6FC}\approx 0.997$ and forms a smooth maximum, despite looking like a cusp on this $(\unicode[STIX]{x1D6FC},{\mathcal{D}}/{\mathcal{D}}_{lam})$ projection.

It should be noted that state (C) does not actually arise in the first, most unstable eigenvalue crossing of the $k_{z}=2$ linear instability but rather through a secondary appearance of a new unstable direction of $u_{lam}$ . Equivalently, the associated values of $k_{x}$ do not correspond to the crossing of the main neutral curves plotted in figure 7. We therefore conjecture that further bifurcations at the first appearance (as $\unicode[STIX]{x1D6FC}$ increases) will either give rise to yet more coherent nonlinear states, or reconnect to the states discussed here. However, our objective is not to conduct an exhaustive bifurcation analysis, but rather to identify the origins of the converged chevron-shaped states which are of interest due to their clear similarity to the previously reported zigzag structures. Practically, we have also demonstrated the utility of continuation in $\unicode[STIX]{x1D6FC}$ as it conveniently establishes the bifurcation structure of the various solution states.

Figure 10. Bifurcation diagram showing solution state branches for flows with $Re=15$ and $B=50$ identified through continuation in $\unicode[STIX]{x1D6FC}$ (or equivalently streamwise wavenumber $k_{x}$ ) against scaled total dissipation $D/D_{lam}$ . The spatial structures of streamwise velocity for states (A)–(C) are shown at the top. As in figure 9 the thick green curve shows the evolution of the primary state (c) (now at $Re=15$ not 8) arising from the linear instability with $k_{z}=0$ , while the thick magenta curve shows the evolution of primary state (C), arising from the $k_{z}=2$ linear instability of the laminar base flow $u_{lam}$ . The thin orange line shows the evolution of the secondary state (B) and the thin blue line shows the evolution of the secondary state (A), both generated by instabilities of the primary state (C). State (A) also attaches back to the primary state (c) with $k_{z}=0$ . Labels are placed adjacent to where the three-dimensional visualisations are made in parameter space.

By definition all of the states in figures 9 and 10, with the exception of state (c), plotted by the green curves, result from the influence of stratification on the flow. The vertical component of the flow for these states is relatively weak compared to the horizontal component, and as such the displacement of isopycnals is small, even though vertical gradients are significant. We omit figures for the sake of brevity but make the observation that the density layering observed in figures 2–5 requires additional shear instabilities to mix the density field; the flow field snapshots show clear Kelvin–Helmholtz-like overturning events (also see movie in the supplementary material available at https://doi.org/10.1017/jfm.2017.661).

In figure 11, we show how the properties (in particular the streamwise velocity structures) of state (A) change under continuation in the stratification parameter $B$ for $\unicode[STIX]{x1D6FC}$ fixed at its original value of $0.9$ and $Re=15$ . After some turning points, the curve traced by the solution state closes on itself. Although it is at least plausible that continuation in $B$ might lead to some smooth variation of the shear inclination as $B$ increases, it is actually apparent that the increasing vertical structure observed in the numerical simulations arises due to the generation through instability of a family of new invariant solution states. It is somewhat surprising that given the complicated behaviour of the solution state curve there is such little variation in the qualitative shape of the flow state, as visualised in figure 11 by the streamwise velocity structure in a $yz$ midplane at $x=\unicode[STIX]{x03C0}$ . Indeed, the direct numerical simulations show that these coherent structure states are subject to additional pattern forming instabilities, generically leading to localisation of the shear in the $z$ -direction and therefore turbulent motions. Typical dynamical behaviour is shown in figure 12, where three-dimensional visualisations of the streamwise velocity $u$ are shown for the extremely strongly stratified and throttled simulation D5 (with $B=2000$ ) from table 1. How such localisation comes about is a question for future research, but we note in passing, that although the notional value of $Re_{B}$ for this simulation is quite low ( $Re_{B}=4.7$ ), it is still above 1, and the shear instability and layered flow structures bear more than a passing qualitative similarity to those shown in the vigorously turbulent and anisotropic simulation D9.6 of Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) (see in particular their figure 5a).

Figure 11. Bifurcation diagram showing solution state branch (A) from figure 10 continued in $B$ against normalised total dissipation $D/D_{lam}$ for flows at fixed $\unicode[STIX]{x1D6FC}=0.9$ and $Re=15$ . Visualisations of the streamwise velocity in the $yz$ midplane at ( $x=\unicode[STIX]{x03C0}$ ) at various labelled locations on the curve are shown above. Note that at large $B$ the solution curve bends back upon itself to close the loop.

Figure 12. Snapshots of streamwise velocity, $u,$ for simulation D5 ( $Re=50,\,D_{0}=250,\,B=2000$ ) showing the localisation in the vertical of the coherent structures leading to inherently anisotropic and layered turbulent motions.