1 Introduction
Horizontal layers of density or buoyancy are commonly identified in observations of stratified turbulence in the ocean (Gregg Reference Gregg1980; Duda & Rehmann Reference Duda and Rehmann2002), the atmosphere (Cho et al.
Reference Cho, Newell, Anderson, Barrick and Lee Thornhill2003), laboratory studies (Manucharyan & Caulfield Reference Manucharyan and Caulfield2015) and simulations (Lucas, Caulfield & Kerswell Reference Lucas, Caulfield and Kerswell2017; Taylor & Zhou Reference Taylor and Zhou2017; Fitzgerald & Farrell Reference Fitzgerald and Farrell2018a
). Several physical mechanisms have been advanced to explain layering phenomena, including the inertial instability (Ménesguen et al.
Reference Ménesguen, Hua, Fruman and Schopp2009), double-diffusive instability (Radko Reference Radko2003) and the Phillips instability (PI) (Phillips Reference Phillips1972; Posmentier Reference Posmentier1977). Of these, PI presents perhaps the greatest challenge to theoretical understanding because the mechanism, in turbulence, is inherently statistical. (In the sense of arising from the response of the mean correlation between velocity and buoyancy perturbations to changes in stratification. We note, however, that the dynamics of buoyancy layering have also been investigated in the context of analytically tractable laminar stratified flows (e.g. Balmforth & Young Reference Balmforth and Young2005).) The basis of PI is an assumption about the flux-gradient relationship between the turbulent buoyancy flux,
$\overline{w^{\prime }b^{\prime }}$
, and the background stratification or mean buoyancy gradient,
$N^{2}\equiv \unicode[STIX]{x2202}\bar{b}/\unicode[STIX]{x2202}z$
. If the flux weakens as stratification increases, small-scale stratification structure is predicted to emerge spontaneously. To see this, suppose that the flux-gradient relation is given by
$$\begin{eqnarray}\overline{w^{\prime }b^{\prime }}={\mathcal{F}}(N^{2}).\end{eqnarray}$$
Absent advection and forcing/dissipation the mean buoyancy and stratification evolve as
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\bar{b}=-\unicode[STIX]{x2202}_{z}\overline{w^{\prime }b^{\prime }},\quad \unicode[STIX]{x2202}_{t}N^{2}=-\unicode[STIX]{x2202}_{zz}^{2}{\mathcal{F}}(N^{2}),\end{eqnarray}$$
where overbars denote the horizontal mean and primes denote perturbations. From (1.1), perturbations,
$\unicode[STIX]{x1D6FF}N^{2}$
, to an equilibrium uniform stratification,
$N_{\star }^{2}$
, modify the buoyancy flux as
$$\begin{eqnarray}{\mathcal{F}}(N_{\star }^{2}+\unicode[STIX]{x1D6FF}N^{2})\approx {\mathcal{F}}(N_{\star }^{2})+(\unicode[STIX]{x2202}_{N^{2}}{\mathcal{F}})|_{N_{\star }^{2}}\unicode[STIX]{x1D6FF}N^{2},\end{eqnarray}$$
to first order in
$\unicode[STIX]{x1D6FF}N^{2}$
. Linearizing (1.2b
) reveals that
$\unicode[STIX]{x1D6FF}N^{2}$
evolves as
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D6FF}N^{2}\approx -(\unicode[STIX]{x2202}_{N^{2}}{\mathcal{F}})|_{N_{\star }^{2}}\unicode[STIX]{x2202}_{zz}^{2}\unicode[STIX]{x1D6FF}N^{2}.\end{eqnarray}$$
If the flux-gradient relation has positive slope, as in the case of downgradient buoyancy fluxes which weaken as stratification increases, then the diffusion coefficient in (1.4) is negative and small-scale structure is predicted to develop.
Identification and prediction of PI has proven to be challenging in practice. Density layer formation has been observed in several laboratory studies in which stratified tanks were stirred by vertical rods (Ruddick, McDougall & Turner Reference Ruddick, McDougall and Turner1989; Park, Whitehead & Gnandadesikan Reference Park, Whitehead and Gnandadesikan1994; Holford & Linden Reference Holford and Linden1999). Although Ruddick et al. (Reference Ruddick, McDougall and Turner1989) and Park et al. (Reference Park, Whitehead and Gnandadesikan1994) attributed layer formation to PI, Thorpe (Reference Thorpe2016) has recently challenged this explanation, arguing that layering in these studies, and in the study of Holford & Linden (Reference Holford and Linden1999), instead results from the interaction of internal shear waves with decaying wake vortices left by the moving rods. Accurate attribution of layering to PI requires the difficult diagnostic task of evaluating the flux-gradient relation in evolving turbulence. Recently, progress has been made by combining idealized simulations with a parameterization of mixing efficiency (Salehipour et al. Reference Salehipour, Peltier, Whalen and MacKinnon2016; Taylor & Zhou Reference Taylor and Zhou2017). An alternate approach based on an implementation of statistical state dynamics (SSD) referred to as S3T (short for stochastic structural stability theory, Farrell & Ioannou (Reference Farrell and Ioannou2003, Reference Farrell and Ioannou2007)) has previously been used to provide a theory for the flux-gradient relation for momentum fluxes in barotropic turbulence (Srinivasan & Young Reference Srinivasan and Young2014). Here we apply the methods of S3T analysis previously developed by Srinivasan & Young (Reference Srinivasan and Young2012, Reference Srinivasan and Young2014) to analyse the novel problem of buoyancy fluxes in stratified turbulence and PI.
In SSD, equations of motion are written directly for statistical variables of the system, such as the cumulants of the state variables. This approach must reckon with the closure problem: the dynamics of cumulants at a given order involves higher-order cumulants. In many turbulent systems of interest, however, the mechanism determining the equilibrium state primarily involves the interaction of the quadratic fluxes arising from the incoherent component of the turbulence with the coherent large-scale structure (Farrell & Ioannou Reference Farrell, Ioannou, Galperin and Read2017). In such cases it turns out that the essential dynamics is contained within the coupled dynamics of the mean and the second cumulant only, with higher cumulants being inessential. This dynamical insight motivates closure of SSD at second order by stochastically parameterizing third-order cumulants. This is the basis of the S3T approach to SSD. S3T has been successful in capturing and providing mechanistic explanation for the dynamics of a variety of turbulent systems exhibiting large-scale organization, including wall-bounded shear flows (e.g. Farrell & Ioannou Reference Farrell and Ioannou2012), rotating magnetohydrodynamics (e.g. Constantinou & Parker Reference Constantinou and Parker2018), and planetary turbulence (e.g. Constantinou, Farrell & Ioannou Reference Constantinou, Farrell and Ioannou2014). In previous work, we applied S3T to analyse the spontaneous emergence of turbulent jets known as vertically sheared horizontal flows in stochastically maintained homogeneous stratified turbulence (Fitzgerald & Farrell Reference Fitzgerald and Farrell2018a ,Reference Fitzgerald and Farrell b ). Jets were found to develop in turbulence without requiring the coherent large scales of the flow to be prescribed or directly forced. In Fitzgerald & Farrell (Reference Fitzgerald and Farrell2018a ), we focused on the finite-amplitude turbulent equilibria, while in Fitzgerald & Farrell (Reference Fitzgerald and Farrell2018b ) we analysed in detail, at the perturbative level, the wave mean flow feedback mechanism responsible for jet emergence. We note that although buoyancy layers emerge in the simulations of Fitzgerald & Farrell (Reference Fitzgerald and Farrell2018a ,Reference Fitzgerald and Farrell b ), these layers form via a mechanism different from PI; namely, organization of buoyancy fluxes by finite-amplitude layered jets producing a buoyancy layer pattern that is phase-synchronized with the jet pattern.
In this work we apply S3T to the two-dimensional (2D) Boussinesq model of stratified turbulence maintained by stochastic excitation in the presence of an imposed background shear flow. The 2D equations are used to reduce the computational burden associated with the application of S3T, but S3T can also be applied to the corresponding 3D system. We show that the flux-gradient relation of 2D Boussinesq sheared stratified turbulence has a region of positive slope consistent with the requirements of PI. We further show that the uniformly stratified state is unstable to layer formation when the stratification is in the positively sloping region. Although this instability has analytical expression only in S3T and not in the dynamics of realizations of the nonlinear Boussinesq system, simulated realizations of the nonlinear Boussinesq system confirm that buoyancy layers also emerge spontaneously in that system corresponding to the predictions made by S3T.
2 Formulation
2.1 Model configuration
We investigate layer formation in 2D Boussinesq sheared stratified turbulence maintained by statistically homogeneous stochastic excitation. The background stratification is uniform and a uniform shear is imposed by a fixed horizontal mean flow,
$U(z)=\unicode[STIX]{x1D6EC}z$
. In a Reynolds decomposition based on the horizontal mean the equations of motion are
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}B=-\unicode[STIX]{x2202}_{z}\overline{w^{\prime }b^{\prime }}-r_{m}B+\unicode[STIX]{x1D708}\unicode[STIX]{x2202}_{zz}^{2}B, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D701}^{\prime }=-(\unicode[STIX]{x1D6EC}z)\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D701}^{\prime }+\unicode[STIX]{x2202}_{x}b^{\prime }-[J(\unicode[STIX]{x1D713}^{\prime },\unicode[STIX]{x1D701}^{\prime })-\overline{J(\unicode[STIX]{x1D713}^{\prime },\unicode[STIX]{x1D701}^{\prime })}]-r\unicode[STIX]{x1D701}^{\prime }+\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\unicode[STIX]{x1D701}^{\prime }+\sqrt{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D709}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}b^{\prime }=-(\unicode[STIX]{x1D6EC}z)\unicode[STIX]{x2202}_{x}b^{\prime }-w^{\prime }(N_{0}^{2}+\unicode[STIX]{x2202}_{z}B)-[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 }. & \displaystyle\end{eqnarray}$$
Here
$J(f,g)\equiv (\unicode[STIX]{x2202}_{x}f)(\unicode[STIX]{x2202}_{z}g)-(\unicode[STIX]{x2202}_{z}f)(\unicode[STIX]{x2202}_{x}g)$
is the Jacobian,
$\unicode[STIX]{x0394}=\unicode[STIX]{x2202}_{xx}^{2}+\unicode[STIX]{x2202}_{zz}^{2}$
is the Laplacian,
$\unicode[STIX]{x1D713}$
is the streamfunction satisfying
$(u,w)=(-\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D713},\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713})$
,
$\unicode[STIX]{x1D701}=\unicode[STIX]{x2202}_{x}w-\unicode[STIX]{x2202}_{z}u=\unicode[STIX]{x0394}\unicode[STIX]{x1D713}$
is the vorticity, overbars denote horizontal averages and primes denote perturbations relative to the horizontal average,
$b$
is the buoyancy relative to the background stratification,
$B\equiv \overline{b}$
denotes the horizontal mean buoyancy relative to the background stratification,
$N_{0}$
is the uniform background buoyancy frequency, and
$\sqrt{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D709}$
is the stochastic excitation applied to the perturbation vorticity field, where
$\unicode[STIX]{x1D700}$
sets the strength of the excitation. Linear damping relaxes the mean stratification towards the background value,
$N_{0}$
. Relaxation provides a simple parameterization of the physical processes that maintain the large-scale stratification in particular examples of stratified turbulence, such as radiative–convective equilibration in the atmosphere or advection–diffusion in the ocean (Munk Reference Munk1966; Manabe & Wetherald Reference Manabe and Wetherald1967). Linear damping is also used to parameterize dissipation of the perturbation fields and serves as a simple model of processes such as turbulent dissipation that are not explicitly contained in our model. The mean state relaxation rate,
$r_{m}$
, is chosen to be weaker than the perturbation damping rate,
$r$
. Weak diffusion, with coefficient
$\unicode[STIX]{x1D708}$
, is included to ensure numerical convergence and is applied equally to all fields. Maintenance of stratified turbulence using stochastic excitation and linear damping allows the turbulent equilibrium state to be established quickly and facilitates investigation of PI. In previous work the sensitivity of large-scale structure emergence in stratified turbulence to the choice of linear or diffusive dissipation was directly examined and it was found that similar structures emerge in both cases (Fitzgerald & Farrell Reference Fitzgerald and Farrell2018a
). We refer to the nonlinear Boussinesq equations (2.1)–(2.3) as NL.
The model domain and excitation structure are shown in figure 1. Our numerical simulations use a doubly periodic rectangular domain
$(x,z)\in [-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0}]\times [-2\unicode[STIX]{x03C0},2\unicode[STIX]{x03C0}]$
, with the additional length in
$z$
allowing for a shear reversal to satisfy periodic boundary conditions. We analyse the interior region
$(x,z)\in [-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0}]\times [-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0}]$
in which the shear is constant (panels
$c$
,
$d$
, shaded). Additional strong mean state relaxation is used in the regions of shear reversal, augmenting the weak mean state relaxation used throughout the domain, to reduce the influence of the boundary conditions and facilitate analysis of buoyancy layering in a constant background shear. The stochastic excitation is white in time and excites a subset of the horizontal wavenumber components of the flow, each with equal energy (panel
$b$
). Details and parameter values are provided in § 3.
Figure 1. (a) A realization of the stochastic vorticity excitation, normalized as
$\unicode[STIX]{x1D709}/\max [\unicode[STIX]{x1D709}]$
. (b) Energy injection spectrum of the excitation, normalized as
$\tilde{\unicode[STIX]{x1D700}}(p,q)/\text{max}[\tilde{\unicode[STIX]{x1D700}}(p,q)]$
.
$\tilde{\unicode[STIX]{x1D700}}(p,q)$
is doubly mirror symmetric and the first quadrant is shown. (c,d) Imposed profiles of horizontal mean velocity,
$U$
, and shear,
$\unicode[STIX]{x2202}U/\unicode[STIX]{x2202}z$
. The central region (
$z\in [-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0}]$
, shaded) is analysed. (e) Profile of the mean state relaxation coefficient,
$r_{m}$
, which is applied to the horizontal mean buoyancy anomaly,
$B$
. Strong relaxation is used at the edges of the analysis region and a small background value
$r_{m}=5\times 10^{-2}$
is used elsewhere. See § 3 for non-dimensionalization details.
2.2 Statistical state dynamics
We analyse layer formation using S3T (Farrell & Ioannou Reference Farrell and Ioannou2003) in the differential form developed by Srinivasan & Young (Reference Srinivasan and Young2012). S3T is a canonical second-order closure incorporating two assumptions. The first is the ergodic assumption that horizontal and ensemble averages can be equated so that
$\overline{f}=\langle f\rangle$
, where angle brackets indicate the ensemble mean. The second assumption is that third cumulants can be parameterized stochastically. The simplest parameterization sets third cumulants to zero, which suffices for our study of PI. To formulate S3T we apply these assumptions to (2.1)–(2.3) to derive equations of motion for the two-point, equal-time covariance functions
$$\begin{eqnarray}Z(\boldsymbol{x}_{1},\boldsymbol{x}_{2},t)\equiv \langle \unicode[STIX]{x1D701}_{1}^{\prime }(t)\unicode[STIX]{x1D701}_{2}^{\prime }(t)\rangle ,\quad T(\boldsymbol{x}_{1},\boldsymbol{x}_{2},t)\equiv \langle b_{1}^{\prime }(t)b_{2}^{\prime }(t)\rangle ,\end{eqnarray}$$
$$\begin{eqnarray}\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}(\boldsymbol{x}_{1},\boldsymbol{x}_{2},t)\equiv \langle \unicode[STIX]{x1D701}_{1}^{\prime }(t)b_{2}^{\prime }(t)\rangle ,\quad \unicode[STIX]{x1D6E4}^{b}(\boldsymbol{x}_{1},\boldsymbol{x}_{2},t)\equiv \langle b_{1}^{\prime }(t)\unicode[STIX]{x1D701}_{2}^{\prime }(t)\rangle ,\end{eqnarray}$$
where
$f_{1}(t)\equiv f(\boldsymbol{x}_{1},t)$
. We write the S3T dynamics using the sum and difference coordinates
$\overline{\boldsymbol{x}}=(\bar{x},\bar{z})\equiv (\boldsymbol{x}_{1}+\boldsymbol{x}_{2})/2$
and
$\tilde{\boldsymbol{x}}=(\tilde{x},\tilde{z})\equiv \boldsymbol{x}_{1}-\boldsymbol{x}_{2}$
. Due to statistical horizontal homogeneity the covariances do not depend on the horizontal sum coordinate,
$\bar{x}$
, but are permitted to depend on the vertical sum coordinate,
$\bar{z}$
, as emergent layers break the homogeneity of turbulence in the vertical.
The S3T equations for the perturbation covariances are
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}Z+\unicode[STIX]{x1D6EC}\tilde{z}\unicode[STIX]{x2202}_{\tilde{x}}Z=-2rZ+\unicode[STIX]{x2202}_{\tilde{x}}(\unicode[STIX]{x1D6E4}^{b}-\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}})+\unicode[STIX]{x1D700}\unicode[STIX]{x1D6EF}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}T+\unicode[STIX]{x1D6EC}\tilde{z}\unicode[STIX]{x2202}_{\tilde{x}}T+N_{0}^{2}\unicode[STIX]{x2202}_{\tilde{x}}(\unicode[STIX]{x0394}_{1}^{-1}\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}-\unicode[STIX]{x0394}_{2}^{-1}\unicode[STIX]{x1D6E4}^{b})+B_{1}^{\prime }\unicode[STIX]{x2202}_{\tilde{x}}\unicode[STIX]{x0394}_{1}^{-1}\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}-B_{2}^{\prime }\unicode[STIX]{x2202}_{\tilde{x}}\unicode[STIX]{x0394}_{2}^{-1}\unicode[STIX]{x1D6E4}^{b}=-2rT,\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}+\unicode[STIX]{x1D6EC}\tilde{z}\unicode[STIX]{x2202}_{\tilde{x}}\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}-(N_{0}^{2}+B_{2}^{\prime })[\tilde{\unicode[STIX]{x0394}}+\unicode[STIX]{x2202}_{\tilde{z}\bar{z}}^{2}+(1/4)\unicode[STIX]{x2202}_{\bar{z}\bar{z}}^{2}]\unicode[STIX]{x2202}_{\tilde{x}}\unicode[STIX]{x0394}_{1}^{-1}\unicode[STIX]{x0394}_{2}^{-1}Z=-2r\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}+\unicode[STIX]{x2202}_{\tilde{x}}T,\qquad & \displaystyle\end{eqnarray}$$
where
$B_{i}^{\prime }$
is the vertical derivative of
$B$
at
$z_{i}$
,
$\unicode[STIX]{x0394}_{i}$
is the Laplacian in the original coordinates
$\boldsymbol{x}_{i}$
,
$\tilde{\unicode[STIX]{x0394}}=\unicode[STIX]{x2202}_{\tilde{x}\tilde{x}}^{2}+\unicode[STIX]{x2202}_{\tilde{z}\tilde{z}}^{2}$
is the Laplacian in difference coordinates, and
$\unicode[STIX]{x1D6EF}$
is the covariance of
$\unicode[STIX]{x1D709}$
such that
$\langle \unicode[STIX]{x1D709}_{1}(t)\unicode[STIX]{x1D709}_{2}(t)\rangle \equiv \unicode[STIX]{x1D6FF}(t_{1}-t_{2})\unicode[STIX]{x1D6EF}(\tilde{x},\tilde{z})$
. The equation for
$\unicode[STIX]{x1D6E4}^{b}$
is obtained from (2.8) by applying
$\boldsymbol{x}_{1}\leftrightarrow \boldsymbol{x}_{2}$
, under which
$\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}\leftrightarrow \unicode[STIX]{x1D6E4}^{b}$
. The buoyancy flux is given by
$$\begin{eqnarray}\overline{w^{\prime }b^{\prime }}=(1/2)\unicode[STIX]{x2202}_{\tilde{x}}(\unicode[STIX]{x0394}_{1}^{-1}\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}-\unicode[STIX]{x0394}_{2}^{-1}\unicode[STIX]{x1D6E4}^{b})|_{\tilde{x}=\tilde{z}=0},\end{eqnarray}$$
so that the mean buoyancy evolves according to the equation
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}B=-r_{m}B-(1/2)\unicode[STIX]{x2202}_{\tilde{x}\bar{z}}^{2}(\unicode[STIX]{x0394}_{1}^{-1}\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}-\unicode[STIX]{x0394}_{2}^{-1}\unicode[STIX]{x1D6E4}^{b})|_{\tilde{x}=\tilde{z}=0}.\end{eqnarray}$$
For simplicity we omit here the diffusion terms included to ensure convergence, but these terms are included in all results. S3T is constituted by (2.6)–(2.8) together with (2.10) and is a closed, autonomous, deterministic dynamical system determining the coupled evolution of the mean state,
$B$
, together with the second-order statistics of the turbulence. For details of the derivation of S3T see Fitzgerald & Farrell (Reference Fitzgerald and Farrell2018b
).
S3T retains nonlinear interactions between the mean and the perturbations because these interactions are expressed in terms of second cumulants. Nonlinear interactions among perturbations are expressed in terms of third cumulants and are not retained (Herring Reference Herring1963). A stochastic approximation to S3T can thus be obtained by removing the bracketed Jacobian terms, which correspond to perturbation–perturbation nonlinearity, from (2.1)–(2.3). The resulting system, which we refer to as QL (for quasilinear), remains stochastic but otherwise has the same dynamical restrictions as S3T. QL can be viewed as employing a finite ensemble approximation to the infinite ensemble second cumulant used in S3T, and provides a bridge between fully nonlinear simulations (NL) and S3T.
3 S3T analysis of the Phillips instability
3.1 Flux-gradient relation
To analyse the possibility of layering via PI we first apply S3T to predict the flux-gradient relation between the buoyancy flux,
$\overline{w^{\prime }b^{\prime }}$
, and the stratification,
$N_{0}^{2}$
, characterizing statistically stationary homogeneous sheared stratified turbulence. Homogeneous turbulence is a fixed point of S3T and is the background in which layers are predicted to emerge if the flux-gradient relation satisfies the conditions required for PI. To obtain the flux-gradient relation we first use (2.6)–(2.8) to determine the covariance structure of the homogeneous turbulence fixed point, corresponding to
$\unicode[STIX]{x2202}_{t}=\unicode[STIX]{x2202}_{\bar{z}}=B=0$
. Defining
$\boldsymbol{p}=(p,q)$
and
$h^{2}=p^{2}+q^{2}$
, and using the Fourier convention
$$\begin{eqnarray}C(x,z,t)=\int \frac{\text{d}p\,\text{d}q}{(2\unicode[STIX]{x03C0})^{2}}\tilde{C}(p,q,t)\text{e}^{\text{i}\boldsymbol{p}\boldsymbol{\cdot }\tilde{\boldsymbol{x}}},\end{eqnarray}$$
where
$C$
denotes any covariance function, we obtain in spectral space
$$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D6EC}p\unicode[STIX]{x2202}_{q}\tilde{Z}_{H}=-2r\tilde{Z}_{H}+\text{i}p(\tilde{\unicode[STIX]{x1D6E4}}_{H}^{b}-\tilde{\unicode[STIX]{x1D6E4}}_{H}^{\unicode[STIX]{x1D701}})+\unicode[STIX]{x1D700}\tilde{\unicode[STIX]{x1D6EF}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D6EC}p\unicode[STIX]{x2202}_{q}\tilde{T}_{H}=\text{i}pN_{0}^{2}h^{-2}(\tilde{\unicode[STIX]{x1D6E4}}_{H}^{\unicode[STIX]{x1D701}}-\tilde{\unicode[STIX]{x1D6E4}}_{H}^{b})-2r\tilde{T}_{H}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D6EC}p\unicode[STIX]{x2202}_{q}\tilde{\unicode[STIX]{x1D6E4}}_{H}^{\unicode[STIX]{x1D701}}=-\text{i}pN_{0}^{2}h^{-2}\tilde{Z}_{H}-2r\tilde{\unicode[STIX]{x1D6E4}}_{H}^{\unicode[STIX]{x1D701}}+\text{i}p\tilde{T}_{H}, & \displaystyle\end{eqnarray}$$
in which the subscript
$H$
indicates homogeneous covariance functions. Defining
$$\begin{eqnarray}{\mathcal{L}}\equiv 2r-\unicode[STIX]{x1D6EC}p\unicode[STIX]{x2202}_{q},\quad {\mathcal{H}}\equiv {\mathcal{L}}h^{2}{\mathcal{L}}+2p^{2}N_{0}^{2},\end{eqnarray}$$
and combining (3.2)–(3.4) with (2.9) we obtain the relations
$$\begin{eqnarray}\displaystyle & \displaystyle [{\mathcal{H}}{\mathcal{L}}^{2}+2p^{2}N_{0}^{2}{\mathcal{H}}h^{-2}-4p^{4}h^{-2}N_{0}^{4}]\tilde{Z}_{H}=\unicode[STIX]{x1D700}{\mathcal{H}}{\mathcal{L}}\tilde{\unicode[STIX]{x1D6EF}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{H}}\tilde{T}_{H}=2p^{2}h^{-2}N_{0}^{4}\tilde{Z}_{H}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{w^{\prime }b^{\prime }}=-\frac{1}{8\unicode[STIX]{x03C0}^{2}N_{0}^{2}}\int \text{d}p\,\text{d}q{\mathcal{L}}\tilde{T}_{H}. & \displaystyle\end{eqnarray}$$
Sequentially inverting the left-hand side operators in (3.6)–(3.7) to obtain
$\tilde{T}_{H}$
and then evaluating (3.8) gives the buoyancy flux as a function of the stratification.
In the case of zero shear,
$\unicode[STIX]{x1D6EC}=0$
, if the excitation has the isotropic ring spectrum
$\tilde{\unicode[STIX]{x1D6EF}}=4\unicode[STIX]{x03C0}k_{e}\unicode[STIX]{x1D6FF}(h-k_{e})$
, the flux-gradient relation is given by the illustrative expression
$$\begin{eqnarray}\overline{w^{\prime }b^{\prime }}=-\frac{\unicode[STIX]{x1D700}}{2}\left(1-\sqrt{\frac{r^{2}}{r^{2}+N_{0}^{2}}}\right).\end{eqnarray}$$
S3T predicts that in this case the buoyancy flux is downgradient and strengthens monotonically as stratification increases, precluding the occurrence of PI. Note that the rate at which the stochastic excitation injects energy into the system is given by
$$\begin{eqnarray}\text{energy injection rate}=\unicode[STIX]{x1D700}\int \frac{\text{d}p\,\text{d}q}{(2\unicode[STIX]{x03C0})^{2}}\frac{\tilde{\unicode[STIX]{x1D6EF}}}{2h^{2}},\end{eqnarray}$$
so that
$\unicode[STIX]{x1D700}$
in (3.9) equals the total energy injection rate (Fitzgerald & Farrell Reference Fitzgerald and Farrell2018b
).
The excitation used in our numerical simulations (figure 1
$b$
) has the spectrum
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D6EF}}(p,q)=\mathop{\sum }_{i=1}^{M}{\mathcal{A}}(p)h^{4}\exp (-\ell _{c}^{2}q^{2}/4)[\unicode[STIX]{x1D6FF}(p-p_{i})+\unicode[STIX]{x1D6FF}(p+p_{i})],\end{eqnarray}$$
and excites with equal energies the horizontal wavenumber components
$p_{i}$
of the perturbation streamfunction with a Gaussian correlation length
$\ell _{c}$
in the vertical. Here
${\mathcal{A}}=2\ell _{c}^{3}\unicode[STIX]{x03C0}^{3/2}/[M(2+\ell _{c}^{2}p^{2})]$
, where
$M$
is the number of included horizontal wavenumbers. The flux-gradient relation for
$\unicode[STIX]{x1D6EC}=0$
when
$\unicode[STIX]{x1D6EF}$
is given by (3.11) is shown in figure 2(a) (dashed curve). The flux-gradient relation is very similar to that obtained for the case of ring excitation (3.9), and in particular is monotonic and does not permit PI.
For
$\unicode[STIX]{x1D6EC}\neq 0$
the flux-gradient relation is modified as shown in figure 2(a) (solid curve). Non-dimensionalizing so that
$r=1$
and the domain width equals
$2\unicode[STIX]{x03C0}$
, the parameters used are
$\unicode[STIX]{x1D6EC}=10$
,
$\unicode[STIX]{x1D700}=2$
,
$\unicode[STIX]{x1D708}=5\times 10^{-4}$
,
$\ell _{c}=10^{-1}$
and the excited wavenumbers are
$p_{i}\in \{5,\ldots ,10\}$
. The flux-gradient relation is non-monotonic and has positive slope over a range of stratification values, suggesting that PI may occur. In § 3.2 we next apply S3T to directly analyse the perturbation stability of homogeneous turbulence.
The mechanism responsible for modifying the flux-gradient relation in the presence of shear can be understood using the equilibrium turbulent kinetic energy (TKE) budget
$$\begin{eqnarray}0=\unicode[STIX]{x1D700}+\overline{w^{\prime }b^{\prime }}-\unicode[STIX]{x1D6EC}\overline{u^{\prime }w^{\prime }}+\text{dissipation}.\end{eqnarray}$$
The dotted curve in figure 2(a) shows the shear production of TKE,
$-\unicode[STIX]{x1D6EC}\overline{u^{\prime }w^{\prime }}$
, as a function of
$N_{0}^{2}$
. Shear production of TKE results from transient growth of 2D gravity waves in shear, the dynamics of which has previously been analysed by Farrell & Ioannou (Reference Farrell and Ioannou1993). Farrell & Ioannou (Reference Farrell and Ioannou1993) showed that transient growth and shear production are associated with enhanced buoyancy fluxes, consistent with the approximate covariation of shear production and buoyancy flux in figure 2(a), and that transient growth is suppressed as stratification increases. These properties of gravity wave shear dynamics suggest the following explanation for the structure of the modified flux-gradient relation in figure 2(a). For very weak stratification (
$N_{0}^{2}\approx 0$
), the buoyancy flux is weak due to the absence of significant buoyancy perturbations. For non-zero stratification that is weak relative to the shear (i.e.
$N_{0}^{2}>0$
but
$N_{0}^{2}\ll \unicode[STIX]{x1D6EC}^{2}$
), shear production of TKE is associated with buoyancy fluxes that are strongly enhanced relative to the unsheared case. As the stratification is further increased (
$N_{0}^{2}\gtrsim \unicode[STIX]{x1D6EC}^{2}$
) the shear production, and correspondingly the buoyancy flux, weakens, resulting in a region of positive slope in the buoyancy flux-gradient relation.
Figure 2. (a) Flux-gradient relation predicted by S3T without shear (
$\unicode[STIX]{x1D6EC}=0$
, dashed) and with shear (
$\unicode[STIX]{x1D6EC}=10$
, solid). For
$\unicode[STIX]{x1D6EC}=0$
the flux varies monotonically with
$N_{0}^{2}$
while for
$\unicode[STIX]{x1D6EC}=10$
the flux varies non-monotonically with
$N_{0}^{2}$
, with a positive slope as required for PI indicated by the tangent line occurring for the example point
$N_{0}^{2}=100$
. The dotted curve shows the shear production of TKE,
$-\unicode[STIX]{x1D6EC}\overline{u^{\prime }w^{\prime }}$
, for
$\unicode[STIX]{x1D6EC}=10$
. Weakening of shear production for
$N_{0}^{2}\gtrsim 30$
is associated with weakening buoyancy fluxes and non-monotonicity of the flux-gradient relation. (b) Layer growth rate,
$s$
, predicted by S3T as a function of the layer vertical wavenumber,
$m$
. Five parameter cases are shown (colours correspond to the parameter values indicated in a). Layers are predicted to form in the case
$\unicode[STIX]{x1D6EC}=10$
,
$N_{0}^{2}=100$
, consistent with the positive slope of the flux-gradient relation in (a). The dashed black curve shows the approximate growth rate from (1.4), modified to account for explicit dissipation. Dashed and dotted red curves show the case
$\unicode[STIX]{x1D6EC}=10$
,
$N_{0}^{2}=100$
, but with the excitation spectrum (3.11) modified to only excite
$p=5$
(dashed) and
$p=10$
(dotted). Parameter details can be found in § 3.
3.2 Linear stability analysis
To analyse the stability of homogeneous turbulence we linearize (2.6)–(2.8) and (2.10) about the S3T fixed point
$B=0$
,
$C=C_{H}$
. Because homogeneous covariance functions do not depend on
$\bar{z}$
, we analyse perturbations using harmonic functions in
$\bar{z}$
. We write the mean buoyancy perturbation and the perturbations to the covariance functions as
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}B=\text{e}^{st}\text{e}^{\text{i}m\bar{z}},\quad \unicode[STIX]{x1D6FF}C=\text{e}^{st}\text{e}^{\text{i}m\bar{z}}\int \frac{\text{d}p\,\text{d}q}{(2\unicode[STIX]{x03C0})^{2}}\unicode[STIX]{x1D6FF}\tilde{C}_{m}(p,q)\text{e}^{\text{i}\boldsymbol{p}\boldsymbol{\cdot }\tilde{\boldsymbol{x}}},\end{eqnarray}$$
plus complex conjugate terms, where
$m$
is the vertical wavenumber of the layers and
$s$
is the growth rate. Substituting (3.13) into the linearized S3T dynamics we obtain
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\begin{array}{@{}cccc@{}}{\mathcal{L}}^{\prime } & 0 & \text{i}p & -\text{i}p\\ 0 & {\mathcal{L}}^{\prime } & -\text{i}pN_{0}^{2}h_{+}^{-2} & \text{i}pN_{0}^{2}h_{-}^{-2}\\ \text{i}pN_{0}^{2}h_{-}^{-2} & -\text{i}p & {\mathcal{L}}^{\prime } & 0\\ -\text{i}pN_{0}^{2}h_{+}^{-2} & \text{i}p & 0 & {\mathcal{L}}^{\prime }\end{array}\right)\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D6FF}\tilde{Z}\\ \unicode[STIX]{x1D6FF}\tilde{T}\\ \unicode[STIX]{x1D6FF}\tilde{\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}}\\ \unicode[STIX]{x1D6FF}\tilde{\unicode[STIX]{x1D6E4}^{b}}\end{array}\right)=mp\left(\begin{array}{@{}c@{}}0\\ -h_{-}^{-2}\tilde{\unicode[STIX]{x1D6E4}}_{H,-}^{\unicode[STIX]{x1D701}}+h_{+}^{-2}\tilde{\unicode[STIX]{x1D6E4}}_{H,+}^{b}\\ h_{+}^{-2}\tilde{Z}_{H,+}\\ -h_{-}^{-2}\tilde{Z}_{H,-}\end{array}\right), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle s=-r_{m}-\frac{m}{2}\int \frac{\text{d}p\,\text{d}q}{(2\unicode[STIX]{x03C0})^{2}}p(h_{+}^{-2}\unicode[STIX]{x1D6FF}\tilde{\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701}}}-h_{-}^{-2}\unicode[STIX]{x1D6FF}\tilde{\unicode[STIX]{x1D6E4}^{b}}), & \displaystyle\end{eqnarray}$$
where we suppress dependencies on
$(p,q,m)$
and define
${\mathcal{L}}^{\prime }=s+2r-\unicode[STIX]{x1D6EC}p\unicode[STIX]{x2202}_{q}$
,
$h_{\pm }^{2}=p^{2}+(q\pm m/2)^{2}$
and
$\tilde{C}_{H,\pm }(p,q)=\tilde{C}_{H}(p,q\pm m/2)$
. Equations (3.14)–(3.15) constitute an eigenvalue problem for
$s$
. We have chosen, without loss of generality, to set the amplitude of
$\unicode[STIX]{x1D6FF}B$
to 1 in (3.13), which provides a normalization of the eigenfunctions. Using an initial guess for
$s$
, (3.14) can be inverted to obtain the perturbation covariance structures including
$\tilde{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D701},b}}$
. An updated guess for
$s$
can then be obtained from (3.15) and the eigenproblem can be solved by iteration.
Figure 2(b) shows
$s$
as a function of
$m$
in cases with and without shear and in regions of positive and negative slope of the flux-gradient relation. For negative slope,
$s$
is negative for all
$m$
. For the positive slope example
$\unicode[STIX]{x1D6EC}=10$
,
$N_{0}^{2}=100$
, which is our primary example case, S3T predicts that
$s>0$
for
$3\leqslant m\leqslant 15$
(bright red, solid curve). The black dashed curve in figure 2(b) shows the prediction of
$s$
based directly on the flux-gradient relation using (1.4), modified to account for explicit mean damping,
$r_{m}$
. This approximation agrees with the S3T prediction for small
$m$
, when the layers are of large scale relative to the turbulence, and begins to deviate as
$m$
increases. (We note that for
$m=0$
, the mean buoyancy perturbation is a perturbative shift of
$B$
with no vertical structure. A uniform perturbation in
$B$
does not induce a perturbation in the turbulent buoyancy flux divergence, and so (2.1) implies that
$s(m=0)=-r_{m}$
, as shown in figure 2(b).) For strong stratification, the flux-gradient relation becomes flat and
$s$
becomes negative. That predictions based on the flux-gradient relation correspond to S3T predictions indicates that the mechanism underlying the S3T layering instability is PI. In § 4 we show that layering is observed in NL simulations, consistent with S3T predictions.
We estimate the conventional characteristic scales and non-dimensional parameters for our primary example case (
$N_{0}^{2}=100$
,
$\unicode[STIX]{x1D6EC}=10$
) as follows, using the homogeneous turbulent equilibrium solution of (3.2)–(3.4). The Richardson number is
$Ri=\unicode[STIX]{x1D6EC}^{2}/N_{0}^{2}=1$
. The characteristic velocity scale is estimated using the square root of the TKE,
$u_{RMS}^{\prime }\equiv \sqrt{2K^{\prime }}\approx 1.1$
. The characteristic wavenumber of the turbulence is estimated using the centre of energy of the homogeneous turbulence spectrum,
$1/L_{0}\approx 19.6$
. The Froude number is estimated as
$Fr\equiv u_{RMS}^{\prime }/(L_{0}N_{0})\approx 2.1$
. We note that we have also obtained examples in which layering is predicted by S3T for
$Fr<1$
. The buoyancy Reynolds number is estimated as
$Re_{b}\equiv \unicode[STIX]{x1D700}/(\unicode[STIX]{x1D708}N_{0}^{2})=40$
. Estimating
$Re_{b}$
based on nonlinear numerical simulations (see § 4) and accounting for the effects of linear damping gives the reduced estimate
$Re_{b}\approx 12.2$
. The Ozmidov and buoyancy wavenumbers are estimated as
$k_{O}\equiv (N_{0}^{3}/\unicode[STIX]{x1D700})\approx 22.4$
and
$k_{b}\equiv N_{0}/u_{RMS}^{\prime }\approx 9.4$
.
Although
$k_{O}$
and
$k_{b}$
are comparable to the small-scale cutoff wavenumber of the layering instability for the primary example case in figure 2(b) (
$m=15$
), we have found no clear relationship between these characteristic scales and the cutoff wavenumber. As
$N_{0}^{2}$
increases, PI is suppressed due to flattening of the flux-gradient curve, whereas
$k_{O}$
and
$k_{b}$
simply move to smaller scale. The cutoff wavenumber of PI does, however, depend in a natural way on the spectrum of the underlying turbulence, as controlled in our model configuration by the excitation spectrum. The dashed and dotted red curves in figure 2(b) show the layer growth rates when the excitation spectrum (3.11), which excites the horizontal wavenumber components
$p\in \{5,\ldots ,10\}$
, is replaced by a spectrum exciting only
$p=5$
(dashed) or only
$p=10$
(dotted), with the correlation length and total energy injection rate held constant. Growth rates agree across the three excitation structures at the largest scales (
$m\rightarrow 0$
), where scale separation between the turbulence and the layers is robust. As
$m$
is increased, larger-scale excitation (
$p=5$
, dashed red curve) deviates rapidly from the flux-gradient prediction (dashed black curve) as scale separation is lost. Smaller-scale excitation (
$p=10$
, dotted red curve) produces a background turbulence that enables PI down to much smaller scales.
The dissipation parameters of the model influence PI in a manner that can be simply understood. Numerical experimentation (not shown) reveals that the primary dependence of the growth rate curve,
$s(m)$
, on the mean relaxation rate,
$r_{m}$
, is that
$s$
shifts towards more negative values as
$r_{m}$
is increased (stronger relaxation) and towards more positive values as
$r_{m}$
is decreased (weaker relaxation). This behaviour reflects the manner in which
$r_{m}$
appears in the dispersion relation (3.15), which is as a shift of the growth rate independent of the other parameters. For sufficiently large values of
$r_{m}$
, PI is suppressed so that
$s(m)<0$
for all
$m$
. This result is expected because rapid relaxation of the stratification towards the uniform background profile prevents the development of perturbations to the stratification. We also note that PI exhibits a similar dependence on the perturbation damping rate,
$r$
. Increasing
$r$
weakens the PI growth rate because the turbulence responsible for the buoyancy fluxes is suppressed by the enhanced dissipation. Reducing
$r$
enhances the shear production of TKE and increases the positive slope of the flux-gradient relation, thereby enhancing the tendency of the system to undergo PI. Diffusive dissipation, with coefficient
$\unicode[STIX]{x1D708}$
, is included only for numerical stability.
4 Emergence of layers in simulations
The evolution of the mean buoyancy,
$B$
, is shown in figure 3(a) in an NL simulation using parameters
$\unicode[STIX]{x1D6EC}=10$
,
$N_{0}^{2}=100$
corresponding to the primary example case shown in figure 2 (bright red, solid curve) for which S3T predicts that PI will occur. The model configuration and parameters are as in §§ 2 and 3. The NL system is integrated using DIABLO (Taylor Reference Taylor2008) with
$256\times 512$
resolution. Approximately four layers emerge and persist for the duration of the simulation. Figure 3(d) shows the time and horizontal mean stratification,
$[N^{2}]=N_{0}^{2}+\unicode[STIX]{x2202}[b]/\unicode[STIX]{x2202}z$
, where square brackets indicate the time and horizontal mean and the time average is calculated over the final 10 time units of the simulation. As described by (2.1), the layers are formed and maintained by turbulent buoyancy fluxes which are balanced at equilibrium by explicit relaxation towards the uniform stratification profile.
Figure 3. Evolution of the mean buoyancy anomaly,
$B$
, in example simulations of the (a) NL, (b) QL and (c) S3T systems. (d) Profiles of
$[N^{2}-N_{0}^{2}]$
, the deviation of the time and horizontal mean stratification from the background stratification, in NL, QL and S3T. (e) Potential energy of the time and horizontal mean buoyancy field,
$[V]$
, in NL, QL and S3T as a function of
$N_{0}^{2}$
. The increase in
$[V]$
as
$N_{0}^{2}$
is increased beyond
$N_{0}^{2}\approx 25$
, and subsequent decrease as
$N_{0}^{2}$
is increased beyond
$N_{0}^{2}\approx 120$
, in NL and QL reflect the bifurcations occurring in S3T at these values of
$N_{0}^{2}$
. (e) Maximum S3T eigenvalue,
$\text{max}_{m}[s(m)]$
, as a function of
$N_{0}^{2}$
. S3T predicts that PI occurs in the shaded region, and that bifurcations occur at the edges of the shaded region. The parameters used are
$\unicode[STIX]{x1D6EC}=10$
and
$N_{0}^{2}=100$
.
Layer evolution in QL and S3T is shown in figure 3(b,c), and stratification profiles are shown in figure 3(d). The time average profile is shown for QL while for S3T the profile is shown after a fixed point has been reached. QL is integrated using a mixed Galerkin finite difference method and S3T is integrated using the matrix formulation described in Fitzgerald & Farrell (Reference Fitzgerald and Farrell2018a
), each with 256 grid points in the vertical direction. In the matrix formulation of S3T the perturbation fields are described by their covariance matrices, which are time-evolved using the time-dependent Lyapunov equation. Identical parameter values are used in NL, QL and S3T. As was found in NL, layers emerge and equilibrate in QL and S3T, with the QL and S3T layers having
${\sim}3$
–
$4$
times larger amplitude and somewhat smaller scale than the NL layers. In S3T, the wavenumber of the emergent layers reflects the
$m$
value at which
$s$
is maximized in figure 2(b). Because S3T is deterministic, this pattern emerges in S3T for any sufficiently small initial perturbation, while the NL and QL patterns differ somewhat between realizations.
The dependence of PI on
$N_{0}^{2}$
is shown in figure 3(e). Simulations of NL, QL and S3T were carried out over the range
$20\leqslant N_{0}^{2}\leqslant 300$
, with the other parameters held fixed. The background state is stable to hydrodynamic modal instability in all cases. To diagnose layering we calculate the time and horizontal mean potential energy,
$[V]=[b]^{2}/(2N_{0}^{2})$
, taking the average over the final 10 time units in NL and QL and using the final state for S3T. As the layers in NL are of smaller amplitude than the QL and S3T layers, the NL value of
$[V]$
is shown multiplied by ten for clarity. We note that the differences in
$[V]$
values between the QL and S3T simulations, which have the same underlying dynamics, are primarily due to the emergence of larger-scale layers in the QL system relative to the S3T system, as illustrated in figure 3(d). S3T equilibria with
$[V]$
values similar to those of the QL system are obtained when the S3T system is initialized with a perturbative layering pattern reflecting that found in the QL system.
In each system
$[V]$
increases rapidly as
$N_{0}^{2}$
is increased and subsequently decreases rapidly for sufficiently large
$N_{0}^{2}$
. This behaviour is consistent with the S3T prediction summarized in figure 3(f), which shows
$\text{max}_{m}[s(m)]$
as a function of
$N_{0}^{2}$
. For weak stratification,
$s<0$
consistent with the negative slope of the flux-gradient relation in figure 2(a). As
$N_{0}^{2}$
is increased into the positively sloping region,
$s$
becomes positive. When
$N_{0}^{2}$
is increased further, the slope of the flux-gradient relation remains positive but becomes very small and
$s$
falls below zero. In S3T the appearance and disappearance of layers as
$N_{0}^{2}$
increases occur through bifurcations which are reflected in NL and QL simulations by rapid changes in
$[V]$
.
The qualitative consistency of the layering behaviour in NL with that in QL and S3T, in which perturbation–perturbation interactions are not retained in the equations of motion, indicates that the dynamics underlying PI in this model of stratified turbulence is the interaction, in the presence of shear, between the mean buoyancy anomaly and the perturbations. Nonlinear perturbation–perturbation advection is revealed to be inessential to the dynamics responsible for layer formation, as layers form robustly when this nonlinearity is not included. We note that this situation differs from what would be expected if layers formed via an upscale cascade of energy from small-scale turbulence to larger-scale layers, in which case perturbation–perturbation nonlinear interactions would be essential to the spectrally local cascade and layers would not form in QL and S3T simulations. However, we emphasize that nonlinearity does influence layering quantitatively, as evidenced by the significant quantitative differences between the NL layers and those in QL and S3T shown in figure 3. Deviation of the NL layers from the structure predicted by S3T likely arises from the modification in NL of the turbulence spectrum by perturbation–perturbation nonlinearity, which is reflected in the associated flux-gradient relation.
5 Discussion
Layering in stratified turbulence via PI has remained enigmatic since its original prediction. One contribution of this work is to demonstrate PI unambiguously in 2D sheared stratified Boussinesq turbulence using S3T. This identification of PI is made in the context of our model of turbulence, which we emphasize includes dynamical simplifications such as linear relaxation and fixed background shear. Another contribution of the present analysis is to complement and extend previous approaches to the analysis of PI based on phenomenological turbulence models. For example, Phillips (Reference Phillips1972) and Posmentier (Reference Posmentier1977) identified that the original formulation of PI was ill-posed due to the existence of an ultraviolet catastrophe. Barenblatt et al. (Reference Barenblatt, Bertsch, Dal Passo, Prostokishin and Ughi1993) and Balmforth, Llewellyn Smith & Young (Reference Balmforth, Llewellyn Smith and Young1998) showed that under more realistic phenomenological assumptions the ultraviolet catastrophe is resolved. Wunsch & Kerstein (Reference Wunsch and Kerstein2001) used a stochastic phenomenological model of turbulent advection to argue that layer formation can result from the action of incoherent turbulent motions at small scales. S3T analysis verifies directly that PI occurs and is physically and mathematically consistent in the turbulence dynamics, as previously predicted using phenomenological arguments.
The present work also places PI into the context of other fundamentally statistical instabilities of turbulence. S3T has previously been used to show that zonal jets in barotropic
$\unicode[STIX]{x1D6FD}$
-plane turbulence (Farrell & Ioannou Reference Farrell and Ioannou2003) and vertically sheared horizontal flows in stratified turbulence (Fitzgerald & Farrell Reference Fitzgerald and Farrell2018a
) also form via S3T instabilities of homogeneous turbulence. These instabilities occur via interaction between the mean state and quadratic turbulent fluxes and have analytical expression only in SSD. Equilibrium statistical mechanics has recently been applied to investigate mixing in inviscid, unforced stratified turbulence (Venaille, Gostiaux & Sommeria Reference Venaille, Gostiaux and Sommeria2017). Venaille et al. (Reference Venaille, Gostiaux and Sommeria2017) found evidence of a non-monotonic relationship between mixing efficiency and Richardson number, consistent with the necessary criterion for the occurrence of the Phillips mechanism. However, the connection between equilibrium statistical mechanics and the behaviour of forced-dissipated systems awaits further clarification.
Our analysis was carried out using an idealized model of stratified turbulence. The configuration was chosen, in the spirit the recent work of Taylor & Zhou (Reference Taylor and Zhou2017), to provide a simple example in which the dynamics can be analysed with clarity using S3T, rather than to faithfully represent geophysical or laboratory turbulence. S3T analysis predicts that, in this model system, PI requires the presence of shear. Layers have been observed in laboratory studies of rod-stirred stratified turbulence without any coherent large-scale shear (Ruddick et al. Reference Ruddick, McDougall and Turner1989; Park et al. Reference Park, Whitehead and Gnandadesikan1994; Holford & Linden Reference Holford and Linden1999), but layer formation in these systems has recently been attributed to a mechanism differing from PI (Thorpe Reference Thorpe2016). The present work provides a theoretical underpinning of PI based on SSD, but is not designed to capture alternative mechanisms that may be important in the case of rod-stirred tanks. That a background shear is required to enable layer formation in our system is consistent with the results of previous studies in which layering occurred in the presence of background shear flows (Manucharyan & Caulfield Reference Manucharyan and Caulfield2015; Lucas et al. Reference Lucas, Caulfield and Kerswell2017; Taylor & Zhou Reference Taylor and Zhou2017; Zhou et al. Reference Zhou, Taylor, Caulfield and Linden2017). We note that although vertical shear is used in our study, Lucas et al. (Reference Lucas, Caulfield and Kerswell2017) observed layer formation in simulations in which a body force was used to drive horizontal shear. In the simulations of Lucas et al. (Reference Lucas, Caulfield and Kerswell2017) the initially horizontal shear spontaneously developed vertical structure and the turbulent equilibrium state consisted of a horizontally and vertically sheared mean flow coexisting with density layers. Examination of the relationship between the background shear and buoyancy layer formation, and other predictions of S3T, in more realistic stratified turbulence is an important avenue of future investigation.
Acknowledgements
The authors thank J. R. Taylor for providing the DIABLO code and three anonymous reviewers whose comments helped to improve the manuscript. J.G.F. was partially supported by a doctoral fellowship from the National Sciences and Engineering Research Council of Canada. B.F.F. was partially supported by the U.S. National Science Foundation under grant nos. NSF AGS-1246929 and NSF AGS-1640989.
