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

TL;DR

This paper investigates how stratified turbulence in a horizontally forced flow leads to layer formation, linking linear instability predictions to the nonlinear layered states observed in simulations.

Contribution

It establishes a novel connection between linear stability analysis and the nonlinear layered structures in stratified turbulence under horizontal shear.

Findings

Layered structures form as stratification increases.

Linear instability theory predicts the formation of layers.

Vertical length scales relate to primary 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 stratifica- tion 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.

Figures

Tables

Table 1. Table 1: Imposed parameters and diagnostic simulation outputs for the four groups of simulations. Simulations in group A are unthrottled simulations, simulations in group B are throttled simulations with R ​ e = 100 𝑅 𝑒 100 Re=100 , simulations in group C are throttled simulations with R ​ e = 75 𝑅 𝑒 75 Re=75 , and simulations in group D are throttled simulations with R ​ e = 50 𝑅 𝑒 50 Re=50 . The Kolmogorov microscale, within the chosen nondimensionalization is given by η = ( R ​ e 3 ​ 𝒟 ¯ ) − 1 / 4 𝜂 superscript 𝑅 superscript 𝑒 3 ¯ 𝒟 1 4 \eta=\left(Re^{3}\bar{\mathcal{D}}\right)^{-1/4} and is scaled by the maximum wavenumber allowable in the simulation k max = N x / 3 = 85 subscript 𝑘 max subscript 𝑁 𝑥 3 85 k_{\mathrm{max}}=N_{x}/3=85 . The nondimensional Ozmidov length scale l O = ( 𝒟 ¯ / B 3 2 ) 1 / 2 subscript 𝑙 𝑂 superscript ¯ 𝒟 superscript 𝐵 3 2 1 2 l_{O}=\left({\bar{\mathcal{D}}}/{B^{\frac{3}{2}}}\right)^{1/2} , horizontal Froude number F h = D / B ​ U 2 , subscript 𝐹 ℎ 𝐷 𝐵 superscript 𝑈 2 F_{h}={D}/{\sqrt{B}U^{2}}, the buoyancy Reynolds number R ​ e B = 𝒟 ¯ ​ R ​ e / B , 𝑅 subscript 𝑒 𝐵 ¯ 𝒟 𝑅 𝑒 𝐵 Re_{B}={\bar{\mathcal{D}}Re}/{B}, a typical horizontal velocity scale given by U = ⟨ | u ¯ | ⟩ V 𝑈 subscript delimited-⟨⟩ ¯ 𝑢 𝑉 U=\langle|\bar{u}|\rangle_{V} and the Taylor microscale Reynolds number R ​ e λ = 𝒦 turb ​ 10 ​ R ​ e / 𝒟 ¯ 𝑅 subscript 𝑒 𝜆 subscript 𝒦 turb 10 𝑅 𝑒 ¯ 𝒟 Re_{\lambda}=\mathcal{K}_{\mathrm{turb}}\sqrt{{10Re}/{\bar{\mathcal{D}}}} are also listed. Here we define 𝒦 t ​ u ​ r ​ b = ( 1 / 2 ) ​ ⟨ ( 𝒖 − 𝒖 ¯ ) 2 ⟩ V subscript 𝒦 𝑡 𝑢 𝑟 𝑏 1 2 subscript delimited-⟨⟩ superscript 𝒖 ¯ 𝒖 2 𝑉 \mathcal{K}_{turb}=(1/2)\langle(\bm{u}-\bar{\bm{u}})^{2}\rangle_{V} and overbars are always time averages over the full T 𝑇 T window. Group A have T = 100 𝑇 100 T=100 and groups B-D T = 50 . 𝑇 50 T=50.

#	$Re$	$B$	$D_0$	$U$	$\eta k_{\max }$	$l_{\overline{u}}$	$l_o$	$F_h$	$R{e}_B$	$R{e}_{\lambda }$
A1	100	1	-	1.3	2.6	2.5	1.05	0.66	110.9	56.6
A2	100	5	-	1.9	2.7	3.9	0.3	0.13	20.2	48.6
A3	100	10	-	1.7	3.0	3.5	0.15	0.074	6.8	41.8
A4	100	50	-	3.1	2.5	1.8	0.06	0.021	2.9	36.8
B1	100	50	50	7.3	1.03	4.5	0.37	0.126	94.9	178.7
B2	100	100	50	5.9	1.1	2.8	0.18	0.099	34.3	204.0
B3	100	200	50	6.9	1.08	2.1	0.12	0.057	19.3	159.8
C1	75	500	100	6.3	1.2	1.6	0.073	0.068	9.0	263.8
C2	75	1000	100	9.8	1.1	1.2	0.052	0.028	6.4	107.0
C3	75	2000	100	11.5	1.1	1.2	0.03	0.014	3.02	55.6
D1	50	500	100	7.8	1.6	1.7	0.079	0.052	7.1	54.3
D2	50	200	250	11.0	1.2	4.4	0.26	0.115	49.4	89.2
D3	50	500	250	10.8	1.3	2.8	0.11	0.057	14.9	114.2
D4	50	1000	250	10.7	1.3	2.4	0.071	0.044	7.9	121.2
D5	50	2000	250	12.1	1.2	1.4	0.046	0.029	4.7	131.2

Equations

$$\frac{\partial {\mathbf{u}}^{\ast }}{\partial t^{\ast }}+{\mathbf{u}}^{\ast }\cdot {\nabla }^{\ast }{\mathbf{u}}^{\ast }+\frac{1}{{\rho }_0}{\nabla }^{\ast }{p}^{\ast }$$

$$\frac{\partial {\rho }^{\ast }}{\partial t^{\ast }}+{\mathbf{u}}^{\ast }\cdot {\nabla }^{\ast }{\rho }^{\ast }+{\mathbf{u}}^{\ast }\cdot {\nabla }^{\ast }{\rho }_B$$

$${\nabla }^{\ast }\cdot {\mathbf{u}}^{\ast }$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}+\nabla p$$

$$\frac{\partial \rho }{\partial t}+\mathbf{u}\cdot \nabla \rho$$

$$\nabla \cdot \mathbf{u}$$

$$Re:=\frac{\sqrt{\chi }}{\nu }{\left(\frac{L_y}{2\pi }\right)}^{3/2},B:=\frac{g\beta L_y}{{\rho }_0\chi 2\pi }=\frac{N_B^2{L}_y}{2\pi \chi }=\frac{1}{4{\pi }^2{F}_{hB}^2},Pr:=\frac{\nu }{\kappa },\alpha =\frac{L_f}{L_x}.$$

$$\mathcal{K}$$

$$\mathcal{B}$$

$${\mathbf{u}}_{lam}=\frac{Re}{n^2}\sin (ny)\hat{\mathbf{x}}.$$

$$R{e}_B=\frac{ϵ}{\nu N_B^2}=\frac{\mathcal{D}Re}{B}$$

$$l_{\bar{u}}=\frac{2\pi \int \hat{\bar{u}}\mathrm{d}{k}_z}{\int k_z\hat{\bar{u}}\mathrm{d}{k}_z},$$

$$F_v=\frac{U}{N_B{l}_u}\simeq 0.78=\mathcal{O}(1),$$

$${\mathbf{U}}_B=\frac{Re}{n^2}\sin (ny)\hat{\mathbf{x}},{\rho }_B=-z,$$

$$F_{hS}=\frac{Re}{n\sqrt{B}}.$$

$$\mathbf{X}=\left[\begin{array}{c}u\\ v\\ w\\ \rho \end{array}\right]$$

$$\mathbf{F}({\mathbf{X}}_0,T_p):=\mathbf{X}({\mathbf{X}}_0,T_p)-{\mathbf{X}}_0=0,$$

$$\frac{\partial \mathbf{F}}{\partial {\mathbf{X}}_0}\delta {\mathbf{X}}_0=-\mathbf{F}({\mathbf{X}}_0,T_p)$$

$$R{i}_G(y,z,t)=\frac{-B{⟨\frac{\partial {\rho }_{tot}}{\partial z}⟩}_x}{{⟨{\left(\frac{\partial u_h}{\partial z}\right)}^2⟩}_x},$$

$$\mathbf{f}=\chi (t)\sin (ny)\hat{\mathbf{x}}.$$

$$\frac{\mathrm{d}\mathcal{K}}{\mathrm{d}t}=\mathcal{I}+\mathcal{B}-\mathcal{D}$$

$$\chi (t)=\frac{{\mathcal{D}}_0-\mathcal{B}}{⟨u\sin (y)⟩}.$$

Full text

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 et al., 1994; Holford & Linden, 1999a, b; Oglethorpe et al., 2013; Thorpe, 2016; Falder et al., 2016; Leclercq et al., 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 et al., 2013). Alternatively, Balmforth et al. (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 Venayagamoorthy & Koseff (2016); Venaille et al. (2016); Maffioli et al. (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 zig-zag 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 et al., 2008; Waite & Smolarkiewicz, 2008; Augier et al., 2015) resulting in the zig-zag instability being a popular explanation for the observation of layers (Thorpe, 2016). A significant question is therefore how generic is a ‘zig-zag’ 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 zig-zag 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 et al. (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’ (Eden & Dengler, 2008; Hua et al., 1997).

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 (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 et al. (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 et al., 2006; Kerswell & Tutty, 2007; Viswanath, 2007; Cvitanović & Gibson, 2010; Kreilos & Eckhardt, 2012; Kawahara et al., 2012; Willis & Short, 2015). 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 zig-zag 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 section 2, we present the formulation of the flow which we are considering. In section 3, we describe the results of our direct numerical simulations, where we do indeed observe the development of layered structures. In section 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 section 5 and demonstrate that they can be connected both to the observed zig-zag structures in the nonlinear simulations and the predicted linear instabilities. We discuss our results and draw conclusions in section 6.

2 Formulation

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

[TABLE]

where $\bm{u}^{*}(x,y,z,t)=u^{*}\mathbf{\hat{x}}+v^{*}\mathbf{\hat{y}}+w^{*}\mathbf{\hat{z}}$ is the three-dimensional velocity field, $n$ is the forcing wavenumber, $\chi$ the forcing amplitude, $\nu$ is the kinematic viscosity, $\kappa$ is the density diffusivity, $p^{*}$ is the pressure, $\rho_{0}$ is an appropriate reference density and $\rho^{*}(x,y,z,t)$ the varying part of the density away from the background linear density profile $\rho_{B}=-\beta z,$ i.e. $\rho_{\mathrm{total}}=\rho_{0}+\rho_{B}(z)+\rho^{*}(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 $\bm{u}^{*}$ and $\rho^{*}.$

For simplicity we set $L_{f}=L_{y}=L_{z}$. The system is naturally non-dimensionalised using the characteristic length scale $L_{f}/2\pi$, characteristic time scale $\sqrt{L_{f}/2\pi\chi}$ and density gradient scale $\beta=-\nabla\rho_{B}\cdot\mathbf{\hat{z}}$ to give

[TABLE]

where we define the Reynolds number $Re$, a buoyancy parameter $B$, the Prandtl number $Pr$ and the aspect ratio of the domain $\alpha$ as

[TABLE]

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\pi^{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\pi$ 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 sub-region 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\pi/\alpha]\times[0,2\pi]^{2}$. We define the diagnostics involved in the energetic budgets as

[TABLE]

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

[TABLE]

where the forcing and dissipation precisely balance and $\langle(\cdot)\rangle_{V}:=\alpha\int\!\!\int\!\!\int(\cdot)\,dxdydz/(2\pi)^{3}$ denotes a volume average. We consider only $n=1$ throughout, i.e. the flow is forced with $\sin(y)\mathbf{\hat{x}}$ and denote a time average with an overbar, i.e. $\bar{(\cdot)}=[\int_{0}^{T}(\cdot){\rm 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 $\bm{\omega}=\nabla\times\bm{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\leq|\bm{k}|\leq 9.5$ such that the total enstrophy $\langle|\bm{\omega}|^{2}\rangle_{V}=1$ and leave the density field unperturbed initially i.e. $\rho=0$. Spatial convergence is checked by comparing the Kolmogorov microscale $\eta=\left(Re^{3}\mathcal{D}\right)^{-1/4}$ with the smallest permitted scale in the system, i.e. the maximum wavenumber in simulation $k_{\mathrm{max}}=N_{x}/3=85,$ in other words we ensure $k_{\mathrm{max}}\eta>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,$ $\alpha=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=\pi$) of $t=50$ snapshots of $\rho$ and $u$ along with $\langle\bar{\rho}\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.

At $B\geq 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.

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 (left panel) where profiles of the total density (including the background linear component ${\rho_{B}}(z)$) are plotted. Figure 4 shows three-dimensional renderings of snapshots of the streamwise velocity, vertical shear $\frac{\partial u}{\partial 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.

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\to 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

[TABLE]

($\epsilon$ 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 stabilizing 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, $\eta=(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 $\eta\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}/\eta)^{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, $\eta,$ 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 (12) 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

[TABLE]

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 (8). 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

[TABLE]

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 ‘zig-zag’ 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 ‘zig-zag’ type. For the case when the forcing scale and streamwise integral scale are the same (i.e. using our notation $\alpha=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).

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

[TABLE]

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 et al., 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 $\sigma_{m}$ (across all streamwise wavenumbers, corresponding to $k_{x}=0.59$) rescaled with $Re$ (since $\sigma_{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 $\bm{U}_{B}$ i.e. $Re/n^{2}$, (half the total velocity jump across the shear layer) while an appropriate length scale is (within our nondimensionalization) $1/n$, and so the stability horizontal Froude number $F_{hS}$ may be defined as

[TABLE]

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 the right panel of figure 6, 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) are 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}=\frac{m}{\alpha}$ (where $m$ is an integer) are admissible. Therefore, any linear instabilities must have $k_{x}\geq 1/\alpha$ to have any chance of occuring within our computational domain. In figure 7 (left panel) 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 $\alpha=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 $\alpha=1$. The right plot in figure 7 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}$.

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;

[TABLE]

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

[TABLE]

where $\bm{X}_{0}$ is the starting condition and $\bm{X}(\bm{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 $\bm{X}_{0}.$ Each Newton iteration updates $\bm{X}_{0}^{new}=\bm{X}_{0}+\delta\bm{X}_{0}$ where $\delta\bm{X}_{0}$ is given by

[TABLE]

The Newton-GMRES-hookstep method then forms the solution to the ensuing linear system of derivatives $\partial\bm{F}/\partial\bm{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 arc-length continuation of the solutions in a certain control parameter (for example $\alpha,$ $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 $\bm{X}$ as an unknown and the solution is extrapolated along the solution branch via a local arc-length 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 $\bm{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 $\alpha<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 $\alpha=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 9 shows the results arising from converging the unstable invariant state (labelled state (a1)) in a flow with $Re=8,$ $B=50$ and $\alpha=0.9$. This state clearly has one chevron-shaped structure in its streamwise velocity as shown in the top left panel. Due to the assumed periodicity in the vertical direction, this state will naturally form a ‘zig-zag’, as clearly observed in the experiments of Billant & Chomaz (2000a).

The lower panel of figure 9 shows the bifurcation diagram mapped out by continuing this solution in $\alpha$, 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 (11). 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. Specifically, as here the streamwise dependence of the steady states is a domain scale long wavelength, the $\alpha$ of the continuation and the $k_{x}$ of the linear instability mode in figure 7 are equivalent as $\mathcal{D}/\mathcal{D}_{lam}\to 1.$ In other words, the bifurcation of the nonlinear steady state (b) from $u_{lam}$ at $\alpha\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 $\alpha=1$, but requires $\alpha<1$.

Moreover when increasing $\alpha$, 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 $\alpha\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 zig-zag instability does not ‘saturate’. To be specific, the nonlinear state labelled (b) has a translational invariance $(x,y,z)\rightarrow(x+\pi/\alpha,y,z+\pi)$ 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 $\alpha$ 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 $\alpha=0.9$ case. The state (denoted (A)) has two chevron-shaped inclined shear layers, as shown in the upper leftmost panel of figure 10, and once again is reminiscent of finite amplitude ‘zig-zag’ 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 the upper righthand panel of 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 $\alpha$ 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+\pi/\alpha,y,z+\pi/2)$ translational symmetry of 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 $\alpha\geq 1.$ Notably, the local maximum in the solution curve is found at $\alpha\approx 0.997$ and forms a smooth maximum, despite looking like a cusp on this $(\alpha,\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 $\alpha$ 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 zig-zag structures. Practically, we have also demonstrated the utility of continuation in $\alpha$ as it conveniently establishes the bifurcation structure of the various solution states.

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).

In figure 11, we show how the properties (in particular the streamwise velocity structures) of the state (A) change under continuation in the stratification parameter $B$ for $\alpha$ 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 $(y,z)$ midplane at $x=\pi$. 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).

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-3) in figure 9) and the other state having two layers ( (A) in figure 10). By constructing a bifurcation diagram using arc-length 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 section 4. The ensuing primary states retain a discrete translational invariance $(x,y,z)\rightarrow(x+\pi/k_{x},y,z+\pi/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\pi/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 (Maffioli et al., 2016; Waite & Bartello, 2004; Rorai et al., 2014).

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 zig-zag 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 nontrivial 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 et al. (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 (Holford & Linden, 1999a, b; Browand et al., 1987) or the extensive literature on wakes, see for example Spedding (2002); Diamessis et al. (2011); Spedding (2014). In such cases it is not uncommon to observe turbulent flows developing layers at relatively large Froude numbers ($\mathcal{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 zig-zag 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, $\alpha=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 toward 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.

Acknowledgements. We extend our thanks, for many helpful and enlightening discussions, to Paul Linden, John Taylor, Stuart Dalziel and the rest of the ‘MUST’ team in Cambridge and Bristol. We also thank the three anonymous referees whose constructive comments have significantly improved the clarity of the manuscript. The source code used in this work is provided at https://bitbucket.org/dan_lucas/psgpu and the associated data including initialisation files and converged states is found at https://doi.org/10.17863/CAM.13090 This work is supported by EPSRC Programme Grant EP/K034529/1 entitled ‘Mathematical Underpinnings of Stratified Turbulence’. The majority of the research presented here was conducted when DL was a postdoctoral researcher in DAMTP as part of the MUST programme grant.

Appendix A DNS

A.1 Recurrent bursts

One important feature of the inclined shear, observed from the DNS and figure 2, is that the mean flow now has a component of vertical shear. We find that for simulations with fixed $Re$ and large $B$ (i.e. for $B\geq 50$) vertical shear instability (of Kelvin-Helmholtz (KH) type) is suppressed and density variations with respect to $z$ remain small. Due to the body-forced nature of this flow, quiescent episodes, when the vertical shear instability is suppressed, are associated with a build up of energy in the mean flow as energy is continually and uniformly injected via the forcing. This results in highly intermittent turbulent bursts. The flow is accelerated until apparently a local gradient Richardson number criterion is overcome, KH shear instability is initiated which then grows to finite amplitude, breaks down to disordered motion and then decays until the mean is regenerated once again by the forcing. In fact this is the generic behaviour observed when the forcing itself provides vertical shear, i.e. $\sin(nz)\mathbf{\hat{x}}$ (figure not shown). The final row of figure 2 shows some turbulent remnant in the $T=50$ snapshots with apparent shear instability on the vertical bands of $u.$ Figure 14 shows the full three-dimensional streamwise velocity for four snapshots in time, spanning the burst of spatiotemporal turbulence around $t\simeq 80-90$ apparent in figure 13.

Defining an appropriate gradient Richardson number

[TABLE]

where the subscript denotes averaging in the streamwise $x-$direction, we plot a time series of the volume average $\langle Ri_{G}\rangle_{V}$ and $y-z$ snapshots at various times for the $B=50$ simulation in figure 15. We find that the suppression of shear instability occurs in two distinct classes of thin strips. One class (at a non-trivial angle to the vertical) is aligned with the minimum of the vertical shear. The other class (close to horizontal) is associated with the maximum of the density gradient, as is apparent from the bottom right two panels in figure 2. In other words, the horizontal structures in figure 15 are caused by large values of density gradient in the numerator of $Ri_{G}$ and the inclined structures are caused by minima of the shear in the denominator of $Ri_{G}$. This lattice structure of the two classes of strips effectively covers the entire flow domain, and eventually manages to stabilise the flow globally. Once the flow is accelerated again sufficiently, we then observe instability nucleating in the regions of minimal $Ri_{G}$ in the gaps of this lattice structure, as is apparent in the frames associated with $t=88$ in figures 14 and 15. Overall, for this flow, the volume average of $Ri_{G}$ reaches a maximum value of $Ri_{G}\approx 0.2$ during the stabilising period, but it is not clear whether there is any significance to this particular value.

From a dynamical systems point of view, one may consider these observed cycles to be similar to the type of homoclinic tangle studied to explain bursting in plane Couette flow by van Veen & Kawahara (2011). We appear to observe a close approach to a coherent structure and a return to it via a complex turbulent trajectory which wraps up the stable and unstable manifolds of the underlying solution. Unfortunately, this particular feature of close approach and complex (near) return makes this flow quite challenging computationally. As is particularly apparent for simulation A4 as shown in figure 13, as the buoyancy parameter $B$ is increased further (i.e. the background buoyancy frequency is increased relative to the horizontal shear forcing) the mean flow required to maintain turbulence becomes stronger with more intense small-scale dissipation. This ‘strengthening’ of the turbulence requires both smaller timesteps to maintain numerical stability as the mean flow increases, and also longer integrations to approach steady statistics and average out the increasingly long intermittent bursting time scale.

A.2 Throttling

To attempt to overcome the combined computational challenge of increasingly long times between increasingly intense intermittent bursts as $B$ is increased, we employ a throttling method similar to that used by Chung & Matheou (2012). This throttling modulates the forcing amplitude with the aim of maintaining a mean dissipation rate in the vicinity of a target value. Therefore, we allow the forcing term in equation (5) to have a time varying amplitude, i.e.

[TABLE]

We then adjust $\chi(t)$ based on the instantaneous kinetic energy budget, i.e.

[TABLE]

where $\mathcal{K},$ $\mathcal{I}$, $\mathcal{B}$ and $\mathcal{D}$ are the total kinetic energy density, volume-averaged energy input, buoyancy flux and dissipation rate respectively, as defined in (9)-(10). Following Chung & Matheou (2012) we set a target dissipation rate $\mathcal{D}_{0}$ and seek $\frac{{\rm d}\mathcal{K}}{{\rm d}t}=0,$ i.e. statistically stationary kinetic energy, so that

[TABLE]

In other words when $\mathcal{D}<\mathcal{D}_{0}$ the energy input is increased, so that the total kinetic energy grows in time, and when $\mathcal{D}>\mathcal{D}_{0}$ the flow is decelerated by reducing $\mathcal{K}$. Notice that the connections between the terms in equation (19) are inherently nonlinear and linked to the turbulent cascade. Energy is input at the largest scale and is ultimately dissipated at much smaller scales. For this reason there is an inevitable lag between the adjustment of $\chi$ and the flow response. Furthermore, there are still undoubted computational challenges. To maintain adequate resolution at large external $Re$ we cannot throttle with large $\mathcal{D}_{0}$, although for statistical stationarity, larger values of $B$ require progressively larger $\mathcal{D}_{0}$. For these reasons we choose the parameters in table 1 to maintain at least some semblance of a turbulent state at progressively stronger stratification.

Figure 3 shows equivalent $yz-$plane $\rho$ and $u$ snapshots and means to those shown in figure 2 for the throttled simulations D1, B3, C2 and C3 as listed in table 1, associated with larger values of the buoyancy parameter $B$. Figure 16 shows the time dependence of the kinetic energy density $\mathcal{K}$ and scaled dissipation rate $\mathcal{D}/\mathcal{D}_{0}$ for these throttled simulations. The throttled simulations exhibit approximately statistical stationarity near the target dissipation rate $\mathcal{D}_{0}$ in contrast to the observed time-lagged bursting in both energy density and dissipation rate shown in figure 13 for the unthrottled simulation A4 with $B=50$. Curiously, the throttling method doesn’t reach the target dissipation rate $\mathcal{D}_{0}$ as closely as reported in Chung & Matheou (2012), with a systematic undershoot in the actually occuring dissipation rate $\mathcal{D}$. This difference is presumably due to the different nature of forcing between the two studies. In particular, Chung & Matheou (2012) enforce a specific background mean linear shear rather than a ‘free’ body force as in our simulations. In spite of this discrepancy, the method serves our purpose by providing a much more stationary flow at large $B$ and the opportunity to investigate the trend in layer scale with $B$.