Layer formation in horizontally forced stratified turbulence: connecting exact coherent structures to linear instabilities

Abstract

We consider turbulence in a stratified ‘Kolmogorov’ flow, driven by horizontal shear in the form of sinusoidal body forcing in the presence of an imposed background linear stable stratification in the third direction. This flow configuration allows the controlled investigation of the formation of coherent structures, which here organise the flow into horizontal layers by inclining the background shear as the strength of the stratification is increased. By numerically converging exact steady states from direct numerical simulations of chaotic flow, we show, for the first time, a robust connection between linear theory predicting instabilities from infinitesimal perturbations to the robust finite-amplitude nonlinear layered state observed in the turbulence. We investigate how the observed vertical length scales are related to the primary linear instabilities and compare to previously considered examples of shear instability leading to layer formation in other horizontally sheared flows.

1 Introduction

Laboratory experiments, numerical simulations and even field measurements of turbulent flows which are strongly stably stratified are regularly observed to exhibit spontaneous layering where the density field is organised into relatively deep, relatively well-mixed regions or ‘layers’ separated by relatively thin ‘interfaces’ with enhanced density gradients (Park, Whitehead & Gnanadeskian 1994; Holford & Linden 1999a , b ; Oglethorpe, Caulfield & Woods 2013; Falder, White & Caulfield 2016; Leclercq et al. 2016; Thorpe 2016). Scaling analyses (Billant & Chomaz 2001; Lindborg 2006) have provided some theoretical basis for the expected behaviour of vertical ‘layer’ scales relative to the basic parameters involved in the distinguished asymptotic limit of extremely ‘strong’ stratification and intense turbulence. At its heart this scaling has the central idea that sufficiently strong stratification (definable in a precise fashion) inevitably introduces anisotropy into the velocity field: vertical velocities are suppressed by the buoyancy force compared to horizontal velocity components, thus leading to pancake-like layering, with characteristic turbulent regions having much larger horizontal extent $l_{h}$ than vertical extent $l_{v}$ . Indeed, following Falder et al. (2016), we refer to this regime as the ‘layered anisotropic stratified turbulence’ (LAST) regime.

Even in a purely one-dimensional model, where there is no characteristic horizontal scale $l_{h}$ , and so the LAST regime is formally not possible, the hypothesis that sufficiently strong stratification suppresses vertical motions can lead to a prediction of layering in the density field. Specifically, if the vertical velocity is suppressed, it is at least plausible that some appropriately averaged vertical (turbulent) buoyancy flux should decrease with sufficiently strong stratification. As originally argued by Phillips (1972), if there is a range of stratifications for which the vertical (turbulent) buoyancy flux decreases with increasing stratification, local perturbations in density gradient will tend to be intensified rather than smoothed out by turbulent mixing, suggesting that uniform density gradients are ‘unstable’, in that they are prone to developing a layer–interface structure (see Park et al. (1994) for a clear discussion). Although in its simplest formulation (where for sufficiently strong stratification the buoyancy flux decreases monotonically with stratification) this Phillips mechanism is ill-posed, corresponding essentially to an ‘anti-diffusive’ problem, various regularisation mechanisms to limit the ‘sharpness’ of the interfacial density gradients have been proposed. For example, Barenblatt et al. (1993) demonstrated that the underlying problem could become well-posed if there was a time-lag between the turbulence and the mixing irreversibly modifying the density distribution, and there is at least some evidence that just such a time-lag exists in transient turbulent mixing driven by shear instabilities (Mashayek, Caulfield & Peltier 2013). Alternatively, Balmforth, Llewellyn Smith & Young (1998) proposed that the relationship between buoyancy flux and stratification should be ‘N-shaped’, with a return to an increase in buoyancy flux with increasing and sufficiently large stratification. Indeed, the possibly non-monotonic dependence of irreversible buoyancy flux on external parameters is a very active area of research controversy (see e.g Maffioli, Brethouwer & Lindborg 2016; Venaille, Gostiaux & Sommeria 2016; Venayagamoorthy & Koseff 2016). Importantly however, the fundamental physical mechanisms leading to either the formation or the maintenance of layered density distributions are still quite open.

One highly promising possible mechanism is suggested by the linear instability of a vortex dipole in a stratified environment, known as the zigzag instability (Billant & Chomaz 2000b ). The linear theory of the instability provides a scaling for layer depth which has been confirmed at finite amplitude numerically and experimentally (Billant & Chomaz 2000a , b , c ; Deloncle, Billant & Chomaz 2008; Waite & Smolarkiewicz 2008; Augier, Billant & Chomaz 2015) resulting in the zigzag instability being a popular explanation for the observation of layers (Thorpe 2016). A significant question is therefore how generic is a ‘zigzag’ mechanism? Specifically, given other horizontally varying base flows, constituting other, less precisely organised distributions of vertical vorticity, do analogous linear instabilities exist to provide vertical structure? In addition, is it possible to make a more robust connection between such a linear stability mechanism and a highly nonlinear, yet identifiable sustained turbulent state? For example Deloncle et al. (2008) did not observe nonlinear saturation in their simulations of the zigzag instability of counter-rotating vortex pairs.

Here we consider the case when forcing provides a horizontal shear which resembles in at least some respects the case of (vertically) stratified Taylor–Couette flow considered by Oglethorpe et al. (2013), Leclercq et al. (2016) (minus rotation and curvature, but where both non-monotonic buoyancy flux with stratification and spontaneous layer formation is known to occur) and the vertically invariant base flows of Billant & Chomaz (2000b ) and Deloncle, Chomaz & Billant (2007). We are further motivated by the results of Basak & Sarkar (2006), who considered the freely decaying case of a horizontal shear layer in a stratified environment. They found that, compared to its vertically sheared counterpart, the flow exhibits more intense turbulence as the stratification does not penalise the initial two-dimensional linear instability of the shear layer. They also observed ‘dislocated pancake vortices’, in that the flow exhibits vertical structure, but low vertical velocity, at strong stratification, precisely as postulated for the Billant & Chomaz scaling and the Phillips mechanism. Of interest here is whether a sustained horizontal shear gives rise to sustained density layering and how the flow is organised in such a situation. Of course, it is important to remember that horizontal shear also has great relevance to oceanic flows, where zonal jets provide horizontal shear but are also observed to develop vertical structure in the form of ‘stacked jets’ (Hua, Moore & LeGentil 1997; Eden & Dengler 2008).

We force our model system via a sinusoidal body-forcing term, leading to what is known as Kolmogorov flow (Arnold & Meshalkin 1960), a flow known to support a rich array of spatiotemporal behaviour (Lucas & Kerswell 2014). Stratified Kolmogorov flow has also been studied in two dimensions when the shear is oriented vertically (see Balmforth & Young (2002, 2005), who consider the linear instabilities, weakly nonlinear theory and direct numerical simulation of the flow). Here some layering connected to a stratified conductive instability is observed when the Prandtl number (i.e. the ratio of the diffusivity of the density field to the kinematic viscosity) is small. The three-dimensional extension of this work was recently discussed by Garaud, Gallet & Bischoff (2015), who, motivated by astrophysical systems, consider the low Péclet number case. They investigate the limits of linear and energy stability and find that strong enough stratification will eventually suppress all instabilities.

Our approach here is to make use of recent developments in (unstratified) shear flow transition where unstable exact coherent structures (ECSs) have been computed and shown to be responsible for organising the transition to and sustenance of turbulence (Kawahara & Kida 2001; van Veen, Kida & Kawahara 2006; Kerswell & Tutty 2007; Viswanath 2007; Cvitanović & Gibson 2010; Kawahara, Uhlmann & van Veen 2012; Kreilos & Eckhardt 2012; Budanur et al. 2017). This is made possible by employing a high-dimensional Newton-GMRES-hookstep algorithm which is able to converge unstable steady and time-periodic solutions efficiently from guesses taken from a chaotic simulation (Chandler & Kerswell 2013; Lucas & Kerswell 2015). This work is, to date and the authors’ best knowledge, only the third work attempting to use such methods to understand the influence of stratification on a shear-driven flow. Olvera & Kerswell (2017) and Deguchi (2017) examine the effect of wall-normal stratification on certain known ECSs in plane Couette flow and find the states are heavily affected by stratification.

The goal in this paper is to use this methodology to extract new solutions which exhibit the relevant coherent structure about which turbulence organises and forms layers in a stratified fluid. Such states are by definition nonlinear and their origins can be investigated by using continuation in parameter space to build a bifurcation diagram and determine if a robust connection can be made to linear instability mechanisms like the zigzag instability. As such this study forms the first comprehensive investigation of layer formation in horizontally driven simulations and the role of nonlinear ECSs in creating such layers. Such a simple basic flow and periodicity in all directions makes the system an efficient setting for the so-called ‘dynamical systems approach’. Fundamentally, our philosophy is to focus on ‘structures not statistics’, as we aim to identify the characteristic, inherently nonlinear structures associated with (robust) layer formation in a stratified flow subjected to horizontal shear.

The rest of the paper is organised as follows. In § 2, we present the formulation of the flow which we are considering. In § 3, we describe the results of our direct numerical simulations, where we do indeed observe the development of layered structures. In § 4, we present the results of a linear stability analysis of our flows, showing that there is a clear connection between the linear instabilities of the horizontally sheared stratified Kolmogorov flow and the linear three-dimensional instability studied by Deloncle et al. (2007). We then identify some exact coherent structures in § 5 and demonstrate that they can be connected both to the observed zigzag structures in the nonlinear simulations and the predicted linear instabilities. We discuss our results and draw conclusions in § 6.

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 (2001) 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. (2016), 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. (2016), 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 (2017). 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 (2001) 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. 2007) 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. (2016), $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 (2016)) 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. (2007), 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 (2001), 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 (2001).

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 (1999a ), Waite & Bartello (2004), Lindborg (2006), Brethouwer et al. (2007), Khani & Waite (2013), Oglethorpe et al. (2013), Augier et al. (2015), Thorpe (2016)). 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 (2000a , b , c ). Since this class of instability is also often invoked as the precursor to layer formation (see e.g. Thorpe 2016), 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 1986; van Veen & Goto 2016). 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. (2007). 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 (2001).

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 2011). We find very good agreement with the results of Deloncle et al. (2007) 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. (2007), 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 (2001) and Deloncle et al. (2007). In particular, see figure 4 in Deloncle et al. (2007), 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 2007). 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 (2007) and the implementation here is identical to that described in detail in Chandler & Kerswell (2013), 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 (2013)) 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 (2000a ).

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. (2008) 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. (2007) (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.

6 Discussion

In this paper we have attempted to lay out the instability and bifurcation mechanisms by which the density field is spontaneously arranged into sustained layers by stratified turbulence driven with a horizontal shear. The key step is that we have been able to connect the coherent structures, observed to organise the mean flow of the turbulence into inclined shear layers, to the stratified linear instabilities of the basic horizontally varying and vertically uniform velocity profile. This is achieved by obtaining, directly from direct numerical simulations, various nonlinear coherent structures as exact, yet unstable, steady states. Specifically, we have found two steady states which are striking representations of the mean flows observed during the turbulence, one state having a single vertical layer ( (a1–a3) in figure 9) and the other state having two layers ( (A) in figure 10). By constructing a bifurcation diagram using arclength continuation in flow parameters, both are found to originate from a sequence of instabilities, the principal or primary of which is a stratified linear instability of the base flow.

We have shown that this sequence of instabilities is responsible for breaking various symmetries of the flow. The bifurcations that we have isolated first break the $z$ and $x$ continuous invariance of the base flow through the linear instability studied in § 4. The ensuing primary states retain a discrete translational invariance $(x,y,z)\rightarrow (x+\unicode[STIX]{x03C0}/k_{x},y,z+\unicode[STIX]{x03C0}/k_{l})$ , where $k_{l}$ is the vertical wavenumber of the instability. It is the breaking of this symmetry which leads to the formation of the chevron pattern with the same vertical wavelength as the primary state, i.e. $2\unicode[STIX]{x03C0}/k_{l}$ observed here. This is the salient distinction between primary states (b) or (C) and the states (a1)–(a3) or (A). Importantly, our results do not exclude the possibility that other points of linear instability of the basic velocity profile could lead directly to a nonlinear saturated chevron state of the same vertical wavelength if $k_{x}=0$ at the bifurcation point or with half the vertical wavelength if $k_{x}\neq 0$ at the bifurcation point. However, we have just not observed such primary bifurcations to chevron states here.

It should be emphasised that while the DNS calculations presented here fall at the periphery of the asymptotic scaling regimes and also do not exhibit the anisotropy between horizontal and vertical length scales traditionally considered when discussing ‘stratified turbulence’, particularly the LAST regime, this flow still shows significant modification of its mean profiles when, e.g. $Re_{B}$ is small and $F_{h}$ is relatively large. For example, case A3, where $Re_{B}=6.8$ and $F_{h}=0.074,$ shows very distinctive inclined layers in its mean flow (figure 2 and table 1). We should therefore stress that even when turbulence is considered to be weakly affected by stratification, the nonlinear flow response to buoyancy may still be significant. Such effects are possibly often masked by complicated stochastic forcing mechanisms (Waite & Bartello 2004; Rorai, Mininni & Pouquet 2014; Maffioli et al. 2016).

We expect the scenario outlined above to be quite general. For example, we conjecture the forcing wavenumber $n$ is not likely to change our analysis greatly, other than to adjust the linear stability regimes slightly and obviously add additional horizontal structure. The linear instability reported here is clearly closely related to the three-dimensional stratified instabilities identified by Billant & Chomaz (2000c ) for dipoles, and Deloncle et al. (2007) for hyperbolic tangent shear layers and Bickley jets. Importantly, we have, for the first time, identified a family of nonlinear steady states which do indeed manifest a chevron-shaped or zigzag vertical structure, very reminiscent of previously observed finite-amplitude structures for a vortex dipole (Billant & Chomaz 2000a , b ). Also in contrast to Billant & Chomaz (2000c ) and Deloncle et al. (2007) we have carried out a viscous analysis, and find that the scalings apparently set by an inviscid instability mechanism are still observable at finite amplitude and in situations where viscosity is playing a significant role. It is clearly of interest to investigate how robust this observation is, as the effects of non-trivial viscosity on stratified turbulence is often complex and subtle (see for example the discussion in Khani & Waite (2016)). Turbulence forced with these other horizontal profiles, therefore, may be expected to behave similarly, being organised about unstable steady states with vertically inclined shear layers. How such structures behave as the imposed horizontal shear becomes progressively more localised in space is an interesting and relevant research problem.

The effect of inclination of initial shear relative to the background stratification has been considered in the freely decaying case by Jacobitz & Sarkar (1998). They found that even a small inclination angle is sufficient to increase turbulent production significantly. The effect of inclination angle on the linear instability of a Bickley jet was considered by Candelier, Le Dizès & Millet (2011), who found that growth rates are proportional to the sine of the angle. We therefore conjecture that a body-forced inclined shear will generate projected versions of the coherent structures discovered here, eventually vanishing as the shear approaches a purely vertical orientation.

A further connection should be made to flows which are not driven by a large-scale mean shear, but are known to exhibit layer formation, including the $U/N$ scaling for layer depth. Examples include oscillating grid or rod turbulence (Browand, Guyomar & Yoon 1987; Holford & Linden 1999a , b ) or the extensive literature on wakes, see for example Spedding (2002), Diamessis, Spedding & Domaradzki (2011), Spedding (2014). In such cases it is not uncommon to observe turbulent flows developing layers at relatively large Froude numbers ( $O(1),$ though of course caution must be exercised in the comparison of specific numerical values of parameters defined in different ways in different studies). Furthermore, at the moment, the precise mechanisms at play are still unclear. Several authors have discussed the possibility of a zigzag mechanism in these flows (Spedding 2002; Thorpe 2016), which, given our findings that nonlinear saturated states underlying the layer pattern formation have their ultimate origins in such instabilities, increases the possibilities that other, more complicated base flows than the one considered here, may also be subject to such behaviours. Note that in the case of wakes some energy injection would be required, for example forcing with a vertically and horizontally varying streamwise flow, to locate steady nonlinear ECSs analogous to those presented here and connect them to a linear stability mechanism.

One open question is how precisely the nonlinear solutions obtained correspond to episodes in the direct numerical simulations at other parameter values, in particular at higher values of the Reynolds number. We do not address this question here, as our principal aim is to demonstrate that flow structures of interest can not only be located but also their source (through an identifiable sequence of instabilities and bifurcations) can be determined. It is reasonable to expect that, given sufficient numerical resources, corresponding states at a desired point in parameter space could also be found and we conjecture that a similar picture will emerge there. However, it is also plausible that the coherent structures observed in one particular part of parameter space are mere ‘ghosts’ of solutions existing in a nearby region of parameter space, and therefore solution efforts would fail to converge at that particular part of parameter space. For example, $\unicode[STIX]{x1D6FC}=1$ appears to be inaccessible to most of the solution state curves described here, yet the coherent structures are still observed.

Significantly, the secondary and tertiary instabilities of the exact states we have identified support turbulence which rearranges the density field into well-defined mixed regions and interfaces which are sustained for long times. This aspect, and the mixing properties of this flow will be the subject of a subsequent paper. A fundamental question in stratified turbulence has been how, or indeed if, a stably stratified buoyancy field becomes organised into layers and interfaces with relatively deep mixed regions separated by relatively thin regions with strong gradients, and in particular how the vertical length scales of such layers might be set. Here we have shown, for the first time, a robust connection which extends all the way from the linear instability of a simple laminar state, to the nonlinear saturated layered state. These nonlinear states which have been isolated are in striking qualitative agreement with the mean flows which are observed in the ‘turbulent flow’ (in the sense that the flow is chaotic and exhibits a range of scales) at larger $Re$ . As such it places the proposed scaling for the layer depth $l_{v}\sim U/N$ (in terms of a characteristic velocity scale $U$ and characteristic buoyancy frequency $N$ ) on a much stronger foundation. The conventional but formally unjustified assumption that linear theory can explain some aspect of a highly nonlinear turbulent flow is actually here given some justification. However, it is still unclear how the basin of attraction for the turbulence is continuously deformed towards the observed layered configurations in phase space as the stratification is increased and how this relates to the appearance of the new unstable invariant manifolds (generated via the instabilities) about which the turbulence organises.