2 VSHF formation in simulations of 2D Boussinesq turbulence

We study VSHF formation using the 2D stochastically excited Boussinesq equations using a unit aspect ratio $(x,z)$ computational domain doubly periodic with length $L$ in the $x$ and $z$ directions. Anticipating the development of horizontal mean structure we use a Reynolds decomposition in which the averaging operator is the horizontal mean. The resulting equations, in terms of the mean velocity, perturbation vorticity, and mean and perturbation buoyancy are

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\overline{u^{\prime }w^{\prime }}-r_{m}U+\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}z^{2}}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\overline{w^{\prime }b^{\prime }}-r_{m}B+\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}B}{\unicode[STIX]{x2202}z^{2}}, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime }}{\unicode[STIX]{x2202}t} & = & \displaystyle -U\frac{\unicode[STIX]{x2202}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime }}{\unicode[STIX]{x2202}x}+w^{\prime }\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}z^{2}}+\frac{\unicode[STIX]{x2202}b^{\prime }}{\unicode[STIX]{x2202}x}-[\mathit{J}(\unicode[STIX]{x1D713}^{\prime },\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime })-\overline{J(\unicode[STIX]{x1D713}^{\prime },\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime })}]\nonumber\\ \displaystyle & & \displaystyle -\,r\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime }+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6E5}^{2}\unicode[STIX]{x1D713}^{\prime }+\sqrt{\unicode[STIX]{x1D700}}S,\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}b^{\prime }}{\unicode[STIX]{x2202}t}=-U\frac{\unicode[STIX]{x2202}b^{\prime }}{\unicode[STIX]{x2202}x}-w^{\prime }\left(N_{0}^{2}+\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}z}\right)-[\mathit{J}(\unicode[STIX]{x1D713}^{\prime },b^{\prime })-\overline{J(\unicode[STIX]{x1D713}^{\prime },b^{\prime })}]-rb^{\prime }+\unicode[STIX]{x1D708}\unicode[STIX]{x0394}b^{\prime }.\end{eqnarray}$$

In these equations an overline indicates a horizontal mean and primes indicate deviations from the horizontal mean so that $f^{\prime }=f-\overline{f}$ . The velocity is $\boldsymbol{u}=(u,w)$ with $u$ and $w$ the horizontal and vertical velocity components, $U=\overline{u}$ is the horizontal mean horizontal velocity, $b$ is the buoyancy with $B=\overline{b}$ the horizontal mean buoyancy, $\unicode[STIX]{x1D713}$ is the streamfunction satisfying $(u,w)=(-\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D713},\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713})$ , and the vorticity is $\unicode[STIX]{x0394}\unicode[STIX]{x1D713}=\unicode[STIX]{x2202}_{x}w-\unicode[STIX]{x2202}_{z}u$ where $\unicode[STIX]{x0394}=\unicode[STIX]{x2202}_{xx}^{2}+\unicode[STIX]{x2202}_{zz}^{2}$ is the Laplacian operator. Perturbation–perturbation advection terms are written using the Jacobian $\mathit{J}(f,g)=(\unicode[STIX]{x2202}_{x}f)(\unicode[STIX]{x2202}_{z}g)-(\unicode[STIX]{x2202}_{x}g)(\unicode[STIX]{x2202}_{z}f)$ . $\sqrt{\unicode[STIX]{x1D700}}S$ denotes the stochastic excitation, which has zero horizontal mean and excites the perturbation vorticity only. $\unicode[STIX]{x1D700}$ controls the strength of the excitation. $N_{0}$ is the constant background buoyancy frequency. Dissipation is provided by Rayleigh drag and diffusion acting on both the buoyancy and vorticity fields. Consistent with previous studies of VSHF formation, dissipation coefficients are assumed equal for vorticity and buoyancy, i.e. the Prandtl numbers associated with the Rayleigh drag and with the diffusive dissipation are each set equal to one. To approximate the effects of diffusive turbulent dissipation, which damps the large scales less rapidly, the Rayleigh drag on the mean fields (with coefficient $r_{m}$ ) is typically taken to be weaker than that on the perturbation fields (with coefficient $r$ ). We refer to (2.1)–(2.4) as the NL equations (for fully nonlinear) to distinguish them from the quasilinear (QL) and S3T systems which we formulate in § 4.

Use of Rayleigh drag in (2.1)–(2.4) departs from the diffusive dissipation commonly used in simulating stratified turbulence (e.g. Smith Reference Smith, Milewski, Smith, Waleffe and Tabak2001). Rayleigh drag provides a simplified parametrization of dissipation that allows the system to reach statistical equilibrium quickly, enabling simulations to obtain the asymptotic state of the VSHF that is difficult to study comprehensively using diffusive dissipation. We emphasize that the essential phenomenon of VSHF formation does not depend crucially on the details of the dissipation, which we demonstrate using examples near the end of the present section.

We choose the stochastic excitation, $\sqrt{\unicode[STIX]{x1D700}}S$ in (2.3), to have the spatial structure of an isotropic ring in Fourier space and to be delta-correlated in time. Figure 1 shows a snapshot of $\sqrt{\unicode[STIX]{x1D700}}S$ (panel a) and its wavenumber power spectrum (panel b), in which $\boldsymbol{k}=(k,m)$ is the vector wavenumber with $k$ and $m$ the horizontal and vertical wavenumber components. The excitation is homogeneous in space and approximately isotropic, with some anisotropy being introduced by the omission of the horizontal mean ( $k=0$ ) and vertical mean ( $m=0$ ) components of the excitation and also by the finite domain size. We set the total wavenumber of the ring, $k_{e}$ , to be global wavenumber six, $k_{e}/(2\unicode[STIX]{x03C0}L)=6$ . As the excitation is delta-correlated in time, the rate at which energy is injected into the flow by the vorticity excitation is a control parameter that is independent of the system state. Here we define the kinetic energy, $K$ , the potential energy, $V$ , and the total energy, $E$ , of the flow as

(2.5a-c ) $$\begin{eqnarray}K=\left[{\textstyle \frac{1}{2}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}\right]\!,\quad V=\left[{\textstyle \frac{1}{2}}N_{0}^{-2}b^{2}\right]\!,\quad E=K+V,\end{eqnarray}$$

in which square brackets indicate the domain average. The energy injection rate as a function of wavenumber, denoted $\unicode[STIX]{x1D700}_{k,m}$ , follows a Gaussian distribution centred at $k_{e}$ so that $\unicode[STIX]{x1D700}_{k,m}=\unicode[STIX]{x1D6FC}\exp [-(|\boldsymbol{k}|-k_{e})^{2}/\unicode[STIX]{x1D6FF}k^{2}]$ , where $\unicode[STIX]{x1D6FF}k=2\unicode[STIX]{x03C0}/L$ sets the ring thickness and $\unicode[STIX]{x1D6FC}$ is a normalization factor chosen so that the total energy injection rate summed over all wavenumbers, $\sum _{k,m}\unicode[STIX]{x1D700}_{k,m}$ , is equal to the value of the parameter $\unicode[STIX]{x1D700}$ appearing in (2.3). With this normalization, $\unicode[STIX]{x1D700}$ corresponds to the rate at which the vorticity excitation injects energy into the system. Global horizontal wavenumbers 1–8 have non-zero excitation and all higher horizontal wavenumbers are omitted from the excitation.

Figure 1. Spatial structure of the stochastic excitation of the vorticity field, $\sqrt{\unicode[STIX]{x1D700}}S$ . (a) A sample realization of the excitation pattern, shown in normalized form as $S(x,z,t)/\text{max}[S(x,z,t)]$ . (b) The wavenumber power spectrum of $S$ , shown in normalized logarithmic form as $\ln (P(k,m)/\text{max}[P(k,m)])$ . Here we define $P(k,m)=\langle |\tilde{S}_{k,m}|^{2}\rangle$ , in which $\tilde{S}_{k,m}(t)$ is the Fourier coefficient of the excitation when $S$ is expanded as $S(x,z,t)=\sum _{k,m}\tilde{S}_{k,m}(t)\text{e}^{\text{i}(kx+mz)}$ . Angle brackets indicate the ensemble average over noise realizations. The parameters of the excitation are $k_{e}/2\unicode[STIX]{x03C0}=6$ and $\unicode[STIX]{x1D6FF}k/2\unicode[STIX]{x03C0}=1$ .

Equations (2.1)–(2.4) are non-dimensionalized by choosing the unit of length to be the domain size, $L$ , and the unit of time to be the Rayleigh damping time of the perturbations, $1/r$ . The non-dimensional parameters of the problem are $\widehat{k_{e}}=Lk_{e}$ , $\widehat{\unicode[STIX]{x1D6FF}k}=L\unicode[STIX]{x1D6FF}k$ , $\widehat{r_{m}}=r_{m}/r$ , $\widehat{\unicode[STIX]{x1D708}}=\unicode[STIX]{x1D708}/(L^{2}r)$ , $\widehat{\unicode[STIX]{x1D700}}=\unicode[STIX]{x1D700}/(r^{3}L^{2})$ , and $\widehat{N_{0}^{2}}=N_{0}^{2}/r^{2}$ . We hold fixed the parameters $\widehat{k_{e}}/2\unicode[STIX]{x03C0}=6$ , $\widehat{\unicode[STIX]{x1D6FF}k}/2\unicode[STIX]{x03C0}=1$ , $\widehat{r_{m}}=0.1$ , $\widehat{\unicode[STIX]{x1D708}}=2.4\times 10^{-5}$ and $\widehat{N_{0}^{2}}=10^{3}$ unless otherwise stated. The choice of $\widehat{k_{e}}$ represents a compromise between providing separation between the excitation scale and the domain scale while minimizing the effects of diffusion on perturbations at the excitation scale. Modelling scale-selective diffusive dissipation motivates setting $\widehat{r_{m}}<1$ and our specific choice to set $\widehat{r_{m}}=0.1$ , so that the mean fields are damped ten times less rapidly than the perturbation fields, is made for computational convenience. We examine the sensitivity of the system to this choice later in figure 6(a). The value of $\widehat{\unicode[STIX]{x1D708}}$ is small and was selected to ensure numerical convergence. The rate of energy injection by the excitation, $\widehat{\unicode[STIX]{x1D700}}$ , is the primary control parameter which is varied to determine the response of the system to changes in excitation.

We choose $\widehat{N_{0}^{2}}=10^{3}$ to place the system in the strongly stratified regime in which VSHFs have previously been found to form (Smith Reference Smith, Milewski, Smith, Waleffe and Tabak2001; Smith & Waleffe Reference Smith and Waleffe2002). The strongly stratified regime is also the regime relevant to the EDJs. For example, taking the equatorial deep stratification as $N_{deep}\sim 2\times 10^{-3}~\text{s}^{-1}$ , a typical gravity wave wavelength of $\unicode[STIX]{x1D706}_{GW}\sim 10~\text{km}$ , and a lateral eddy viscosity of $\unicode[STIX]{x1D708}_{eddy}\sim 100~\text{m}^{2}~\text{s}^{-1}$ gives an effective Rayleigh drag coefficient of $r_{eff}\sim (2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}_{GW})^{2}\unicode[STIX]{x1D708}_{eddy}\sim 4\times 10^{-5}~\text{ s}^{-1}$ and so $\widehat{N_{0}^{2}}_{,EDJ}=N_{deep}^{2}/r_{eff}^{2}\sim 2500$ . Although we do not attempt in this work to model the EDJs, which have 3D structure and are influenced by rotation and boundaries, this estimate suggests that the presently studied idealized turbulence is in the appropriate parameter regime to allow comparison between our VSHF dynamics and EDJ phenomena. For the remainder of the paper we work exclusively in terms of non-dimensional parameters and drop hats in our notation.

We now summarize the behaviour of an NL simulation exhibiting VSHF formation in which the system was integrated from rest over $t\in [0,60]$ with $\unicode[STIX]{x1D700}=0.25$ and the other parameters as described above, which we refer to as the standard parameter case. The standard case value of $\unicode[STIX]{x1D700}$ places the system in the parameter regime in which strong VSHF formation occurs; the sensitivity of the system to $\unicode[STIX]{x1D700}$ is examined in §§ 6, 7 and 9. In § 5 we compare the first- and second-order statistical features of NL simulations with the results of the QL and S3T simulations. To perform the numerical integration we use a 2D finite-difference configuration of DIABLO (Taylor Reference Taylor2008) with 512 gridpoints in both the $x$ and $z$ directions.

To estimate the canonical scales and non-dimensional parameters of the standard case simulation we use the estimates $U_{0}\sim \sqrt{\unicode[STIX]{x1D700}}$ for the velocity scale and $L_{0}\sim 1/k_{e}$ for the length scale. The velocity scale is estimated based on the approximate energy balance in the absence of a VSHF, ${\dot{E}}\approx -2E+\unicode[STIX]{x1D700}$ , together with the estimate $U_{0}\sim \sqrt{E}$ . Using these estimates, the Froude number of the standard parameter case is $Fr\equiv U_{0}/(L_{0}N_{0})\approx 0.6$ , the Ozmidov wavenumber is $k_{O}/(2\unicode[STIX]{x03C0})\approx 56$ , where $k_{O}\equiv (N_{0}^{3}/\unicode[STIX]{x1D700})^{1/2}$ , and the buoyancy wavenumber is $k_{b}/(2\unicode[STIX]{x03C0})\approx 10$ , where $k_{b}\equiv N_{0}/U_{0}\sim N_{0}/\sqrt{\unicode[STIX]{x1D700}}$ . The buoyancy Reynolds number is conventionally defined as $Re_{b}\equiv \unicode[STIX]{x1D700}/(\unicode[STIX]{x1D708}N_{0}^{2})$ and is used to estimate the ratio of the vertical advection term to the viscous damping term in the horizontal momentum equation in 3D stratified turbulence (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). Using this definition, the value of $Re_{b}$ in the standard parameter case is $Re_{b}\approx 10.4$ . Although our system is 2D and includes Rayleigh drag, this estimate of $Re_{b}$ is consistent with the time average value in the standard case simulation of the ratio of interest, $(w^{\prime }\unicode[STIX]{x2202}_{z}u^{\prime })_{RMS}/(-u^{\prime }+\unicode[STIX]{x1D708}\unicode[STIX]{x0394}u^{\prime })_{RMS}\approx 10.7$ , where the time average is calculated over the final 15 time units of the simulation and the subscript RMS denotes the root mean square average over space.

Indicative example snapshots and time series of the NL system are shown in figures 2–4. Near the start of the integration (figure 2 a,b), the structure of the flow reflects the structure of the stochastic excitation and is incoherent with a dominant length scale corresponding to the stochastic excitation scale, $1/k_{e}$ . By $t=60$ (figure 2 c,d) the system has evolved into a state in which the flow is dominated by the VSHF, $U$ , which manifests as horizontal ‘stripes’ in both the vorticity and streamfunction fields with vertical wavenumber $m_{U}/(2\unicode[STIX]{x03C0})=6$ . Simulated realizations of the NL system in the standard parameter case are always found to form a VSHF, but the VSHF wavenumber, $m_{U}$ , differs slightly between simulations when the system is initialized from rest. We focus, in this section, on an example in which $m_{U}/(2\unicode[STIX]{x03C0})=6$ to facilitate comparison with SSD results in § 5. However, VSHFs with $m_{U}/(2\unicode[STIX]{x03C0})=7$ form somewhat more frequently, which we discuss in § 6. We analyse how the VSHF wavenumber, $m_{U}$ , is related to the parameters in §§ 6 and 7, but presently note that in the standard parameter case $m_{U}$ is closely related to the excitation wavenumber, $k_{e}/(2\unicode[STIX]{x03C0})=6$ , and that $m_{U}$ differs from the Ozmidov wavenumber, $k_{O}/(2\unicode[STIX]{x03C0})\approx 56$ , and from the buoyancy wavenumber, $k_{b}/(2\unicode[STIX]{x03C0})\approx 10$ .

The time evolution of $U$ is shown in figure 3(a). The VSHF forms by $t\approx 15$ and persists until the end of the integration. Figure 4 shows the time evolution of the kinetic energy of the VSHF, $\overline{K}$ , and of the perturbations, $K^{\prime }$ , where these energies are defined as

(2.6a,b ) $$\begin{eqnarray}\overline{K}=\left[{\textstyle \frac{1}{2}}U^{2}\right]\!,\quad K^{\prime }=\left[{\textstyle \frac{1}{2}}\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}^{\prime }\right]\!.\end{eqnarray}$$

The VSHF is the energetically dominant feature of the statistically steady flow, containing approximately six times more kinetic energy than the perturbations. In the statistical equilibrium state, the kinetic energy that is injected into the perturbation field by the stochastic excitation is transferred both into the mean flow, thereby maintaining the VSHF, and into the buoyancy field. Energetic balance is maintained by dissipation of the mean and perturbation energies at large scales by Rayleigh drag, with viscosity contributing only weakly to the total dissipation.

Figure 2. Snapshots of the vorticity, streamfunction, and velocity fields for the standard case NL simulation showing the development of the VSHF in turbulence. Just after initialization ( $t=2.5$ ), the vorticity field (a) and the streamfunction and associated velocity field (b) are characterized by perturbations at the scale of the excitation. The system evolves into a statistical equilibrium state by $t=60$ in which the vorticity field (c) is dominated by horizontal stripes with alternating sign indicative of a strong VSHF. The streamfunction and velocity field at $t=60$ (d) show that the VSHF is the dominant feature of the instantaneous flow. Parameters are set to the standard values $r_{m}=0.1$ , $N_{0}^{2}=10^{3}$ , $k_{e}/2\unicode[STIX]{x03C0}=6$ , $\unicode[STIX]{x1D6FF}k/2\unicode[STIX]{x03C0}=1$ , $\unicode[STIX]{x1D708}=2.4\times 10^{-5}$ and $\unicode[STIX]{x1D700}=0.25$ . The buoyancy Reynolds number is $Re_{b}=10.4$ and the Froude number is $Fr=0.6$ .

Figure 3. Development of the VSHF and associated density layers in the standard case NL simulation. (a) Time evolution of the horizontal mean flow, $U$ , which develops from zero at $t=0$ into a persistent VSHF pattern with vertical wavenumber $m_{U}/2\unicode[STIX]{x03C0}=6$ by $t\approx 15$ . (b) Time evolution of the horizontal mean stratification, $\overline{N^{2}}$ , which develops into a pattern with vertical wavenumber $m_{B}/2\unicode[STIX]{x03C0}=12$ that is phase-aligned with $U$ so that regions of weak stratification coincide with the shear regions of the VSHF structure.

Figure 4. Kinetic energy evolution in the standard case NL simulation. In statistically steady state, the kinetic energy of the VSHF (dotted line) is approximately six times that of the perturbations (solid line).

Although the phenomenon of VSHF emergence in stratified turbulence is well known, the concurrent development of coherent horizontal mean structure in the buoyancy field has not been emphasized in the literature. Figure 3(b) shows the time evolution of the horizontal mean stratification $\overline{N^{2}}=N_{0}^{2}+\unicode[STIX]{x2202}_{z}B$ . Although $\overline{N^{2}}$ exhibits more temporal variability than $U$ , it is clear that for these parameter values the turbulent fluxes systematically weaken the stratification ( $\overline{N^{2}}<N_{0}^{2}$ ) in the shear regions of the VSHF. Association of mean stratification anomalies with the mean shear produces a vertical wavenumber in $\overline{N^{2}}$ of $m_{B}/2\unicode[STIX]{x03C0}=12$ , twice that of the $m_{U}/2\unicode[STIX]{x03C0}=6$ structure of the VSHF.

The statistical equilibrium horizontal mean state, obtained by averaging the flow subsequent to a spin-up period of 30 time units, is shown in figure 5. Panels (a,b) show that, for these parameters, the VSHF has a vertical structure that deviates somewhat from harmonic, with flattened shear regions resulting in a profile resembling a sawtooth structure. Comparison of panels (b,c) reveals that the shear extrema coincide with the minima of $\overline{N^{2}}$ . These $\overline{N^{2}}$ minima correspond to narrow density layers in which $\overline{N^{2}}$ is reduced by approximately 40 % relative to $N_{0}^{2}$ . Similar density layers have been reported in observations and simulations of the EDJs (Ménesguen et al. Reference Ménesguen, Hua, Fruman and Schopp2009). As the vertical integral of $\overline{N^{2}}-N_{0}^{2}$ must vanish due to the periodicity of the boundary conditions in the vertical direction by (2.2), the narrow density layers are compensated by regions of enhanced stratification. These regions of enhanced stratification have a characteristic structure in which the $\overline{N^{2}}$ maxima occur just outside the extrema of $U$ , with weak local minima of $\overline{N^{2}}$ occurring at the locations of the VSHF peaks.

The locations of strongest shear and weakest stratification correspond to the local minima of the horizontal mean Richardson number, $\overline{Ri}=\overline{N^{2}}/(\unicode[STIX]{x2202}_{z}U)^{2}$ , as shown in figure 5(d). The minimum value of $\overline{Ri}$ is near $\overline{Ri}\approx 0.8>0.25$ , indicating that the time mean VSHF structure would be free of modal instabilities in the absence of excitation and dissipation by the Miles–Howard (MH) criterion. Although the MH criterion is formally valid only for steady unforced inviscid flows, it remains useful in our stochastically maintained turbulent flow to guide intuition about the maximum stable shear attainable by the VSHF for a given stratification. We note that this usage of the MH criterion differs from an alternative usage in which $Ri$ is used to distinguish between regions of a flow that are likely to become laminar and regions that are likely to maintain turbulence. This alternative interpretation of the implication of $Ri$ is based on the fact that large perturbation growth is obtained by optimal perturbations in shear flows for which $Ri>1/4$ although modal instability is not permitted (Farrell & Ioannou Reference Farrell and Ioannou1993b ). In accord with this result, turbulence is observed to be supported in shear flows with $Ri>1$ (Galperin, Sukoriansky & Anderson Reference Galperin, Sukoriansky and Anderson2007).

Figure 5. Vertical structure of the time average horizontal mean state in the standard case NL simulation. (a) Mean flow, $U$ . (b) Mean shear, $\unicode[STIX]{x2202}U/\unicode[STIX]{x2202}z$ . (c) Mean stratification, $\overline{N^{2}}$ ; the vertical dashed line indicates $N_{0}^{2}$ . (d) Mean Richardson number, $\overline{Ri}=\overline{N^{2}}/(\unicode[STIX]{x2202}U/\unicode[STIX]{x2202}z)^{2}$ ; the vertical dashed line indicates $\overline{Ri}=1/4$ . Profiles are time averages over $t\in [30,60]$ of the structures shown in figure 3.

To demonstrate that VSHF formation is robust to changes in the control parameters, we show in figure 6 the time evolution of $U$ in four additional cases. Panels (a,b) show the response of the system to changes in the dissipation parameters. Panel (a) shows the development of $U$ when the Rayleigh drag on the mean fields is increased by a factor of five ( $r_{m}=0.5$ ). The mean fields in this case are damped half as rapidly as the perturbations, rather than ten times less rapidly as in the standard case. The excitation strength is $\unicode[STIX]{x1D700}=0.5$ and other parameters are as in the standard case, so that $Re_{b}=20.8$ and $Fr=0.84$ . The VSHF has $m_{U}/(2\unicode[STIX]{x03C0})=7$ and is similar to that seen in the standard case in figure 3(a). Panel (b) shows the effect of removing Rayleigh drag entirely ( $r=r_{m}=0$ ), so that all dissipation is provided by diffusion. In this case, some ambiguity arises regarding how the other parameters should be set, as we non-dimensionalize time by the perturbation damping time, $1/r$ , in examples other than this figure. For simplicity we choose to retain all parameters as they are set in the standard case as if Rayleigh drag were still present with $r=1$ , which gives $Re_{b}=10.4$ and $Fr=0.23$ , where for this example only we use the definition $Fr=(\unicode[STIX]{x1D700}k_{e}^{2})^{1/3}/N_{0}$ due to the absence of Rayleigh drag. The VSHF in this example initially emerges with $m_{U}/(2\unicode[STIX]{x03C0})\approx 6$ before transitioning to larger scale (smaller $m_{U}$ ) as the integration is continued. Transition of the VSHF to smaller values of $m_{U}$ for weaker damping or stronger excitation is consistent with previous studies of VSHF emergence (Herring & Métais Reference Herring and Métais1989; Smith Reference Smith, Milewski, Smith, Waleffe and Tabak2001; Smith & Waleffe Reference Smith and Waleffe2002) and is expected on the basis of analysis of the SSD system in the case of strong excitation, as we show in § 7. Figure 6(c,d) show the response of the system to reductions in stratification. In these examples we reduce the excitation strength to $\unicode[STIX]{x1D700}=1.5\times 10^{-2}$ for ease of comparison because VSHFs form more rapidly at these stratification values than they do in the standard case. Panel (c) shows the development of $U$ when the stratification is reduced by a factor of ten relative to the standard case ( $N_{0}^{2}=100$ rather than $N_{0}^{2}=1000$ , corresponding to $Re_{b}=6.3$ , $Fr=0.05$ ) and panel (d) shows the effect of reducing the stratification by a factor of 25 relative to the standard case ( $N_{0}^{2}=40$ , $Re_{b}=15.6$ , $Fr=0.12$ ). As in the case of modified dissipation, the VSHFs in these examples develop with similar structures as in the standard case shown in figure 3(a). We note (not shown) that VSHF formation ceases for sufficiently weak stratification (Smith Reference Smith, Milewski, Smith, Waleffe and Tabak2001; Kumar et al. Reference Kumar, Verma and Sukhatme2017). We return to the dependence of the VSHF on stratification in § 6.

Figure 6. Time evolution of the VSHF in four additional cases. Unless otherwise stated all parameters are as in figure 2. (a) An example with enhanced Rayleigh drag on the mean fields, $r_{m}=0.5$ , and excitation strength $\unicode[STIX]{x1D700}=0.5$ ( $Re_{b}=20.8$ , $Fr=0.8$ ). (b) An example with zero Rayleigh drag on both the mean and perturbations, $r=r_{m}=0$ ( $Re_{b}=10.4$ , $Fr=0.23$ ). Dissipation is provided solely by diffusion. (c,d) Two examples with reduced stratification and with excitation strength $\unicode[STIX]{x1D700}=1.5\times 10^{-2}$ : (c) $N_{0}^{2}=100$ ( $Re_{b}=6.3$ , $Fr=0.05$ ) and (d) $N_{0}^{2}=40$ , ( $Re_{b}=15.6$ , $Fr=0.12$ ). This figure demonstrates that VSHFs form robustly when the dissipation and stratification are varied.

3 Mechanism of horizontal mean structure formation

In a statistically steady state the VSHF, $U$ , and the associated buoyancy structure, $B$ , must be supported against dissipation by perturbation fluxes of momentum and buoyancy as expressed in (2.1)–(2.2). In the absence of any horizontal mean structure (i.e. if $U=B=0$ ), isotropy of the stochastic excitation implies that the statistical mean perturbation momentum flux vanishes ( $\langle \overline{u^{\prime }w^{\prime }}\rangle =0$ , where angle brackets indicate the ensemble average over realizations of the stochastic excitation) and that the statistical mean perturbation buoyancy flux is constant ( $-\unicode[STIX]{x2202}_{z}\langle \overline{w^{\prime }b^{\prime }}\rangle =0$ ). For the observed horizontal mean structures to emerge and persist, their presence must modify the fluxes so that the fluxes reinforce these structures. In this section we analyse the interaction between the turbulence and the horizontal mean state and demonstrate that the horizontal mean structures do influence the turbulent fluxes in this way.

We analyse turbulence–mean state interactions by applying two modifications to (2.1)–(2.4). The first modification is to hold the mean fields constant as $U=U_{test}$ , $B=B_{test}$ . The second modification is to discard the perturbation–perturbation nonlinear terms $[\mathit{J}(\unicode[STIX]{x1D713}^{\prime },\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime })-\overline{\mathit{J}(\unicode[STIX]{x1D713}^{\prime },\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime })}]$ and $[\mathit{J}(\unicode[STIX]{x1D713}^{\prime },b^{\prime })-\overline{\mathit{J}(\unicode[STIX]{x1D713}^{\prime },b^{\prime })}]$ from (2.3)–(2.4). The resulting equations are

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime }}{\unicode[STIX]{x2202}t}=-U_{test}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime }}{\unicode[STIX]{x2202}x}+w^{\prime }\frac{\unicode[STIX]{x2202}^{2}U_{test}}{\unicode[STIX]{x2202}z^{2}}+\frac{\unicode[STIX]{x2202}b^{\prime }}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime }+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6E5}^{2}\unicode[STIX]{x1D713}^{\prime }+\sqrt{\unicode[STIX]{x1D700}}S, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}b^{\prime }}{\unicode[STIX]{x2202}t}=-U_{test}\frac{\unicode[STIX]{x2202}b^{\prime }}{\unicode[STIX]{x2202}x}-w^{\prime }\overline{N^{2}}_{test}-b^{\prime }+\unicode[STIX]{x1D708}\unicode[STIX]{x0394}b^{\prime }, & \displaystyle\end{eqnarray}$$

in which $\overline{N^{2}}_{test}=N_{0}^{2}+\unicode[STIX]{x2202}_{z}B_{test}$ . Equations (3.1)–(3.2) are a system of linear differential equations for the perturbation fields. For this system the time mean fluxes are identical to the ensemble mean fluxes averaged over noise realizations, and either method of averaging can be used to calculate the average fluxes in the presence of the imposed horizontal mean state ( $U=U_{test}$ and $\overline{N^{2}}=\overline{N^{2}}_{test}$ ). We refer to the calculation of perturbation fluxes from (3.1)–(3.2) as test function analysis, as it allows us to probe the turbulent dynamics by imposing chosen test functions for the mean flow and buoyancy, $U_{test}$ and $\overline{N^{2}}_{test}$ . This approach has been applied to estimate perturbation fluxes in the midlatitude atmosphere (Farrell & Ioannou Reference Farrell and Ioannou1993a ) and in wall-bounded shear flows (Farrell & Ioannou Reference Farrell and Ioannou2012; Farrell et al. Reference Farrell, Ioannou and Nikolaidis2017b ), and we will evaluate its effectiveness in the 2D Boussinesq system in § 5.

That the modified perturbation equations are capable of producing realistic perturbation fluxes given the observed mean flow is related to the non-normality of the perturbation dynamics in the presence of shear (Farrell & Ioannou Reference Farrell and Ioannou1996). The modified equations correctly capture the non-normal dynamics, which produce both the positive and negative energetic perturbation–mean flow interactions. The non-normal dynamics of perturbations in stratified shear flow have been analysed in 2D (Farrell & Ioannou Reference Farrell and Ioannou1993b ) and in 3D (Bakas, Ioannou & Kefaliakos Reference Bakas, Ioannou and Kefaliakos2001; Kaminski, Caulfield & Taylor Reference Kaminski, Caulfield and Taylor2014).

As an illustrative example we show in figure 7 the results of test function analysis in the case of an imposed mean state comprising a Gaussian jet peaked in the centre of the domain, $U_{test}=\exp (-50(z-(1/2))^{2})$ , and an unmodified background stratification, $\overline{N^{2}}_{test}=N_{0}^{2}=10^{3}$ . Panel (a) shows the imposed jet, $U_{test}$ , while panel (b) shows the induced perturbation momentum flux divergence, $-\unicode[STIX]{x2202}_{z}\langle \overline{u^{\prime }w^{\prime }}\rangle$ , alongside the negative of the jet dissipation, $(r_{m}-\unicode[STIX]{x1D708}\unicode[STIX]{x2202}_{zz})U_{test}$ . The core of the jet is clearly being supported against dissipation by the perturbation momentum fluxes resulting from its modification of the turbulence. This organization of turbulence producing up-gradient momentum fluxes in the presence of a background shear flow is the essential mechanism of VSHF emergence: an initially perturbative VSHF that arises randomly from turbulent fluctuations modifies the turbulence to produce fluxes reinforcing the initial VSHF. This wave–mean flow mechanism is consistent with the results of rapid distortion theory for stratified shear flow (Galmiche & Hunt Reference Galmiche and Hunt2002) and has been identified in simulations of decaying sheared and stratified turbulence (Galmiche, Thual & Bonneton Reference Galmiche, Thual and Bonneton2002). Wave–mean flow interaction has also been hypothesized to be the mechanism responsible for the formation and maintenance of the EDJs (Muench & Kunze Reference Muench and Kunze1999; Ascani et al. Reference Ascani, Firing, McCreary, Brandt and Greatbatch2015).

Consistent with the results of the NL system shown in § 2, the buoyancy fluxes are also modified by imposing a test function horizontal mean state. Figure 7(c) shows the imposed stratification, $\overline{N^{2}}_{test}$ , which is equal to $N_{0}^{2}$ in this example. Figure 7(d) shows the driving by perturbation fluxes of the stratification anomaly, $-\unicode[STIX]{x2202}_{zz}\langle \overline{w^{\prime }b^{\prime }}\rangle$ , alongside the negative of the dissipation of the stratification anomaly, $(r_{m}-\unicode[STIX]{x1D708}\unicode[STIX]{x2202}_{zz})(\overline{N^{2}}_{test}-N_{0}^{2})$ , which is zero in this Gaussian jet example as $\overline{N^{2}}=N_{0}^{2}$ . The vertical structure of $-\unicode[STIX]{x2202}_{zz}\langle \overline{w^{\prime }b^{\prime }}\rangle$ is complex. For these parameter values the fluxes act to enhance $\overline{N^{2}}$ most strongly at the jet maximum, which departs from the NL results in which $\overline{N^{2}}$ has weak local minima at the locations of the VSHF peaks.

Figure 7. Test function analysis showing the perturbation flux divergences that develop in response to an imposed horizontal mean state consisting of a Gaussian jet and an unmodified background stratification. (a) Imposed jet, $U_{test}$ . (b) The resulting ensemble mean perturbation momentum flux divergence, $-\unicode[STIX]{x2202}_{z}\langle \overline{u^{\prime }w^{\prime }}\rangle$ , and the negative of the dissipation of the jet, $(r_{m}-\unicode[STIX]{x1D708}\unicode[STIX]{x2202}_{zz})U_{test}$ . (c) Imposed stratification, $\overline{N^{2}}_{test}$ , which is equal to $N_{0}^{2}$ in this example. (d) Ensemble mean driving by perturbation fluxes of the stratification anomaly, $-\unicode[STIX]{x2202}_{zz}\langle \overline{w^{\prime }b^{\prime }}\rangle$ , and the negative of the dissipation of the stratification anomaly, $(r_{m}-\unicode[STIX]{x1D708}\unicode[STIX]{x2202}_{zz})(\overline{N^{2}}_{test}-N_{0}^{2})$ , which is zero in this example. This example shows that a Gaussian jet organizes the turbulence so that the perturbation momentum fluxes generally accelerate the jet. The buoyancy fluxes are also organized by the jet in such a way as to drive a stratification anomaly with a complex vertical structure. Parameters are as in figure 2.

This simple example demonstrates the general physical mechanism of horizontal mean structures modifying turbulent fluxes so as to modify the mean state. However, the results of this example indicate that a Gaussian jet together with an unmodified background stratification does not constitute a steady state, as neither the jet acceleration nor the driving of the stratification anomaly due to the perturbation fluxes reflect the specific structure of the imposed mean state ( $U=U_{test}$ and $\overline{N^{2}}=\overline{N^{2}}_{test}$ ). Although the perturbation fluxes generally act to strengthen $U$ , they also distort its structure by sharpening the jet core and driving retrograde jets on the flanks. Similarly, the $\overline{N^{2}}=N_{0}^{2}$ structure is not in equilibrium with the buoyancy fluxes. To maintain a statistically steady mean state as seen in the NL simulations, the turbulence and the mean state must be adjusted by their interaction to produce mean structure which the corresponding fluxes precisely support against dissipation.

Figure 8. Test function analysis showing the perturbation flux divergences that develop in response to an imposed horizontal mean state corresponding to that which emerges in the standard case NL simulation shown in § 2, with $U_{test}$ and $\overline{N^{2}}_{test}$ smoothed and symmetrized. Panels are as in figure 7, with the additional vertical dashed line in panel (c) indicating $N_{0}^{2}$ . This example shows that the horizontal mean structure that emerges in the NL system, consisting of the VSHF and associated density layers, organizes the turbulent fluxes so that these fluxes support the specific structure of the horizontal mean state against dissipation. Parameters are as in figure 2.

To demonstrate how such cooperative equilibria are established, we show in figure 8 the results of test function analysis applied to the case in which $U_{test}$ (panel (a)) and $\overline{N^{2}}_{test}$ (panel (c)) are taken to be the time average profiles from the standard case NL integration discussed in § 2, smoothed and symmetrized so that the sixfold symmetry of the VSHF and twelvefold symmetry of $\overline{N^{2}}$ are made exact. As in the Gaussian jet example, the perturbation momentum fluxes support the jet against dissipation (panel (b)). However, unlike the results obtained in the case of a Gaussian jet, the approximately harmonic VSHF that emerges in the NL system leads to flux divergences that are precisely in phase with $U$ . This provides an explanation for the structure of the emergent VSHF: its approximately harmonic $U$ profile is a structure which the associated statistical equilibrium fluxes precisely support. Similarly, the structure of $\overline{N^{2}}$ is supported against dissipation by the perturbation buoyancy fluxes. Some differences between the structure of the perturbation driving and that of the dissipation are seen in panels (b,d). In particular, the perturbation driving of the jet is slightly too strong, and the perturbation driving of the stratification anomaly is too strongly negative at the local stratification minima that coincide with the VSHF extrema. These differences arise because the VSHF and horizontal mean stratification anomaly tend to strengthen when perturbation–perturbation nonlinearities are discarded (see § 5) and also because the smoothed and symmetrized stratification anomaly has a somewhat weaker local minimum than that which is found in snapshots of the NL system.

This analysis demonstrates that the linear dynamics of the stochastically excited Boussinesq equations produces fluxes consistent with the emergent VSHF and density layers seen in the NL system. In this sense, test function analysis provides a ‘mechanism denial study’ that demonstrates that spectrally local perturbation–perturbation interactions associated with a cascade of energy to large scales are not required to produce VSHFs in stratified turbulence. That VSHF formation does not occur via such a cascade has previously been noted by Smith & Waleffe (Reference Smith and Waleffe2002). The analysis in this section has been conducted using an imposed, constant horizontal mean state. In the next section we extend (3.1)–(3.2) by coupling the dynamics of the mean fields to the linearized perturbation equations to formulate the S3T implementation of SSD for this system.

4 Formulating the QL and S3T equations of motion

The QL system is obtained by combining the perturbation equations (3.1)–(3.2) with the NL equations for the horizontal mean state (2.1)–(2.2). The resulting QL equations of motion are

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\overline{u^{\prime }w^{\prime }}-r_{m}U+\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}z^{2}}, & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\overline{w^{\prime }b^{\prime }}-r_{m}B+\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}B}{\unicode[STIX]{x2202}z^{2}}, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime }}{\unicode[STIX]{x2202}t}=-U\frac{\unicode[STIX]{x2202}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime }}{\unicode[STIX]{x2202}x}+w^{\prime }\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}z^{2}}+\frac{\unicode[STIX]{x2202}b^{\prime }}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime }+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6E5}^{2}\unicode[STIX]{x1D713}^{\prime }+\sqrt{\unicode[STIX]{x1D700}}S, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}b^{\prime }}{\unicode[STIX]{x2202}t}=-U\frac{\unicode[STIX]{x2202}b^{\prime }}{\unicode[STIX]{x2202}x}-w^{\prime }\left(N_{0}^{2}+\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}z}\right)-b^{\prime }+\unicode[STIX]{x1D708}\unicode[STIX]{x0394}b^{\prime }. & \displaystyle\end{eqnarray}$$

This system can also be obtained directly from the NL system (2.1)–(2.4) by discarding the perturbation–perturbation nonlinearities $[\mathit{J}(\unicode[STIX]{x1D713}^{\prime },\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime })-\overline{J(\unicode[STIX]{x1D713}^{\prime },\unicode[STIX]{x0394}\unicode[STIX]{x1D713}^{\prime })}]$ and $[\mathit{J}(\unicode[STIX]{x1D713}^{\prime },b^{\prime })-\overline{\mathit{J}(\unicode[STIX]{x1D713}^{\prime },b^{\prime })}]$ . The QL dynamics is a coupled system that, while simplified, retains the dynamics of the consistent evolution of the horizontal mean state together with the stochastically excited turbulence. The 2D Boussinesq equations in the QL approximation have previously been applied to analyse mean flow formation in the case of an unstable background stratification (Fitzgerald & Farrell Reference Fitzgerald and Farrell2014).

Because (4.3)–(4.4) are linear in perturbation quantities, the QL system does not retain the transfer by perturbation–perturbation interaction of perturbation energy into horizontal wavenumber components that are not stochastically excited. We choose to excite only global horizontal wavenumbers 1–8. The QL system will therefore not exhibit the full range of small-scale motions seen in the NL system. However, in § 5 we compare the results of QL simulations with those of the NL system, and show that the QL system reproduces the large-scale structure formation observed in the NL system. This implies that the small-scale structures produced by perturbation–perturbation interaction in the NL system do not strongly influence the horizontal mean state and that a faithful representation of the turbulence at all scales is inessential for understanding the statistical structure of the turbulence to second order.

The energetics of the QL system, with respect to both the mean and perturbation kinetic and potential energies, is identical to that of the NL system, with the exception that the terms originating from perturbation–perturbation interaction, which redistribute energy within the perturbation field but do not change the domain-averaged kinetic or potential energies, are not retained in the QL system. The QL system thus possesses identical energetics to the NL system in the domain-averaged sense.

Although the QL system constitutes a substantial mathematical and conceptual simplification compared to the NL system, QL dynamics remains stochastic and exhibits significant turbulent fluctuations. These fluctuations obscure the statistical relationships between the horizontal mean structure and the turbulent fluxes discussed in § 3. To understand the mechanism underlying these statistical relationships it is useful to formulate a dynamics directly in terms of statistical quantities, which we refer to as a statistical state dynamics (SSD). We now formulate the S3T dynamics, which is the SSD we use to study our system. S3T is a closure that retains the interactions between the horizontal mean state and the ensemble mean two-point covariance functions of the perturbation fields which determine the turbulent fluxes. For readers unfamiliar with S3T, appendix A provides a derivation of the S3T equations for a reduced model of stratified turbulence illustrating the conceptual utility of this closure in the context of this reduced model.

Derivation of the S3T dynamics begins with the QL equations (4.1)–(4.4). We expand the perturbation fields in horizontal Fourier series as

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}^{\prime }(x,z,t)=\text{Re}\left[\mathop{\sum }_{n=1}^{N_{k}}\tilde{\unicode[STIX]{x1D713}}_{n}(z,t)\text{e}^{\text{i}k_{n}x}\right], & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle b^{\prime }(x,z,t)=\text{Re}\left[\mathop{\sum }_{n=1}^{N_{k}}\tilde{b}_{n}(z,t)\text{e}^{\text{i}k_{n}x}\right]. & \displaystyle\end{eqnarray}$$

Here $N_{k}$ is the number of retained Fourier modes ( $N_{k}=8$ for our choice of stochastic excitation) and $k_{n}=2\unicode[STIX]{x03C0}n$ . Considering the Fourier coefficients as vectors in the discretized numerical system (e.g. $\tilde{\unicode[STIX]{x1D713}}_{n}(z,t)\rightarrow \unicode[STIX]{x1D74D}_{n}(t)$ ), the QL equations (4.3)–(4.4) can be combined into the vector equation

(4.7) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D74D}_{n}\\ \boldsymbol{b}_{n}\end{array}\right)=\unicode[STIX]{x1D63C}_{n}(\boldsymbol{U},\boldsymbol{B})\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D74D}_{n}\\ \boldsymbol{b}_{n}\end{array}\right)+\left(\begin{array}{@{}c@{}}\sqrt{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D743}_{n}\\ 0\end{array}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D743}_{n}=\unicode[STIX]{x1D6E5}_{n}^{-1}\boldsymbol{S}_{n}$ is the $n$ th horizontal Fourier component of the stochastic excitation of the streamfunction. Here $\unicode[STIX]{x1D6E5}_{n}=-k_{n}^{2}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D63F}^{2}$ in which $\unicode[STIX]{x1D644}$ is the identity matrix and $\unicode[STIX]{x1D63F}$ is the discretized vertical derivative operator. The linear dynamical operator $\unicode[STIX]{x1D63C}_{n}$ is given by the expression

(4.8) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D63C}_{n}(\boldsymbol{U},\boldsymbol{B})\nonumber\\ \displaystyle & & \displaystyle \quad =\left(\begin{array}{@{}cc@{}}-\text{i}k_{n}\unicode[STIX]{x1D6E5}_{n}^{-1}\text{diag}(\boldsymbol{U})\unicode[STIX]{x1D6E5}_{n}+\text{i}k_{n}\unicode[STIX]{x1D6E5}_{n}^{-1}\text{diag}(\unicode[STIX]{x1D63F}^{2}\boldsymbol{U})-\unicode[STIX]{x1D644}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6E5}_{n} & \text{i}k_{n}\unicode[STIX]{x1D6E5}_{n}^{-1}\\ -\text{i}k_{n}N_{0}^{2}\unicode[STIX]{x1D644}-\text{i}k_{n}\text{diag}(\unicode[STIX]{x1D63F}\boldsymbol{B}) & -\text{i}k_{n}\text{diag}(\boldsymbol{U})-\unicode[STIX]{x1D644}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6E5}_{n}\end{array}\right)\!,\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

in which diag $(\boldsymbol{v})$ denotes the diagonal matrix for which the non-zero elements are given by the entries of the column vector  $\boldsymbol{v}$ .

We now make use of the ergodic assumption that horizontal averages and ensemble averages are equal, so that, for example, $U=\overline{u}=\langle u\rangle$ and $\overline{u^{\prime }w^{\prime }}=\langle u^{\prime }w^{\prime }\rangle$ . For our system, which is statistically horizontally uniform, this assumption is justified in a domain large enough so that several approximately independent perturbation structures are found at each height – as seen, for example, in figure 2(b). It can then be shown (using the fact that $\sqrt{\unicode[STIX]{x1D700}}S$ is delta-correlated in time) that the ensemble mean covariance matrix, defined as

(4.9) $$\begin{eqnarray}\unicode[STIX]{x1D63E}_{n}=\left\langle \left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D74D}_{n}\\ \boldsymbol{b}_{n}\end{array}\right)\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D74D}_{n}^{\dagger } & \boldsymbol{b}_{n}^{\dagger }\end{array}\right)\right\rangle =\left(\begin{array}{@{}cc@{}}\langle \unicode[STIX]{x1D74D}_{n}\unicode[STIX]{x1D74D}_{n}^{\dagger }\rangle & \langle \unicode[STIX]{x1D74D}_{n}\boldsymbol{ b}_{n}^{\dagger }\rangle \\ \langle \boldsymbol{b}_{n}\unicode[STIX]{x1D74D}_{n}^{\dagger }\rangle & \langle \boldsymbol{b}_{n}\boldsymbol{b}_{n}^{\dagger }\rangle \end{array}\right)=\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D63E}_{\unicode[STIX]{x1D713}\unicode[STIX]{x1D713},n} & \unicode[STIX]{x1D63E}_{\unicode[STIX]{x1D713}b,n}\\ \unicode[STIX]{x1D63E}_{\unicode[STIX]{x1D713}b,n}^{\dagger } & \unicode[STIX]{x1D63E}_{bb,n}\end{array}\right),\end{eqnarray}$$

in which daggers indicate Hermitian conjugation, evolves according to the time-dependent Lyapunov equation

(4.10) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D63E}_{n}=\unicode[STIX]{x1D63C}_{n}(\boldsymbol{U},\boldsymbol{B})\unicode[STIX]{x1D63E}_{n}+\unicode[STIX]{x1D63E}_{n}\unicode[STIX]{x1D63C}_{n}(\boldsymbol{U},\boldsymbol{B})^{\dagger }+\unicode[STIX]{x1D700}\unicode[STIX]{x1D64C}_{n}, & \displaystyle\end{eqnarray}$$

(4.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64C}_{n}=\left[\begin{array}{@{}cc@{}}\langle \unicode[STIX]{x1D743}_{n}\unicode[STIX]{x1D743}_{n}^{\dagger }\rangle & 0\\ 0 & 0\end{array}\right], & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D64C}_{n}$ is the ensemble mean covariance matrix of the stochastic excitation and has non-zero entries only in the upper-left block matrix because we apply excitation only to the vorticity field. Equation (4.10) constitutes the perturbation dynamics of the S3T system and is the S3T analogue of the QL equations (4.3)–(4.4).

To complete the derivation of the S3T system it remains to write the mean equations (4.1)–(4.2) in terms of the covariance matrix. The ensemble mean perturbation flux divergences can be written as functions of the covariance matrix as

(4.12) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\langle \overline{u^{\prime }w^{\prime }}\rangle =\mathop{\sum }_{n=1}^{N_{k}}\frac{k_{n}}{2}\text{Im}[\text{vecd}(\unicode[STIX]{x1D6E5}_{n}\unicode[STIX]{x1D63E}_{\unicode[STIX]{x1D713}\unicode[STIX]{x1D713},n})], & \displaystyle\end{eqnarray}$$

(4.13) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\langle \overline{w^{\prime }b^{\prime }}\rangle =\mathop{\sum }_{n=1}^{N_{k}}\frac{k_{n}}{2}\text{Im}[\text{vecd}(\unicode[STIX]{x1D63F}\unicode[STIX]{x1D63E}_{\unicode[STIX]{x1D713}b,n})], & \displaystyle\end{eqnarray}$$

in which vecd $(\unicode[STIX]{x1D648})$ denotes the vector comprising the diagonal elements of the matrix $\unicode[STIX]{x1D648}$ . The mean state dynamics then become

(4.14) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\boldsymbol{U}=\mathop{\sum }_{n=1}^{N_{k}}\frac{k_{n}}{2}\text{Im}[\text{vecd}(\unicode[STIX]{x1D6E5}_{n}\unicode[STIX]{x1D63E}_{\unicode[STIX]{x1D713}\unicode[STIX]{x1D713},n})]-r_{m}\boldsymbol{U}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D63F}^{2}\boldsymbol{U}, & \displaystyle\end{eqnarray}$$

(4.15) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\boldsymbol{B}=\mathop{\sum }_{n=1}^{N_{k}}\frac{k_{n}}{2}\text{Im}[\text{vecd}(\unicode[STIX]{x1D63F}\unicode[STIX]{x1D63E}_{\unicode[STIX]{x1D713}b,n})]-r_{m}\boldsymbol{B}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D63F}^{2}\boldsymbol{B}. & \displaystyle\end{eqnarray}$$

Equations (4.10), (4.14) and (4.15) together constitute the S3T SSD closed at second order.

The S3T system is deterministic and autonomous and so provides an analytic description of the evolving relationships between the statistical quantities of the turbulence up to second order, including fluxes and horizontal mean structures, without the turbulent fluctuations inherent in the dynamics of particular realizations of turbulence, such as those present in the NL and QL systems. Although some previous attempts to formulate turbulence closures have been found to have inconsistent energetics (Kraichnan Reference Kraichnan1957; Ogura Reference Ogura1963), the dynamics of the S3T system are QL, and so S3T inherits the consistent energetics of the QL system. That the second-order S3T closure is capable of capturing the dynamics of VSHF formation, as we demonstrate in § 5, is due to the appropriate choice of mean. In the present formulation of S3T, we have chosen the mean to be the horizontal average structure so that the second-order closure retains the QL mean–perturbation interactions that account for the mechanism of VSHF formation in turbulence.

Before proceeding to analysis of the QL and S3T systems, we wish to make two remarks regarding the mathematical structure and physical basis of S3T. First, we note that the ergodic assumption used in deriving S3T is formally valid in the limit that the horizontal extent of the domain tends to infinity and the number of independent perturbation structures at each height correspondingly tends to infinity. In this ideal limit described by S3T, the statistical homogeneity of the turbulence is only broken by the initial state of the horizontal mean structure (which in the examples is perturbatively small), which then determines the phase of the emergent VSHF in the vertical direction. In simulations of the QL and NL systems in a finite domain, the initial mean structure instead results from random Reynolds stresses arising from fluctuations in the perturbation fields. Second, we note that S3T is a canonical closure of the turbulence problem at second order, in that it is a truncation of the cumulant expansion at second order achieved by setting the third cumulant to zero. The mathematical structure of the cumulant expansion determines which nonlinearities are retained and discarded in the QL system, which is a stochastic approximation to the ideal S3T closure. Wave–mean flow coupling enters the equations through second-order cumulants and so is retained, while perturbation–perturbation nonlinearities enter as third-order cumulants and so are not retained.

In the next section we demonstrate that the QL and S3T systems reproduce the major statistical phenomena observed in the NL system.

5 Comparison of the NL, QL and S3T systems

The most striking feature of the standard case NL simulation discussed in § 2 is the spontaneous development of a VSHF, $U$ , with $\text{max}(U)\approx 1$ and vertical wavenumber $m_{U}/2\unicode[STIX]{x03C0}=6$ . The horizontal mean stratification, $\overline{N^{2}}$ , is also modified by the turbulence and develops a structure with vertical wavenumber $m_{B}/2\unicode[STIX]{x03C0}=12$ in phase with $U$ such that weakly stratified density layers develop in the regions of strongest mean shear. The VSHF in the NL system is approximately steady in time, while the horizontal mean stratification is more variable. In this section we compare these NL results to the behaviour of the QL and S3T systems for the same parameter values. We initialize the QL system from rest, matching the procedure used for the NL system. We initialize the S3T system with $\unicode[STIX]{x1D63E}_{n}$ corresponding to homogeneous turbulence together with a small VSHF perturbation (amplitude 0.1) with $m_{U}/2\unicode[STIX]{x03C0}=6$ that is slightly modified by additional small perturbations (amplitude 0.005). We note that the details of the S3T initialization are unimportant in this example because, as we will show in § 6, the VSHF emerges via a linear instability of the homogeneous turbulence and so any sufficiently small initial perturbation to the S3T system will evolve into a $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF for these parameter values.

Figure 9. Development of the VSHF and associated density layers in the QL and S3T systems. Panels show the time evolution of (a) $U$ in the QL system, (b) $\overline{N^{2}}$ in the QL system, (c) $U$ in the S3T system, and (d) $\overline{N^{2}}$ in the S3T system. This figure demonstrates that the QL and S3T systems reproduce the phenomenon of spontaneous VSHF and density layer formation shown in figure 3 for the NL system. Parameters are as in figure 2.

Figure 10. Comparison of the kinetic energy evolution and equilibrium VSHF profiles in the NL, QL and S3T systems. (a) Mean and perturbation kinetic energy evolution. (b) Aligned VSHF profiles. The NL and QL profiles are averaged over $t\in [30,60]$ and the S3T profile is taken to be the state after the S3T system reaches a fixed point. This figure demonstrates that VSHF emergence in the S3T and QL systems occurs with similar structure and energy evolution to that which occurs in the NL system. Parameters are as in figure 2.

Figure 9(a,c) shows the time evolution of the VSHF in the QL and S3T systems (see figure 3 for the corresponding evolution in the NL system). The QL and S3T systems develop VSHF structures with $m_{U}/2\unicode[STIX]{x03C0}=6$ and the $U$ profiles in the NL, QL and S3T systems are compared in figure 10(b). For the NL and QL systems the profiles are time-averaged over $t\in [30,60]$ , while for the S3T system we show the $U$ state after the S3T system has reached a fixed point. The aligned VSHF structures agree well across the three systems. The time evolution of the horizontal mean stratification, $\overline{N^{2}}$ , in the QL and S3T systems is shown in figure 9(b,d). Like the NL stratification, the QL profile of $\overline{N^{2}}$ develops a $m_{B}/2\unicode[STIX]{x03C0}=12$ structure that is more variable in time than $U$ and is phase-aligned with $U$ so that $\overline{N^{2}}$ is weakest in the regions of strongest shear. The S3T system behaves similarly but is free of fluctuations. The evolution of $\overline{N^{2}}$ in the S3T system also reveals that the vertical structure of $\overline{N^{2}}$ changes over time. During the development of the VSHF ( $t\lesssim 8$ ), the stratification is enhanced in the regions of strongest shear. As the VSHF begins to equilibrate at finite amplitude, the $\overline{N^{2}}$ profile reorganizes such that the shear regions are the most weakly stratified. Such reorganization may also occur in the NL and QL systems but is difficult to identify due to the fluctuations present in these systems.

Figure 10(a) shows the evolution of the mean and perturbation kinetic energies of the NL, QL and S3T systems. The growth rate of mean kinetic energy is similar in these three systems. The equilibrium mean energies differ somewhat among the systems, with the VSHFs in the S3T and QL systems having more energy than the NL VSHF. The relative weakness of the NL VSHF is consistent with the scattering of perturbation energy to small scales by the perturbation–perturbation advection terms that are included in NL but not in QL or S3T. The temporal variability of the NL and QL VSHFs, as indicated by the fluctuations in $\overline{K}$ , is similar in the stochastic NL and QL systems. The VSHF in S3T is time-independent once equilibrium has been reached, as the S3T VSHF corresponds to a fixed point of the S3T dynamics.

Figure 11. Vertical structure of the horizontal mean states of the QL and S3T systems. Panels are as in figure 5, with solid lines showing the S3T state and dotted lines showing the QL state. This figure demonstrates that the QL and S3T systems capture the structure of the horizontal mean state in the NL system, including the phase relationship between $U$ and $\overline{N^{2}}$ . Parameters are as in figure 2.

The relationship between the $U$ and $\overline{N^{2}}$ structures is shown in figure 11 for the QL (dotted curves) and S3T (solid curves) systems (see figure 5 for the corresponding structures in the NL system). The equilibrium horizontal mean structures in the QL and S3T systems agree well with those of the NL system. The $U$ profiles (panel (a)) are approximately harmonic with somewhat flattened shear regions and, remarkably, the detailed structure of $\overline{N^{2}}$ seen in the NL integration is reproduced by the QL and S3T systems (panel (c)), which discard perturbation–perturbation nonlinear interactions. In particular, the presence of weak local stratification minima at the locations of the VSHF peaks is captured by the QL and S3T systems.

Figure 12. Comparison of the wavenumber power spectra of kinetic ( $K$ ) and potential ( $V$ ) energy in the NL, QL and S3T systems. (a–c) 2D $K$ spectra of the (a) NL, (b) QL and (c) S3T systems as functions of $(k,m)$ . (d–f) 2D $V$ spectra of those corresponding systems. 2D spectra are shown in terms of their natural logarithms and no normalization is performed. (g,h) Kinetic energy spectra in the conventional 1D form as functions of (g) vertical wavenumber, $m$ , and (h) horizontal wavenumber, $k$ . In panel (g), the contributions to the spectra from the VSHFs in each system are also shown. This figure demonstrates that the QL and S3T systems reproduce structural details of the turbulence beyond the horizontal mean state, including the wavenumber distribution of perturbation energy at large scales. Parameters are as in figure 2.

The above comparisons demonstrate that the horizontal mean structures and domain mean kinetic energies of the QL and S3T systems show good agreement with those of the NL system. Figure 12 compares the energy spectra in the three systems. In panels (a–f), the 2D spectra of kinetic and potential energy are compared as functions of $(k,m)$ , while in panels (g,h) the kinetic energy spectra are compared in the more traditional 1D integrated forms as functions of $k$ and $m$ separately. Panel (a) shows the kinetic energy spectrum of the NL system. The dominant and most important feature of the $K$ spectrum is the concentration of energy at $(k,m)=2\unicode[STIX]{x03C0}(0,6)$ , which corresponds to the $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF structure. This feature is also evident in panel (g), which shows that the peak in the vertical wavenumber spectrum of NL kinetic energy is dominated by the VSHF component of the flow. The energy of the VSHF is also spread across the neighbouring vertical wavenumbers, reflecting both the deviation of the structure of the VSHF from a pure harmonic and also that fluctuations in the VSHF structure project onto nearby vertical wavenumbers. Away from the $k=0$ axis, the $K$ spectrum reveals the expected concentration of energy on the ring of excited wavenumbers $k^{2}+m^{2}=k_{e}^{2}$ , and the spread of this ring to higher $m$ values. This spread is due to the shearing of the ring by the $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF, which produces the sum and difference wavenumber components. The quantitative structure of the spectrum associated with this spreading can be seen more clearly in panel (g), in which selected power law slopes are provided for reference.

The most important features of the 2D $K$ spectrum of the NL system are captured by the QL and S3T systems. The QL $K$ spectrum (figure 12 b) reproduces the primary feature of the energetic dominance of the VSHF over the perturbation field, as well as the minor features of concentration of energy at $k_{e}$ and the spread of the excited ring structure to higher vertical wavenumbers. The 1D spectrum in panel (g) shows that the QL system quantitatively captures the spectrum of perturbation kinetic energy in the wavenumber range $m/(2\unicode[STIX]{x03C0})\lesssim 80$ . We note that the stochastic excitation directly influences the energy spectrum only near $m/(2\unicode[STIX]{x03C0})\approx 6$ , and so the agreement seen in panel (g) is not a direct result of the structure of the excitation. The primary difference between the $K$ spectra of the NL and QL systems is that the NL system scatters some kinetic energy into the unexcited part of the horizontal wavenumber spectrum ( $|k|/2\unicode[STIX]{x03C0}>8$ ), whereas these unexcited wavenumber components have no energy in the QL system, as can also be seen in panel (h). The vertical wavenumber spectrum of $K$ in the NL system, shown in panel (g), also contains small-scale structure for $m\gtrsim 80$ that is not present in the QL system and so can be attributed to perturbation–perturbation nonlinearity. The S3T $K$ spectrum (figure 12 c) also captures the most important features of the NL spectrum, but some differences between the S3T spectrum and those of the NL and QL systems are also visible. In the S3T system the VSHF energy is more strongly concentrated in the $m_{U}/2\unicode[STIX]{x03C0}=6$ harmonic than it is in the NL and QL systems. Additionally, the concentration of energy at the excited ring and the spread of energy to higher $m$ are more distinct in the S3T system than in the NL and QL systems, in which the gaps are filled in by a broad background spectrum. These features are also visible in panel (g). These minor differences between the spectrum of S3T and those of the NL and QL systems are due in part to the absence of fluctuations in the S3T system that are present in the stochastic NL and QL systems. Noise in the stochastic systems produces VSHF fluctuations that spread mean flow energy into $k=0$ modes neighbouring the $m_{U}/2\unicode[STIX]{x03C0}=6$ harmonic. These transient VSHF fluctuations also contribute to producing the broad background spectrum seen in the NL and QL systems by shearing the ring of excited wavenumbers.

The spectrum of potential energy in the NL system is shown in figure 12(d). Unlike the $K$ spectrum, which is dominated by the horizontal mean flow, $U$ , the $V$ spectrum is not dominated by the horizontal mean buoyancy, $B$ , although a peak is evident at the $m_{B}/2\unicode[STIX]{x03C0}=12$ component. In this sense, the VSHF is a ‘manifest’ structure, whereas the horizontal mean density layers are ‘latent’ structures (Berloff, Kamenkovich & Pedlosky Reference Berloff, Kamenkovich and Pedlosky2009). Other features of the spectrum are the expected concentration of potential energy at the ring wavenumber and the spread of the ring to higher vertical wavenumbers as was found for the $K$ spectrum. The $V$ spectra for the QL and S3T systems are shown in figure 12(e,f). The QL and S3T spectra capture the peak associated with the horizontal mean buoyancy layers, the concentration of potential energy at the excitation scale, and the spread of the ring to higher $m$ . The differences between the three $V$ spectra are similar to those identified when comparing the three $K$ spectra.

Agreement between the NL, QL and S3T systems indicates that the QL dynamics of the horizontal mean state interacting with the perturbation field accounts for the physical mechanisms responsible for determining the most important aspects of the energy spectra. We emphasize that the QL and S3T systems involve no free parameters, and that the demonstrated agreement between the three systems is not the result of parameter tuning. Because the QL and S3T systems do not contain the perturbation–perturbation interactions required to produce a turbulent cascade, agreement between these systems and the NL system indicates that, in the present model configuration, such a cascade is not essentially involved in determining the equilibrium turbulent state, including the large-scale spectrum. Nonlinear cascades in the NL system only weakly influence the large-scale dynamics and energetics, and we note that for this reason the parameter $\unicode[STIX]{x1D700}$ , which is the energy injection rate of the stochastic excitation, should not be interpreted as a rate of turbulent dissipation in the sense of a classical cascade. In the NL system we employ to model VSHF formation the turbulent dissipation rate in the classical cascade sense is close to zero and it is exactly zero in the QL and S3T systems.

Similar results have been obtained for barotropic turbulence characterized by strong zonal jets. In the presence of such jets the meridional wavenumber ( $\ell$ ) spectrum of the zonal flow obtains a well-known $\ell ^{-5}$ structure at large scales (Huang, Galperin & Sukoriansky Reference Huang, Galperin and Sukoriansky2001; Galperin, Sukoriansky & Dikovskaya Reference Galperin, Sukoriansky and Dikovskaya2010). Constantinou (Reference Constantinou2015) showed that the barotropic S3T system, with underlying QL dynamics, captures this $\ell ^{-5}$ spectrum, indicating that the large-scale structure of the spectrum in barotropic turbulence with strong zonal jets is primarily determined by perturbation–mean interaction rather than by a turbulent cascade.

QL dynamics does not account for the spectrum at very small scales, which is produced by perturbation–perturbation nonlinearity and is inessential to the dynamics of VSHFs. We note that these small-scale features may be important to the stirring of passive tracers at small scales (Sukoriansky, Dikovskaya & Galperin Reference Sukoriansky, Dikovskaya and Galperin2009), which the QL and S3T systems would not be expected to accurately capture in the present model turbulence.

Motivated by these results we proceed in the rest of this paper to exploit the S3T system to analyse the mechanisms underlying the organization of structure in stratified turbulence.

6 Linear stability analysis of the S3T system

In the previous section we showed that the S3T system reproduces the essential statistical features, up to second order, of the NL system, including both the structure of the horizontal mean state as well as the spectral characteristics of the perturbation field. The S3T system can be understood and analysed with much greater clarity than the NL system because the S3T system is a deterministic and autonomous dynamical system and is amenable to the usual techniques of dynamical systems analysis. In this section we show that the emergence of VSHFs in 2D stratified turbulence can be traced to a linear instability in the SSD of the stationary state of homogeneous turbulence that has analytic expression in the S3T SSD while lacking analytic expression in the dynamics of single realizations. To determine the properties of this instability, and in particular to understand how the vertical scale of the initially emergent VSHF is selected, we now perform a linear stability analysis of the S3T system.

Before linearizing the S3T system it is first necessary to obtain the fixed point statistical state that is unstable to VSHF formation. As shown in § 5, the equilibrium state with a finite-amplitude VSHF and modified horizontal mean stratification is a fixed point of the S3T system. However, the fixed point of the S3T system whose stability we wish to analyse is the state of homogeneous turbulence that is excited by the stochastic excitation and equilibrated by dissipation. This homogeneous state is obscured in the NL and QL systems, both by noise fluctuations and (in examples for which it is SSD unstable) by the development of a VSHF, but roughly corresponds to the interval of nearly constant perturbation kinetic energy at early times ( $t\lesssim 5$ ) in figure 10(a). If homogeneous turbulence is unstable we obtain an explanation for the observed VSHF formation, since the alternative possibility of sustained homogeneous turbulence is not possible in the presence of small perturbations.

For homogeneous turbulence $\boldsymbol{U}=\boldsymbol{B}=0$ and from (4.10) the steady-state perturbation covariance matrix at wavenumber $k_{n}$ obeys

(6.1) $$\begin{eqnarray}\unicode[STIX]{x1D63C}_{n}^{\star }\unicode[STIX]{x1D63E}_{n}^{\star }+\unicode[STIX]{x1D63E}_{n}^{\star }\unicode[STIX]{x1D63C}_{n}^{\star ,\dagger }+\unicode[STIX]{x1D700}\unicode[STIX]{x1D64C}_{n}=0,\end{eqnarray}$$

where the $\unicode[STIX]{x1D63C}_{n}^{\star }$ operator is given by

(6.2) $$\begin{eqnarray}\unicode[STIX]{x1D63C}_{n}^{\star }=\left(\begin{array}{@{}cc@{}}-\unicode[STIX]{x1D644}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6E5}_{n} & \text{i}k_{n}\unicode[STIX]{x1D6E5}_{n}^{-1}\\ -\text{i}k_{n}N_{0}^{2}\unicode[STIX]{x1D644} & -\unicode[STIX]{x1D644}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6E5}_{n}\end{array}\right).\end{eqnarray}$$

Equation (6.1) can be solved analytically for $\unicode[STIX]{x1D63E}_{n}^{\star }$ , and we show details of the solution in appendix B. Figure 13 shows the kinetic and potential energy spectra for this fixed point homogeneous turbulent state.

Figure 13. Spectral structure of the homogeneous S3T fixed point. (a) Kinetic energy ( $K$ ) spectrum. (b) Potential energy ( $V$ ) spectrum. The spectra are shown in terms of their natural logarithms and no normalization is performed. The $K$ and $V$ spectra are nearly identical to one another, even though only the vorticity field is stochastically excited, due to the strong stratification. This figure shows that the homogeneous turbulence from which the VSHF emerges inherits its structure directly from the stochastic excitation whose structure is shown in figure 1. Parameters are as in figure 2.

To analyse the linear stability of this homogeneous turbulent state we perturb the S3T state, $(\unicode[STIX]{x1D63E}_{n},\boldsymbol{U},\boldsymbol{B})$ , about the fixed point, $(\unicode[STIX]{x1D63E}_{n}^{\star },0,0)$ , as

(6.3a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D63E}_{n}=\unicode[STIX]{x1D63E}_{n}^{\star }+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63E}_{n},\quad \boldsymbol{U}=\unicode[STIX]{x1D6FF}\boldsymbol{U},\quad \boldsymbol{B}=\unicode[STIX]{x1D6FF}\boldsymbol{B},\end{eqnarray}$$

where the $\unicode[STIX]{x1D6FF}$ notation indicates that the first-order terms are treated as infinitesimal perturbations. The operator $\unicode[STIX]{x1D63C}_{n}$ in (4.8) may then be written as

(6.4) $$\begin{eqnarray}\unicode[STIX]{x1D63C}_{n}=\unicode[STIX]{x1D63C}_{n}^{\star }+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63C}_{n},\end{eqnarray}$$

where $\unicode[STIX]{x1D63C}_{n}^{\star }$ is given by (6.2) and

(6.5) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63C}_{n}=\left(\begin{array}{@{}cc@{}}-\text{i}k_{n}\unicode[STIX]{x1D6E5}_{n}^{-1}\text{diag}(\unicode[STIX]{x1D6FF}\boldsymbol{U})\unicode[STIX]{x1D6E5}_{n}+\text{i}k_{n}\unicode[STIX]{x1D6E5}_{n}^{-1}\text{diag}(\unicode[STIX]{x1D63F}^{2}\unicode[STIX]{x1D6FF}\boldsymbol{U}) & 0\\ -\text{i}k_{n}\text{diag}(\unicode[STIX]{x1D63F}\unicode[STIX]{x1D6FF}\boldsymbol{B}) & -\text{i}k_{n}\text{diag}(\unicode[STIX]{x1D6FF}\boldsymbol{U})\end{array}\right).\end{eqnarray}$$

The linearized equations of motion are

(6.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D6FF}\boldsymbol{U}=\mathop{\sum }_{n=1}^{N_{k}}\frac{k_{n}}{2}\text{Im}[\text{vecd}(\unicode[STIX]{x1D6E5}_{n}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63E}_{\unicode[STIX]{x1D713}\unicode[STIX]{x1D713},n})]-r_{m}\unicode[STIX]{x1D6FF}\boldsymbol{U}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D63F}^{2}\unicode[STIX]{x1D6FF}\boldsymbol{U}, & \displaystyle\end{eqnarray}$$

(6.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D6FF}\boldsymbol{B}=\mathop{\sum }_{n=1}^{N_{k}}\frac{k_{n}}{2}\text{Im}[\text{vecd}(\unicode[STIX]{x1D63F}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63E}_{\unicode[STIX]{x1D713}b,n})]-r_{m}\unicode[STIX]{x1D6FF}\boldsymbol{B}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D63F}^{2}\unicode[STIX]{x1D6FF}\boldsymbol{B}, & \displaystyle\end{eqnarray}$$

(6.8) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63E}_{n}=\unicode[STIX]{x1D63C}_{n}^{\star }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63E}_{n}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63E}_{n}\unicode[STIX]{x1D63C}_{n}^{\star ,\dagger }+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63C}_{n}\unicode[STIX]{x1D63E}_{n}^{\star }+\unicode[STIX]{x1D63E}_{n}^{\star }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63C}_{n}^{\dagger }. & \displaystyle\end{eqnarray}$$

As usual in linear stability analysis, we express the solutions of (6.6)–(6.8) in terms of the eigenvectors and eigenvalues of the system. The natural matrix form of the S3T equations obscures the operator-vector structure of the linearized system. The most direct technique for conducting the eigenanalysis is to rewrite the equations in superoperator form by unfolding the matrices $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D63E}_{n}$ (Farrell & Ioannou Reference Farrell and Ioannou2002). This technique results in linear operators of very high dimension for which eigenanalysis is expensive. We use an alternative method to obtain the eigenstructures in which the linearized equations are rewritten as coupled Sylvester equations (see appendix B in Constantinou et al. (Reference Constantinou, Farrell and Ioannou2014)).

We note that, for our choice of stochastic excitation, equations (6.6)–(6.8) decouple into two separate eigenproblems: one determining the eigenmodes involving mean flow perturbations $\unicode[STIX]{x1D6FF}\boldsymbol{U}$ , which have $\unicode[STIX]{x1D6FF}\boldsymbol{B}=0$ , and a separate eigenproblem determining the eigenmodes involving mean buoyancy perturbations $\unicode[STIX]{x1D6FF}\boldsymbol{B}$ , which have $\unicode[STIX]{x1D6FF}\boldsymbol{U}=0$ . The eigenproblem involving $\unicode[STIX]{x1D6FF}\boldsymbol{U}$ gives unstable eigenmodes associated with growing VSHFs for the parameter regime we address in this work, while the mean buoyancy eigenproblem has only stable eigenmodes in our parameter regime. The mean buoyancy eigenproblem is therefore irrelevant, in our parameter range, to VSHF formation and we focus on the eigenproblem concerning $\unicode[STIX]{x1D6FF}\boldsymbol{U}$ .

We now describe the results of the eigenanalysis of (6.6)–(6.8). As the fixed point underlying the linearization corresponds to homogeneous turbulence, the eigenfunctions have harmonic structure in $z$ so that $\unicode[STIX]{x1D6FF}\boldsymbol{U}$ and $\unicode[STIX]{x1D6FF}\langle \overline{u^{\prime }w^{\prime }}\rangle$ are both proportional to $\text{e}^{st}\text{e}^{\text{i}m_{U}z}$ . For each $m_{U}$ permitted by the periodic domain there is a dominant eigenmode with eigenvalue $s(m_{U})$ . For the parameter range we study, we find that these eigenvalues are real, corresponding to structures for which the perturbation momentum flux divergence, $-\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D6FF}\langle \overline{u^{\prime }w^{\prime }}\rangle$ , and the mean flow, $\unicode[STIX]{x1D6FF}\boldsymbol{U}$ , are aligned in phase.

Figure 14. Growth rates of the eigenmodes responsible for VSHF formation in the S3T system. (a) Growth rate as a function of the VSHF wavenumber $m_{U}$ for $\unicode[STIX]{x1D700}=0.25$ and two different excitation structures: $k_{e}/2\unicode[STIX]{x03C0}=6$ (dotted, $Fr=0.6$ ) and $k_{e}/2\unicode[STIX]{x03C0}=12$ (solid, $Fr=1.2$ ). (b) Growth rate as a function of $m_{U}$ and $\unicode[STIX]{x1D700}$ for $k_{e}/2\unicode[STIX]{x03C0}=6$ . Note the logarithmic $\unicode[STIX]{x1D700}$ axis. (c) Growth rate as a function of $m_{U}$ for $k_{e}/2\unicode[STIX]{x03C0}=12$ and four values of $N_{0}^{2}$ . (d) Growth rate of the fastest growing VSHF structure as a function of $N_{0}^{2}$ for $k_{e}/2\unicode[STIX]{x03C0}=6$ . This figure shows that the vertical wavenumber, $m_{U}$ , of the initially emergent VSHF is very sensitive to changes in the spectral structure of the excitation, and also that $m_{U}\rightarrow 0$ as $N_{0}^{2}\rightarrow \infty$ so that the initially emergent VSHF takes on the largest scale permitted by the domain if the stratification is sufficiently strong. Unless otherwise specified, parameters are as in figure 2.

Figure 14 summarizes how the dominant eigenvalue, $s$ , which is the VSHF growth rate, depends on $m_{U}$ and on the parameters $k_{e}$ , $N_{0}^{2}$ and $\unicode[STIX]{x1D700}$ . The dotted curve in panel (a) shows the VSHF growth rate as a function of $m_{U}$ for the standard parameter case with $\unicode[STIX]{x1D700}=0.25$ ( $Fr=0.6$ ). VSHFs with $1\leqslant m_{U}/2\unicode[STIX]{x03C0}\leqslant 10$ have positive growth rates, with the $m_{U}/2\unicode[STIX]{x03C0}=6$ structure having the fastest growth rate. This eigenvalue problem thus predicts that a VSHF with vertical wavenumber $m_{U}/2\unicode[STIX]{x03C0}=6$ will initially emerge from the turbulence, consistent with the structure of the VSHF discussed in § 2 and § 5. The solid curve in panel (a) shows the VSHF growth rates for the standard parameter case except with $k_{e}/2\unicode[STIX]{x03C0}=12$ ( $Fr=1.2$ ) so that the excitation is at a smaller scale. Increasing $k_{e}$ shifts the peak of the growth rate curve towards larger $m_{U}$ values, resulting in smaller scale VSHFs, and also enhances the VSHF growth rates.

Figure 14(b) shows the dependence of $s$ on $m_{U}$ and $\unicode[STIX]{x1D700}$ for the standard parameter case with $k_{e}/2\unicode[STIX]{x03C0}=6$ . For VSHFs with $1\leqslant m_{U}/2\unicode[STIX]{x03C0}\leqslant 11$ the growth rate becomes positive for sufficiently large $\unicode[STIX]{x1D700}$ . For VSHFs in this wavenumber band, the perturbation fluxes reinforce the infinitesimal VSHF and $s$ increases with increasing $\unicode[STIX]{x1D700}$ . For VSHFs outside this band, with $m_{U}/2\unicode[STIX]{x03C0}\geqslant 12$ , the perturbation fluxes oppose the VSHF so that the growth rate becomes increasingly negative as $\unicode[STIX]{x1D700}$ increases. The dashed line shows the stability boundary, $s=0$ . The homogeneous turbulence first becomes unstable near $\unicode[STIX]{x1D700}\approx 0.042$ ( $Fr=0.25$ ) to a VSHF with $m_{U}/2\unicode[STIX]{x03C0}=5$ . As $\unicode[STIX]{x1D700}$ increases, the growth rate of the $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF structure exceeds that of the $m_{U}/2\unicode[STIX]{x03C0}=5$ VSHF, so that for the standard parameter case with $\unicode[STIX]{x1D700}=0.25$ ( $Fr=0.6$ ) a VSHF with vertical wavenumber $m_{U}/2\unicode[STIX]{x03C0}=6$ initially emerges in the turbulence. We note that the emergent VSHF wavenumber varies between NL simulations depending on the particular realization of the stochastic excitation, with the $m_{U}/(2\unicode[STIX]{x03C0})=7$ structure occurring somewhat more frequently than $m_{U}/(2\unicode[STIX]{x03C0})=6$ . The existence of multiple turbulent equilibria in this system is predicted by S3T, as discussed in § 8. That the NL system often forms a VSHF with a slightly different scale than that predicted by linear stability analysis of the S3T system is likely due to the modification of the background spectrum of turbulence by perturbation–perturbation interactions. The influence of perturbation–perturbation interactions on the S3T stability of homogeneous turbulence has been analysed in detail by Constantinou et al. (Reference Constantinou, Farrell and Ioannou2014) in the context of barotropic beta-plane turbulence.

Figure 14(c) shows how $N_{0}^{2}$ influences the scale selection of the initially emergent VSHF. In panel (c) $s$ is shown as a function of $m_{U}$ for four $N_{0}^{2}$ values with $k_{e}/2\unicode[STIX]{x03C0}=12$ and $\unicode[STIX]{x1D700}=0.25$ . As $N_{0}^{2}$ increases, $s$ decreases and the peak (indicated by the vertical lines) shifts towards smaller $m_{U}$ . For very large $N_{0}^{2}$ the largest values of $s$ occur for VSHFs at the domain scale with $m_{U}/2\unicode[STIX]{x03C0}=1$ . However, unless $\unicode[STIX]{x1D700}$ is also very large the homogeneous turbulent state will remain stable and a domain-scale VSHF will not emerge, because $s$ decreases as $N_{0}^{2}$ becomes large. We note that the decrease of the VSHF wavenumber as $N_{0}^{2}$ increases demonstrates that $m_{U}$ is not directly related to either the Ozmidov wavenumber, $k_{O}=(N_{0}^{3}/\unicode[STIX]{x1D700})^{1/2}$ , or the buoyancy wavenumber, $k_{b}=N_{0}/\sqrt{\unicode[STIX]{x1D700}}$ , both of which increase as $N_{0}^{2}$ is increased, and also that the S3T prediction of $m_{U}$ depends on the parameter values and is not always equal to the excitation wavenumber, $k_{e}$ . As $N_{0}^{2}$ is decreased towards moderate and weak stratification (not shown), the wavenumber of the initially emergent VSHF remains near $k_{e}$ , consistent with the results of NL simulations shown in figure 6(c,d).

The dependence of $s$ on $N_{0}^{2}$ is shown directly in panel (d), which shows $\text{max}[s(m_{U})]$ , where the maximum is taken over $m_{U}$ , as a function of $N_{0}^{2}$ . For small $N_{0}^{2}$ , all modes have negative growth rates. This result provides an explanation for the frequent observation in numerical simulations that the VSHF ceases to emerge when the stratification becomes sufficiently weak. However, this result depends on the details of the stochastic excitation. In appendix A we describe a reduced model in which the excitation is anisotropic and which has the property that $s$ remains positive as $N_{0}^{2}\rightarrow 0$ , a result which was also obtained for similarly anisotropic excitation by Bakas & Ioannou (Reference Bakas and Ioannou2011). As $N_{0}^{2}$ increases from zero, $s$ increases to a maximum near $N_{0}^{2}\approx 10^{2}$ . This increase in growth rate is associated with the strengthening of the feedback between the VSHF and the turbulence described in § 3. The dependence of the S3T wave–mean flow feedback for harmonic mean structures on the parameter that sets the wave restoring force has been explained analytically in terms of wave dynamics by Bakas & Ioannou (Reference Bakas and Ioannou2013) for the case of barotropic beta-plane turbulence. For $N_{0}^{2}\gtrsim 10^{3}$ the growth rate falls off as ${\sim}1/N_{0}^{2}$ and approaches a constant asymptotic value as $N_{0}^{2}\rightarrow \infty$ that is set by the dissipation parameters.

Figure 15. Time evolution of the VSHF in two example simulations illustrating the correspondence between the behaviour of the NL system and the predictions of the linear stability analysis of the S3T system as the parameters are varied. Unless otherwise stated all parameters are as in figure 2. (a) An example with smaller scale excitation, $k_{e}/(2\unicode[STIX]{x03C0})=12$ ( $Re_{b}=10.4$ , $Fr=1.2$ ). The VSHF forms with $m_{U}/(2\unicode[STIX]{x03C0})\approx 12$ . (b) An example with smaller scale excitation, $k_{e}/(2\unicode[STIX]{x03C0})=12$ , and also stronger stratification, $N_{0}^{2}=5\times 10^{3}$ ( $Re_{b}=2.1$ , $Fr=0.53$ ). The VSHF forms with $m_{U}/(2\unicode[STIX]{x03C0})\approx 10$ . This figure demonstrates that, in the NL system, the VSHF forms at smaller scale when the turbulence is excited at smaller scale and that the VSHF forms at larger scale when the stratification is increased, consistent with the predictions of S3T.

Figure 15 shows the time evolution of the VSHF, $U$ , in two example simulations that illustrate the correspondence of the behaviour of the NL system with the predictions of linear stability analysis of the S3T system shown in figure 14. Panel (a) shows the VSHF evolution in an example in which the parameters are as in the standard case simulation but with the excitation scale modified to $k_{e}/(2\unicode[STIX]{x03C0})=12$ ( $Fr=1.2$ ). Consistent with the shift of peak of the VSHF growth rate curve to $m_{U}/(2\unicode[STIX]{x03C0})=12$ in figure 14(a), the emergent VSHF in the NL system has $m_{U}/(2\unicode[STIX]{x03C0})\approx 12$ . Panel (b) shows the VSHF evolution in an example with the same parameters as in panel (a), but with the stratification increased to $N_{0}^{2}=5\times 10^{3}$ ( $Fr=0.53$ ). Consistent with the shift of the peak of the VSHF growth rate curve towards lower values of $m_{U}$ under increased stratification in figure 14(c), the emergent VSHF in the NL system has $m_{U}/(2\unicode[STIX]{x03C0})\approx 10$ . Although the peak of $s(m_{U})$ in the S3T system occurs at the slightly larger scale $m_{U}/(2\unicode[STIX]{x03C0})=8$ for $N_{0}^{2}=5\times 10^{3}$ , the growth rates of the nearby VSHF wavenumbers are very similar to the peak value, as shown in figure 14(c), so that the VSHF wavenumber that is observed in a given realization of the stochastic system for these parameter values is likely to depend on the particular noise realization.

Results of this section demonstrate that the scale of the initially emergent VSHF, $m_{U}$ , is strongly influenced by the spectral structure of the perturbation field, which in our problem is set by $k_{e}$ . As the stratification becomes very strong, the VSHF scale is modified from the scale set by the excitation and tends towards the largest scale allowed by the domain. In realistic turbulence, the implication of this result is that we expect the spectral characteristics of the background turbulence to imprint strongly on the VSHF scale if the turbulence is sufficiently close to the stability boundary and the stratification is not too strong. We emphasize that the linear stability analysis conducted in this section provides a prediction of the scale of the initially emergent VSHF, rather than of the scale of the statistical equilibrium VSHF. For excitation strengths sufficiently near the stability boundary, the prediction based on linear stability analysis is expected to agree with the equilibrium VSHF structure. As the excitation strength is increased beyond the stability boundary, the structure of the VSHF may be modified from the initially emergent structure, as suggested by the NL simulation shown in figure 6(b). Previous studies of VSHF emergence have primarily been conducted using weak dissipation or strong excitation, so that the excitation strength lies well beyond the stability boundary, and these studies have also observed VSHFs with larger scale than that of the excitation. For example, Herring & Métais (Reference Herring and Métais1989) obtained a $m_{U}=6$ VSHF in 3D stratified turbulence maintained with $k_{e}\approx 11$ excitation, Smith (Reference Smith, Milewski, Smith, Waleffe and Tabak2001) obtained a VSHF with energy concentrated near $m_{U}\approx 10-15$ in the 2D system using $k_{e}\approx 96$ , and Smith & Waleffe (Reference Smith and Waleffe2002) obtained a VSHF with energy concentrated near $m_{U}\approx 9{-}11$ in the 3D system maintained with $k_{e}\approx 24$ . Although precise parameter correspondence between our study and these previous studies is difficult to establish due to differences in model formulation, these previous examples demonstrate that VSHFs with larger scale than that of the excitation are often observed when the system is strongly excited. This behaviour is expected based on analysis of the S3T system, which we discuss in § 7.

7 Equilibration of horizontal mean structure

In § 5 we showed that the S3T system initialized with a perturbative VSHF with $m_{U}/2\unicode[STIX]{x03C0}=6$ in the standard parameter case evolves into an equilibrium state with the same VSHF wavenumber (figures 9 and 10). We now analyse how the structure of this finite-amplitude equilibrium depends on the control parameters.

Figure 16(a, solid curve) shows the maximum value of the $m_{U}/2\unicode[STIX]{x03C0}=6$ equilibrium $U$ structure, maximized over $z$ , as a function of $\unicode[STIX]{x1D700}$ . The dotted curve shows an estimate of $U$ from a simple momentum balance model which we will explain later in this section. As suggested by the stability analysis in § 6, the $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF forms near $\unicode[STIX]{x1D700}\approx 0.044$ ( $Fr=0.25$ ) when the growth rate of the corresponding eigenmode crosses zero. The bifurcation is supercritical, with the VSHF equilibrating at weak amplitude just beyond the bifurcation point. Near the bifurcation point, $U$ increases rapidly with $\unicode[STIX]{x1D700}$ , with this rate of increase slowing as $\unicode[STIX]{x1D700}$ increases.

Figure 16. Equilibrium structure diagnostics for S3T and the simple momentum balance model in the case of the $m_{U}/2\unicode[STIX]{x03C0}=6$ horizontal mean state as a function of $\unicode[STIX]{x1D700}$ . (a) Maximum value of $U$ , maximized over $z$ , for the stable S3T fixed point with $m_{U}/2\unicode[STIX]{x03C0}=6$ (solid). This fixed point becomes secondarily unstable near $\unicode[STIX]{x1D700}=0.55$ ( $Fr=0.88$ ) and the dashed continuation shows the amplitude of the unstable solution. The dotted curve shows the estimate of the amplitude of $U$ from the simple momentum balance model (see text). (b–e) The vertical structure of the horizontal mean state as in figure 5 with dotted curves indicating the $\unicode[STIX]{x1D700}=0.08$ ( $Fr=0.34$ ) state and solid curves indicating the $\unicode[STIX]{x1D700}=0.54$ ( $Fr=0.88$ ) state. This figure shows that weak equilibration of the VSHF is captured by the simple momentum balance model and that the $U$ and $\overline{N^{2}}$ structures, and their phase relationship to one another, vary as $\unicode[STIX]{x1D700}$ is increased. Unless otherwise specified, parameters are as in figure 2.

The structure of the horizontal mean state depends on $\unicode[STIX]{x1D700}$ . Figure 16(b–e, dotted curves) shows the horizontal mean structure of the marginally supercritical equilibrium at $\unicode[STIX]{x1D700}=0.08$ ( $Fr=0.34$ ). The VSHF structure is similar to that of the unstable $m_{U}/2\unicode[STIX]{x03C0}=6$ harmonic eigenmode. The phase relationship between $U$ and $\overline{N^{2}}$ differs from that found in the more strongly excited $\unicode[STIX]{x1D700}=0.25$ ( $Fr=0.6$ ) case discussed in §§ 2 and 5. In the $\unicode[STIX]{x1D700}=0.08$ ( $Fr=0.34$ ) case of weak equilibration the stratification is enhanced in the shear regions, rather than weakened, and $\overline{Ri}$ is large for all $z$ due to the weak shear. The solid curves in figure 16(b–e) show the horizontal mean structure of the $\unicode[STIX]{x1D700}=0.54$ ( $Fr=0.88$ ) equilibrium. For this more strongly supercritical equilibrium, the shear regions are characterized by weakened stratification and $\overline{Ri}<1/4$ . This structure is similar to that shown in figure 11 for the $\unicode[STIX]{x1D700}=0.25$ case, but with stronger shear and smaller $\overline{Ri}$ values. The VSHF remains hydrodynamically stable (i.e. all eigenvalues of $\unicode[STIX]{x1D63C}_{n}$ have negative real parts) despite having $\overline{Ri}<1/4$ due to the dissipation acting on the perturbation fields. When $\unicode[STIX]{x1D700}$ is further increased the $m_{U}/2\unicode[STIX]{x03C0}=6$ fixed point becomes secondarily unstable, indicated by the dashed continuation of the solid curve in figure 16(a).

The changing phase relationship between $U$ and $\overline{N^{2}}$ shown in figure 16 that occurs as a function of $\unicode[STIX]{x1D700}$ mirrors the change in this relationship shown in figure 9 that occurs as a function of time. Comparison of figure 9(c,d) shows that, when the developing VSHF is weak, $\overline{N^{2}}$ is enhanced where the shear is strongest. When the VSHF becomes strong, the stratification is reorganized by the turbulent fluxes so that $\overline{N^{2}}$ is weakest where the shear is strongest.

The mechanism of VSHF equilibration at weak amplitude can be understood using a simple momentum balance model based on the test function analysis of § 3. To construct the simple model we first approximate the horizontal mean state as $U=U_{0}\sin (m_{U}z)$ and $B=0$ , where $U_{0}$ is the equilibrium VSHF amplitude that we will estimate. We then estimate the acceleration of the VSHF produced by the induced perturbation momentum fluxes as a function of $U_{0}$ using (3.1)–(3.2). Our estimate of the equilibrium VSHF amplitude is the value of $U_{0}$ for which this acceleration is balanced by dissipation. As $\unicode[STIX]{x1D700}\rightarrow \unicode[STIX]{x1D700}_{c}$ this simple model becomes exact because both $U$ and the perturbation flux divergence become exactly harmonic and $B\rightarrow 0$ . For $\unicode[STIX]{x1D700}>\unicode[STIX]{x1D700}_{c}$ , the structure of $-\unicode[STIX]{x2202}_{z}\langle \overline{u^{\prime }w^{\prime }}\rangle$ deviates from harmonic and we estimate the equilibrium VSHF amplitude by projecting the acceleration onto the assumed harmonic VSHF structure.

Figure 17. Illustration of the simple momentum balance model for weakly supercritical VSHF equilibration for $\unicode[STIX]{x1D700}=0.08$ and $m_{U}/2\unicode[STIX]{x03C0}=6$ , with other parameters as in figure 2. The solid curve shows the projection of the perturbation momentum flux divergence, calculated using the test function analysis of § 3, onto the assumed harmonic VSHF structure. The dashed line shows the dissipation acting on the VSHF, given by $(r_{m}+\unicode[STIX]{x1D708}m_{U}^{2})U_{0}$ . The simple model estimate of the equilibrium VSHF amplitude is the value of $U_{0}$ at which these terms balance one another. The vertical dotted line indicates the equilibrium VSHF amplitude obtained from the full S3T system. This figure demonstrates that the dynamics of weakly supercritical VSHF equilibration is captured by the simple balance model.

We illustrate the simple momentum balance model for $\unicode[STIX]{x1D700}=0.08$ ( $Fr=0.34$ ) and $m_{U}/2\unicode[STIX]{x03C0}=6$ in figure 17, which shows the estimated acceleration (solid) and dissipation (dashed) of the VSHF as functions of the VSHF amplitude, $U_{0}$ . The dissipation, $(r_{m}+\unicode[STIX]{x1D708}m_{U}^{2})U_{0}$ , increases linearly with $U_{0}$ . For small $U_{0}$ the acceleration is stronger than the dissipation, consistent with spontaneous VSHF formation as a linear instability for these parameters. Due to the negative curvature of the acceleration as a function of $U_{0}$ the two terms balance near $U_{0}\approx 0.65$ , which gives the simple model estimate of the equilibrium VSHF amplitude. The vertical dotted line indicates the equilibrium VSHF strength in the full S3T system. For this value of $\unicode[STIX]{x1D700}$ the simple model captures the equilibration dynamics, implying that modification of $\overline{N^{2}}$ and changes in $U$ structure do not play important roles in the weak equilibration process. The simple model estimate of $\text{max}(U)$ as a function of $\unicode[STIX]{x1D700}$ is shown in figure 16(a) as the dotted curve. The model estimate matches the results of the full calculation as $\unicode[STIX]{x1D700}\rightarrow \unicode[STIX]{x1D700}_{c}$ and diverges from the full solution as $\unicode[STIX]{x1D700}$ increases.

Figure 18. Secondary instability of the S3T fixed point corresponding to the $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF for $\unicode[STIX]{x1D700}=1$ ( $Fr=1.2$ ), with other parameters as in figure 2. (a,b) The time evolution of (a) $U$ and (b) $\overline{N^{2}}$ . (c) The kinetic energy evolution. This figure shows that for strong excitation the $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF state is unstable to the development of a global vertical wavenumber 2 pattern in $U$ that is superimposed on the initial wavenumber 6 pattern, strengthening the VSHF and modifying its structure to produce wider shear regions.

As $\unicode[STIX]{x1D700}$ is increased, the global minimum of $\overline{Ri}$ falls further below $1/4$ and the $m_{U}/2\unicode[STIX]{x03C0}=6$ state becomes secondarily unstable just beyond $\unicode[STIX]{x1D700}=0.54$ ( $Fr=0.88$ ). Although this instability occurs when $U$ is near the laminar stability boundary, which is modified from $\overline{Ri}=1/4$ by the presence of dissipation and by our choice of a finite periodic domain which quantizes the permitted perturbation horizontal wavenumbers, we emphasize that this secondary instability is a property of the S3T dynamics, rather than of the perturbation dynamics determined by the operator $\unicode[STIX]{x1D63C}_{n}$ . In particular, the $m_{U}/2\unicode[STIX]{x03C0}=6$ state remains hydrodynamically stable at all times during the instability development. Figure 18(a,b) shows the time evolution of $U$ and $\overline{N^{2}}$ during the development of the secondary instability for $\unicode[STIX]{x1D700}=1$ ( $Fr=1.2$ ). The VSHF structure near $t=30$ reveals that the sixfold symmetry of the $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF is spontaneously broken by the instability. As the instability develops, the positive VSHF peaks near $z=0.2$ and $z=0.7$ contract and weaken while their neighbouring negative peaks strengthen and expand. Similarly, the negative VSHF peaks near $z=0.5$ and $z=0.9$ contract and weaken while their neighbouring positive VSHF peaks strengthen and expand. The particular locations of the strengthening and weakening features are the result of the symmetry breaking and so depend on the small perturbations included in the initialization. Figure 18(c) shows the evolution of kinetic energy during the instability. The changes in the VSHF structure are associated with an increase in the mean kinetic energy consistent with the broadening of the VSHF pattern allowing $U$ to strengthen while maintaining a hydrodynamically stable shear. Similar behaviour is shown to occur in the NL system in figure 6(b), which shows an example in which the effective excitation strength has been increased relative to the standard case integration by removing the Rayleigh drag terms from the dynamics. In this example, a VSHF with $m_{U}/(2\unicode[STIX]{x03C0})\approx 6$ initially emerges from the turbulence and this VSHF transitions to lower wavenumber as the integration is continued. Secondary instabilities of finite-amplitude mean shear flows that result in broader shear patterns also occur in the barotropic beta-plane system and have been analysed using S3T by Constantinou et al. (Reference Constantinou, Farrell and Ioannou2014).

The structure of the horizontal mean state before and after the development of the secondary instability for $\unicode[STIX]{x1D700}=1$ ( $Fr=1.2$ ) is shown in figure 19. Prior to the instability development ( $t=10$ , dotted) the structure is similar to that shown for $\unicode[STIX]{x1D700}=0.54$ ( $Fr=0.88$ ) in figure 16(b–e) and is characterized by a VSHF with $m_{U}/2\unicode[STIX]{x03C0}=6$ and weakened stratification in the shear extrema. The $U$ profile of the final equilibrium structure ( $t=50$ , solid) contains shear regions with two distinct widths which are associated with distinct phase relationships between $U$ and $\overline{N^{2}}$ . For the wider shear regions, the $U$ profile inflects in the centre of the shear region and $\overline{N^{2}}$ is locally maximized there, resulting in $\overline{Ri}\gtrsim 1/4$ . The narrower shear regions are similar to those that precede the secondary instability development and have $\overline{Ri}<1/4$ .

Figure 19. Vertical structure of the horizontal mean state in the S3T system before and after the development of the secondary instability for $\unicode[STIX]{x1D700}=1$ ( $Fr=1.2$ ), with other parameters as in figure 2. Panels are as in figure 5 with dotted curves showing the structure for $t=10$ and solid curves showing the structure for $t=50$ . This figure shows how the structure of the horizontal mean state is reorganized by the secondary instability. The unstable equilibrium state at $t=10$ has $\overline{Ri}<1/4$ in regions of strongest shear and weakest stratification. The final equilibrium state has shear regions of two different widths in which the broader shear regions have $\overline{Ri}>1/4$ due to enhanced stratification and weakened shear in the cores of the shear regions.

8 Multiple turbulent equilibria in stratified turbulence

In § 6 we showed an NL simulation in which a $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF emerges, corresponding to the eigenmode of the linearized S3T system that has the fastest growth rate in the standard parameter case. However, figure 14(a) shows that, for the standard parameter case, all VSHF structures in the wavenumber band $1\leqslant m_{U}/2\unicode[STIX]{x03C0}\leqslant 10$ have positive growth rates. The subdominant eigenmodes (i.e. those with $m_{U}/2\unicode[STIX]{x03C0}\neq 6$ ) continue to finite-amplitude VSHF equilibria at the corresponding wavenumbers. These equilibria may or may not be stable. In this section we demonstrate that multiple turbulent equilibrium states are possible in 2D Boussinesq turbulence by providing an example of such an alternate stable equilibrium in the S3T and NL systems.

Figure 20. Time evolution of the horizontal mean structure of the $m_{U}/2\unicode[STIX]{x03C0}=4$ equilibrium state in the NL and S3T systems. Panels show the time evolution of (a) $U$ and (b) $\overline{N^{2}}$ in the NL system and (c) $U$ and (d) $\overline{N^{2}}$ in the S3T system. This figure shows that when initialized with a finite-amplitude VSHF with wavenumber $m_{U}/2\unicode[STIX]{x03C0}=4$ the NL system maintains this structure, resulting in a turbulent equilibrium state different from that discussed in § 2 for the same parameter values, and that this alternate equilibrium state is also a fixed point of the S3T system. Parameters are as in figure 2.

In figure 20(a,c) we show the development of $U$ in the NL and S3T systems in an example in which the parameters are set to the standard values ( $Fr=0.6$ as in figure 2) but the initial conditions are chosen to initiate a VSHF with wavenumber $m_{U}/2\unicode[STIX]{x03C0}=4$ . The NL system is initialized with a mean flow $U\propto \sin (m_{U}z)$ for $m_{U}/2\unicode[STIX]{x03C0}=4$ and the S3T system is initialized with the same $U$ profile and $\unicode[STIX]{x1D63E}_{n}=\unicode[STIX]{x1D63E}_{n}^{\star }$ . In the S3T system this initial condition evolves into a stable $m_{U}/2\unicode[STIX]{x03C0}=4$ fixed point. We note that, as shown in § 6, the S3T system will evolve, in the standard parameter case, towards the $m_{U}/2\unicode[STIX]{x03C0}=6$ fixed point for any sufficiently small initial perturbation, but that in this example the system evolves towards the $m_{U}/2\unicode[STIX]{x03C0}=4$ fixed point as a result of the finite initial perturbation. In the NL system the $m_{U}/2\unicode[STIX]{x03C0}=4$ turbulent equilibrium is maintained for the length of the integration. Due to noise in the NL system, the turbulence may eventually transition to another equilibrium state, such as the $m_{U}/2\unicode[STIX]{x03C0}=6$ state discussed in § 2. The development of $\overline{N^{2}}$ for this example is shown in figure 20(b,d). As in the previous examples of equilibria, $\overline{N^{2}}$ has a doubled vertical wavenumber relative to $U$ and is more variable than $U$ in the NL system. The vertical structure of the horizontal mean state is shown in figure 21. The VSHF structure (panels a,b) resembles a sawtooth in both the NL (dotted) and S3T (solid) systems. The phase relationship between $U$ and $\overline{N^{2}}$ (panel c) shares some features with that shown in figure 16 for the $m_{U}/2\unicode[STIX]{x03C0}=6$ equilibrium with $\unicode[STIX]{x1D700}=0.54$ ( $Fr=0.88$ ). In particular, the weakest values of $\overline{N^{2}}$ occur in the centres of the shear regions. The excitation strength $\unicode[STIX]{x1D700}=0.25$ falls in a transitional range for the $m_{U}/2\unicode[STIX]{x03C0}=4$ equilibrium in which $\overline{N^{2}}$ has an approximately $m_{B}/2\unicode[STIX]{x03C0}=16$ structure. As $\unicode[STIX]{x1D700}$ is increased (not shown), the stratification in the shear centres continues to weaken, producing density layers at these locations, and the stratification near the VSHF peaks is enhanced.

Figure 21. Vertical structure of the horizontal mean state of the $m_{U}/2\unicode[STIX]{x03C0}=4$ equilibrium in the NL and S3T systems. Panels are as in figure 5 with dotted curves showing the time-averaged structure over $t\in [22,45]$ for the NL system and solid curves showing the final fixed point structure for the S3T system. Parameters are as in figure 2.

Figure 22 shows the evolution of kinetic energy in the NL and S3T systems. Consistent with the results for the $m_{U}/2\unicode[STIX]{x03C0}=6$ equilibrium in § 5, the equilibrium value of $\overline{K}$ in the S3T system exceeds that of NL system. In both systems, the broader VSHF in the $m_{U}/2\unicode[STIX]{x03C0}=4$ equilibrium is more energetic than the VSHF in the $m_{U}/2\unicode[STIX]{x03C0}=6$ equilibrium. This is consistent with the behaviour shown in figure 18(c), in which the broadened VSHF resulting from the secondary instability is more energetic than the $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF that precedes the instability.

Figure 22. Kinetic energy evolution in the NL and S3T systems initialized with a $m_{U}/2\unicode[STIX]{x03C0}=4$ VSHF. This figure shows that, as for the $m_{U}/2\unicode[STIX]{x03C0}=6$ equilibrium, the VSHF in the S3T system is more energetic than the VSHF in the NL system, and comparison with figure 10 shows that in both the NL and S3T systems the $m_{U}/2\unicode[STIX]{x03C0}=4$ VSHF is more energetic than the $m_{U}/2\unicode[STIX]{x03C0}=6$ VSHF. Parameters are as in figure 2.

9 Reflection of the S3T bifurcation in the NL and QL systems

In § 5 we compared the behaviour of the NL, QL and S3T systems with all parameter values fixed. Comparing the three systems in this way allows for a detailed comparison of the structures of the mean state and of the turbulent spectra to be made. However, our analysis of the S3T system has revealed phenomena, including the bifurcation associated with the initial formation of the VSHF, that can be analysed only by allowing variation of the control parameters. We now compare the behaviour of the three systems as a function of the excitation strength, $\unicode[STIX]{x1D700}$ , in terms of the fraction of the total kinetic energy of the flow that is associated with the VSHF. We define this fraction as $\text{zmf}=\overline{K}/(\overline{K}+K^{\prime })$ , for zonal mean flow (zmf) index, borrowing this definition from studies of barotropic beta-plane turbulence (Srinivasan & Young Reference Srinivasan and Young2012; Constantinou et al. Reference Constantinou, Farrell and Ioannou2014). We note that in the context of barotropic turbulence an alternative approach based on regime diagrams has also been used to characterize the transition of turbulence to states dominated by zonal jets (Galperin et al. Reference Galperin, Sukoriansky and Dikovskaya2010).

Figure 23. Equilibrium zmf indices in the NL, QL and S3T systems as functions of $\unicode[STIX]{x1D700}$ , with other parameters as in figure 2. The zmf index measures the fraction of the total kinetic energy of the flow that is associated with the VSHF. This figure shows that the bifurcation through which the VSHF forms in the deterministic S3T system is reflected in the behaviour of the NL and QL systems, which show an abrupt increase in the fraction of the total kinetic energy contained in the VSHF near the S3T bifurcation point.

Figure 23 shows the equilibrium zmf value as a function of $\unicode[STIX]{x1D700}$ for the NL, QL and S3T systems. As all three systems possess multiple equilibria, there is some ambiguity as to the meaning of the equilibrium energies. For the S3T system we show the maximum zmf obtained when the system is initialized with a perturbative VSHF at each unstable wavenumber $m_{U}$ . For the NL and QL systems we initialize from rest and show the time average of zmf over the final 10 time units of a $t\in [0,450]$ integration. As many long integrations are required for this comparison, the simulations shown in this section are spun up at low resolution and the resulting turbulence is interpolated to the standard resolution of $512^{2}$ grid points to initialize a simulation of the final 10 time units.

As discussed in § 6, the S3T system passes through a bifurcation near $\unicode[STIX]{x1D700}\approx 0.04$ ( $Fr=0.24$ ). This bifurcation is reflected in the zmf indices of the QL and NL systems. For $\unicode[STIX]{x1D700}\lesssim 0.04$ the VSHF accounts for only a few percent of the total kinetic energy of the flow. As $\unicode[STIX]{x1D700}$ increases beyond the S3T bifurcation point, the zmf index increases rapidly and the VSHF becomes energetically dominant. As was found in § 5, the S3T VSHF is the most energetic and the QL VSHF tends to be more energetic than the NL VSHF. The eventual decrease of the zmf indices in the QL and NL systems as $\unicode[STIX]{x1D700}$ is increased may be due to the tendency of those systems to maintain VSHF structures with $m_{U}/2\unicode[STIX]{x03C0}=6$ , even when this is not the most energetic VSHF structure. The S3T curve does not show this decrease as we choose the most energetic VSHF equilibrium to define the S3T equilibrium zmf. This maximally energetic VSHF equilibrium often has a lower vertical wavenumber than that of the fastest growing eigenmode, as discussed in § 8. Developing a complete understanding of the behaviour of, and correspondence between, the QL, NL and S3T systems in the limit of strong excitation is beyond the scope of the present study but is an important avenue for future investigation. Note also in figure 23 the characteristic increase in fluctuating VSHF amplitude in the NL and QL cases as the bifurcation point is approached. This results from excitation of the reflection in QL and NL of the stable modes of the S3T system, which have no analytical expression in the QL and NL systems themselves. These modes are excited by the noise inherent in the QL and NL systems, while no such excitation is seen in the noise-free S3T system (Constantinou et al. Reference Constantinou, Farrell and Ioannou2014).