# Analysis of scale-dependent kinetic and potential energy in sheared, stably
stratified turbulence

Xiaolong Zhang  Rohit Dhariwal  Gavin Portwood  Stephen M. de Bruyn Kops
Andrew D. Bragg  Xiaolong Zhang1  Rohit Dhariwal2  Gavin Portwood3  Stephen M.
de Bruyn Kops 4  Andrew D. Bragg1<EMAIL_ADDRESS>1Department of Civil
and Environmental Engineering, Duke University, Durham, NC 27708, USA 2Center
for Institutional Research Computing, Washington State University, Pullman, WA
99164, USA 3Lawrence Livermore National Laboratory, Livermore, CA 94550, USA
4Department of Mechanical and Industrial Engineering, University of
Massachusetts Amherst, Amherst, MA 01003, USA

###### Abstract

Budgets of turbulent kinetic energy (TKE) and turbulent potential energy (TPE)
at different scales $\ell$ in sheared, stably stratified turbulence are
analyzed using a filtering approach. Competing effects in the flow are
considered, along with the physical mechanisms governing the energy fluxes
between scales, and the budgets are used to analyze data from direct numerical
simulation (DNS) at buoyancy Reynolds number $Re_{b}=O(100)$. The mean TKE
exceeds TPE by an order of magnitude at the large scales, with the difference
reducing as $\ell$ is decreased. At larger scales, buoyancy is never observed
to be positive, with buoyancy always converting TKE to TPE. As $\ell$ is
decreased, the probability of locally convecting regions increases, though it
remains small at scales down to the Ozmidov scale. The TKE and TPE fluxes
between scales are both downscale on average and their instantaneous values
are positively correlated, but not strongly so, and this occurs due to the
different physical mechanisms that govern these fluxes. Moreover, the
contribution to these fluxes arising from the sub-grid fields are shown to be
significant, in addition to the filtered scale contributions associated with
the processes of strain-self amplification, vortex stretching, and density
gradient amplification. Probability density functions (PDFs) of the $Q,R$
invariants of the filtered velocity gradient are considered. Unlike isotropic
turbulence, as $\ell$ increases the sheared-drop shape of the PDF disappears
and the PDF becomes symmetric about $R=0$, meaning regions of vortex
stretching and compression become equi-probable, as well as regions of strain
amplification or suppression.

##  1 Introduction

When turbulence occurs in environmental flows, it is often affected by both
mean-shear and stable stratification (Vallis, 2006; Wyngaard, 2010; Ferrari &
Wunsch, 2009; Zorzetto et al., 2018; Ayet et al., 2020), leading to Sheared,
Stably Stratified Turbulence (SSST). Mean-shear (vertical gradient of
horizontal flow) produces turbulence, with hairpin vortices, streaks, and
strong fluctuations in all three directions (Lee et al., 1990; Davidson, 2004;
Pope, 2000). Stable stratification suppresses fluctuations in the vertical
direction, and if sufficiently strong generates a quasi two-dimensional flow
behavior (Riley & Lelong, 2000; Riley & Lindborg, 2012). Stable stratification
also provides a restoring buoyancy force that enables the propagation of
internal waves and the formation of quasi-horizontal “pancake” vortical
structures in the flow (Davidson, 2004). Since mean-shear and stable
stratification have competing effects, and produce different flow structures,
the dynamics of SSST are rich and complex. Moreover, the flow is often
“patchy”, with turbulent and non-turbulent regions interspersed, depending on
the local competition between shear and buoyancy. Indeed, when the local shear
in the flow is strong enough for the local Richardson number to be
sufficiently small, the local flow may undergo Kelvin-Helmholtz instabilities
which can evolve due to nonlinearity into turbulent motion (Riley & de Bruyn
Kops, 2003).

Understanding and modeling SSST is an active area of research with many open
questions. One such issue concerns understanding the mixing efficiencies in
SSST and their parametric dependence (Peltier & Caulfield, 2003; Portwood et
al., 2019), which are vital for predicting mixing in oceans (Jayne, 2009;
Gregg et al., 2018). Another vital area is to understand the properties of
SSST across its range of dynamical scales, and the physical mechanisms that
govern the fluxes of turbulent kinetic and potential energy between scales.
Not only is this important for a basic understanding of the flow, it is also
of crucial importance for developing large eddy simulation (LES) models for
SSST. A number of studies have considered the effect of stable stratification
on the multiscale properties of turbulence using Fourier analysis and
considering the average behavior of the flow in terms of the energy spectrum,
the mean buoyancy and mean inter-scale energy transfer terms (Riley & de Bruyn
Kops, 2003; Lindborg, 2006; Almalkie & de Bruyn Kops, 2012). The study of
Riley & de Bruyn Kops (2003) showed that the horizonal energy spectrum
exhibits a $k_{h}^{-5/3}$ scaling for wavenumbers smaller than the Ozmidov
wavenumber $k_{O}$ (corresponding to the wavenumber at which inertial and
buoyancy forces are of the same order), where $k_{h}$ is the horizontal
wavenumber, and the results indicated a downscale energy transfer of kinetic
energy in the flow. This motivated Lindborg (2006) who confirmed that strongly
stratified turbulent flows exhibit a downscale cascade of turbulent kinetic
and potential energy on average, and developed phenomenological predictions
similar in spirit to Kolmogorov’s 1941 theory (Kolmogorov, 1941). The
observation of a downscale energy cascade was contrary to predictions that had
been made in Gage (1979); Lilly (1983) based on the assumption that in the
limit of strong stratification, the flow should behave as two-dimensional
turbulence. The basic reason why the predictions of Gage (1979); Lilly (1983)
failed is that, as shown in Billant & Chomaz (2001), as the strength of the
stratification increases and the vertical velocity of the flow is suppressed,
the vertical lengthscale of the flow also reduces in such a way that the terms
in the dynamical equations associated with vertical motion always remain
$O(1)$, and hence the flow never becomes two-dimensional.

Many questions remain, however, regarding the multiscale properties of SSST.
For example, how do fluctuations of the SSST flow about its mean-field state
behave? What are the mechanisms of the turbulent kinetic and potential energy
transfers among scales, and to what extent are these transfers coupled to each
other, and to fluctuations in the local buoyancy? How do the sub-grid scale
terms in SSST contribute to the flow dynamics across scales? To address these
and other questions, in this paper we will use a filtering-based approach
wherein the velocity and density fields are considered for different filtering
lengths $\ell$. Using these we explore the turbulent kinetic and potential
energies at different scales, and the various processes that determine their
behavior. Using a filtering approach allows us to consider spatially local
couplings between the processes that control the flow energetics at different
scales, and to consider fluctuations in the flow as well as the mean-field
behavior. Moreover, such an analysis of the multiscale properties of SSST can
provide insights for developing large eddy simulation (LES) models of SSST,
since LES models are often developed in the physical-space, rather than
Fourier space context. LES modeling of stratified turbulence is an under-
developed area Khani & Waite (2015), and new insights into the dynamics of
stratified turbulence across scales could would aid in the development of
appropriate sub-grid models.

The outline of the paper is as follows. In §2, we introduce the equations
governing the turbulent kinetic and potential energy at different scales,
consider the behavior of the mean-field state of the flow at different scales,
and discuss the physical mechanisms governing the energy transfer among scales
in SSST. In §3 we describe the DNS dataset used in this study, and in §4 we
present and discuss results from the DNS for the quantities introduced in §2
that govern the turbulent kinetic and potential energy at different scales.
Finally, in §5 we draw conclusions from the study and discuss future areas for
investigation.

##  2 Theoretical Considerations

###  2.1 Governing equations for scale-dependent energy fields

We consider the case where the mean velocity gradient, $\gamma$, and the mean
density gradient, $\zeta$, are constant in space and time, such that the total
fluid velocity vector can be written as
$\boldsymbol{U}=z\gamma\boldsymbol{e}_{x}+\boldsymbol{u}$, where
$\boldsymbol{u}$ is the fluctuating component of the velocity, and the total
density as $\rho=\rho_{r}+z\zeta+\rho^{\prime}$, where $\rho_{r}$ is a
constant reference density and $\rho^{\prime}$ is the fluctuating component of
the density. While $\boldsymbol{e}_{x}$ is the unit vector in the direction of
the mean velocity $\langle\boldsymbol{U}\rangle=z\gamma\boldsymbol{e}_{x}$,
$\boldsymbol{e}_{z}$ is the unit vector pointing in the direction opposite to
the gravitational acceleration $\boldsymbol{g}$ (i.e. the vertical direction),
and $\boldsymbol{e}_{y}\equiv\boldsymbol{e}_{z}\times\boldsymbol{e}_{x}$.

The filtering operator to be used is defined for an arbitrary field
$\boldsymbol{a}(\boldsymbol{x},t)$ as
$\widetilde{\boldsymbol{a}}(\boldsymbol{x},t)\equiv\int_{\mathbb{R}^{3}}\boldsymbol{a}(\boldsymbol{x}-\boldsymbol{y},t)\mathcal{G}_{\ell}(\boldsymbol{y})\,d\boldsymbol{y}$,
with $\mathcal{G}_{\ell}$ an isotropic filtering kernel with lengthscale
$\ell$. In order for the energy fields introduced below to be strictly non-
negative, which is a physical requirement, then a necessary and sufficient
condition is that the filtering kernel $\mathcal{G}_{\ell}$ be non-negative
(Vreman et al., 1994). We will use the isotropic Gaussian filter
$\mathcal{G}_{\ell}(\boldsymbol{y})\equiv(2\pi\ell^{2})^{-3/2}\exp(-\|\boldsymbol{y}\|^{2}/2\ell^{2})$
which satisfies this condition.

Assuming that $|\rho^{\prime}|/\rho\ll 1$, the governing equation for
$\boldsymbol{u}$ is then the Boussinesq-Navier-Stokes equation coupled with an
advection-diffusion equation for $\rho^{\prime}$ (Vallis, 2006). The filtered
versions of these equations are

$\displaystyle D_{t}\widetilde{\boldsymbol{u}}$
$\displaystyle=-(1/\rho_{r})\boldsymbol{\nabla}\widetilde{p}-\boldsymbol{\nabla\cdot\tau}+2\nu\boldsymbol{\nabla\cdot}\widetilde{\boldsymbol{s}}-N\widetilde{\phi}\boldsymbol{e}_{z}+\widetilde{\boldsymbol{F}},$
(1) $\displaystyle D_{t}\widetilde{\phi}$
$\displaystyle=-\boldsymbol{\nabla\cdot\Sigma}+\kappa\nabla^{2}\widetilde{\phi}+N\widetilde{u}_{z}+\widetilde{f},$
(2)

where
${D_{t}\equiv\partial_{t}+\widetilde{\boldsymbol{u}}\boldsymbol{\cdot\nabla}}$,
$\widetilde{\phi}\equiv g\widetilde{\rho^{\prime}}/N\rho_{r}$ is the scaled
density field (with dimensions of a velocity), $g\equiv\|\boldsymbol{g}\|$,
$\nu$ and $\kappa$ are the kinematic viscosity and thermal diffusivity,
respectively, $N\equiv\sqrt{-g\zeta/\rho_{r}}$ is the buoyancy frequency,
$\boldsymbol{\tau}\equiv\widetilde{\boldsymbol{uu}}-\widetilde{\boldsymbol{u}}\widetilde{\boldsymbol{u}}$
is the sub-grid stress tensor,
$\widetilde{\boldsymbol{s}}\equiv(\boldsymbol{\nabla}\widetilde{\boldsymbol{u}}+[\boldsymbol{\nabla}\widetilde{\boldsymbol{u}}]^{\top})/2$
is the filtered strain-rate, and
$\boldsymbol{\Sigma}\equiv\widetilde{\boldsymbol{u}\phi}-\widetilde{\boldsymbol{u}}\widetilde{\phi}$
is the sub-grid vector for the scaled density field. The terms
$\boldsymbol{F}$ and $f$, whose filtered forms appear in the equations above,
are the forcing terms associated with mean velocity gradient applied to the
flow, and for the SSST under consideration are given by

$\displaystyle\boldsymbol{F}$
$\displaystyle=-\gamma(\boldsymbol{e}_{z}\boldsymbol{\cdot
x})(\boldsymbol{e}_{x}\boldsymbol{\cdot\nabla})\boldsymbol{u}-\gamma\boldsymbol{e}_{x}(\boldsymbol{e}_{z}\boldsymbol{\cdot
u}),$ (3) $\displaystyle f$
$\displaystyle=-\gamma(\boldsymbol{e}_{z}\boldsymbol{\cdot
x})(\boldsymbol{e}_{x}\boldsymbol{\cdot\nabla})\phi.$ (4)

Following Germano (1992), the turbulent kinetic energy (TKE) in a given region
of size $\ell$ is $\widetilde{\boldsymbol{u\cdot u}}/2$ and this may be split
up into the TKE at scales $\geq\ell$ (“large-scales”) denoted by
$E_{K}\equiv\|\widetilde{\boldsymbol{u}}\|^{2}/2$ and the TKE at scales
$<\ell$ (“small-scales”) denoted by
$e_{K}\equiv(\widetilde{\|\boldsymbol{u}\|^{2}}-\|\widetilde{\boldsymbol{u}}\|^{2})/2=\mathrm{tr}[\boldsymbol{\tau}]/2$.
The equations governing $E_{K}$ and $e_{K}$ may be derived from (1) and are
given by

$\displaystyle D_{t}E_{K}$
$\displaystyle=-(1/\rho_{r})\boldsymbol{\nabla\cdot}(\widetilde{\boldsymbol{u}}\widetilde{p})-\boldsymbol{\nabla\cdot}(\widetilde{\boldsymbol{u}}\boldsymbol{\cdot\tau})-\Pi_{K}+2\nu\boldsymbol{\nabla\cdot}(\widetilde{\boldsymbol{u}}\boldsymbol{\cdot}\widetilde{\boldsymbol{s}})-2\nu\|\widetilde{\boldsymbol{s}}\|^{2}-N\widetilde{{u}}_{z}\widetilde{\phi}+\widetilde{\boldsymbol{F}}\boldsymbol{\cdot}\widetilde{\boldsymbol{u}},$
(5)
$\displaystyle\begin{split}D_{t}e_{K}&=\boldsymbol{\nabla\cdot}\boldsymbol{T}_{K}+\mathcal{B}+\Pi_{K}-\varepsilon_{K}+\mathcal{F}_{K}.\end{split}$
(6)

In these equations, the scale-to-scale TKE flux is defined as
${\Pi_{K}\equiv-\boldsymbol{\tau:}\widetilde{\boldsymbol{s}}}$, such that
$\Pi_{K}>0$ corresponds to a transfer of TKE from the large to the small
scales. The small-scale buoyancy term is
$\mathcal{B}\equiv-N(\widetilde{u_{z}\phi}-\widetilde{{u}}_{z}\widetilde{\phi})$,
the small-scale TKE dissipation-rate is $\varepsilon_{K}\equiv
2\nu(\widetilde{\|\boldsymbol{s}\|^{2}}-\|\widetilde{\boldsymbol{s}}\|^{2})$,
the small-scale forcing is $\mathcal{F}_{K}\equiv\widetilde{\boldsymbol{F\cdot
u}}-\widetilde{\boldsymbol{F}}\boldsymbol{\cdot}\widetilde{\boldsymbol{u}}$,
and the small-scale TKE transport term involves

$\boldsymbol{T}_{K}\equiv-(1/\rho_{r})(\widetilde{\boldsymbol{u}p}-\widetilde{\boldsymbol{u}}\widetilde{p})-(1/2)(\widetilde{\boldsymbol{u}\|\boldsymbol{u}\|^{2}}-\widetilde{\boldsymbol{u}}\widetilde{\|\boldsymbol{u}\|^{2}})+(\widetilde{\boldsymbol{u}}\boldsymbol{\cdot\tau})+2\nu(\widetilde{\boldsymbol{u\cdot
s}}-\widetilde{\boldsymbol{u}}\boldsymbol{\cdot}\widetilde{\boldsymbol{s}}).$

The turbulent potential energy (TPE) in a given region of size $\ell$ is
$\widetilde{\phi\phi}$, and this may be split up into the amount contained in
the large scales denoted by $E_{P}\equiv\widetilde{\phi}\widetilde{\phi}/2$
and the amount contained in the small-scales denoted by
$e_{P}\equiv(\widetilde{\phi\phi}-\widetilde{\phi}\widetilde{\phi})/2$. The
equations governing $E_{P}$ and $e_{P}$ may be derived from (2) and are given
by

$\displaystyle D_{t}E_{P}$
$\displaystyle=-\boldsymbol{\nabla\cdot}(\widetilde{\phi}\boldsymbol{\Sigma})-\Pi_{P}+\frac{1}{2}\kappa\nabla^{2}\widetilde{\phi}\widetilde{\phi}-\kappa\|\boldsymbol{\nabla}\widetilde{\phi}\|^{2}+N\widetilde{\phi}\widetilde{u}_{z}+\widetilde{\phi}\widetilde{f},$
(7)
$\displaystyle\begin{split}D_{t}e_{P}&=\boldsymbol{\nabla\cdot}\boldsymbol{T}_{P}-\varepsilon_{P}-\mathcal{B}+\Pi_{P}+\mathcal{F}_{P}.\end{split}$
(8)

In these equations, the scale-to-scale TPE flux is defined as
$\Pi_{P}\equiv-\boldsymbol{\Sigma\cdot\nabla}\widetilde{\phi}$, such that
$\Pi_{P}>0$ corresponds to a transfer of TPE from the large to the small
scales. The small-scale TPE dissipation-rate is
$\varepsilon_{P}\equiv\kappa(\widetilde{\|\boldsymbol{\nabla}\phi\|^{2}}-\|\boldsymbol{\nabla}\widetilde{\phi}\|^{2})$,
the small-scale scalar forcing is $\mathcal{F}_{P}\equiv\widetilde{\phi
f}-\widetilde{\phi}\widetilde{f}$, and the small-scale scalar transport
involves

$\boldsymbol{T}_{P}\equiv(\kappa/2)\boldsymbol{\nabla}(\widetilde{\phi\phi}-\widetilde{\phi}\widetilde{\phi})+\widetilde{\phi}\boldsymbol{\Sigma}.$

In the following, attention will be given to the behavior of the terms that
contribute to $D_{t}e_{K}$ and $D_{t}e_{P}$. Note that while $e_{K}$ and
$e_{P}$ are referred to as the small-scale TKE and TPE, respectively, they
correspond to the TKE and TPE contained in all scales $<\ell$, so that when
$\ell$ exceeds the integral legnthscale of the flow, $e_{K}$ and $e_{P}$
actually contain the contributions from the largest dynamical scales in the
flow.

###  2.2 Length scales in SSST

As discussed in Portwood et al. (2019), for the velocity field there are four
important length scales for the SSST flow under consideration. The first is
the large-eddy lengthscale which may be characterized by
$L=\mathcal{E}_{K}^{3/2}/\langle\epsilon_{K}\rangle$, where
$\mathcal{E}_{K}\equiv\lim_{\ell\to\infty}\langle e_{K}\rangle$ is the total
(i.e. involving contributions from all scales) mean TKE in the flow, and
$\langle\epsilon_{K}\rangle\equiv\lim_{\ell\to\infty}\langle\varepsilon_{K}\rangle$
is the total mean TKE dissipation-rate. The second is the Kolmogorov
lengthscale $\eta=(\nu^{3}/\langle\epsilon_{K}\rangle)^{1/4}$ which
characterizes the scale at which viscous and inertial scales in the flow are
of the same order. The third is the Ozmidov lengthscale
$\ell_{O}=(\langle\epsilon_{K}\rangle/N^{3})^{1/2}$, which characterizes the
scale at which buoyancy and inertial forces are of the same order. The fourth
is the Corrsin lengthscale
$\ell_{C}=(\langle\epsilon_{K}\rangle/\gamma^{3})^{1/2}$, which characterizes
the scale at which mean-shear and inertial forces are of the same order.
Different nondimensional flow parameters may be related to these lengthcales
(assuming $\nu=\kappa$)

$\displaystyle Ri$
$\displaystyle\equiv\frac{N^{2}}{\gamma^{2}}=\Big{(}\ell_{C}/\ell_{O}\Big{)}^{4/3},$
(9) $\displaystyle Fr$
$\displaystyle\equiv\frac{\langle\epsilon_{K}\rangle}{N\mathcal{E}_{K}}=\Big{(}\ell_{O}/L\Big{)}^{2/3},$
(10) $\displaystyle Re_{b}$
$\displaystyle\equiv\frac{\langle\epsilon_{K}\rangle}{\nu
N^{2}}=\Big{(}\ell_{O}/\eta\Big{)}^{4/3},$ (11) $\displaystyle Re_{s}$
$\displaystyle\equiv\frac{\langle\epsilon_{K}\rangle}{\nu\gamma^{2}}=\Big{(}\ell_{C}/\eta\Big{)}^{4/3}.$
(12)

In order for the flow to become turbulent, the Richardson number must satisfy
$Ri<O(1)$ (and therefore $\ell_{C}<O(\ell_{O})$), and in order for
stratification to have an impact on the flow the Froude number must satisfy
$Fr<O(1)$. In order for there is to exist a range scales $\ell\ll\ell_{O}$
that are not affected by buoyancy, the buoyancy Reynolds number must satisfy
$Re_{b}\gg 1$. Correspondingly, in order for there is to exist a range scales
$\ell\ll\ell_{C}$ that are not affected by the mean-shear the shear Reynolds
number must satisfy $Re_{s}\gg 1$. These Reynolds numbers are, however,
related since $Re_{s}=Re_{b}Ri$.

For the scalar field, the large scale is
$L_{P}\equiv\mathcal{E}_{P}^{3/2}\langle\epsilon_{K}\rangle^{1/2}/\langle\epsilon_{P}\rangle^{3/2}$,
where $\mathcal{E}_{P}\equiv\lim_{\ell\to\infty}\langle e_{P}\rangle$ is the
total mean TPE in the flow, and
$\langle\epsilon_{P}\rangle\equiv\lim_{\ell\to\infty}\langle\varepsilon_{P}\rangle$
is the total mean TPE dissipation-rate. The small scale is the Batchelor scale
$\eta_{B}\equiv(\nu\kappa^{2}/\langle\epsilon_{K}\rangle)^{1/4}$. In this
paper we focus on flows with $\nu=\kappa$ such that $\eta_{B}=\eta$.

###  2.3 Mean-field behavior

We now turn to consider the contributions to the equations governing the mean
energy fields $\langle e_{K}\rangle$ and $\langle e_{P}\rangle$. For SSST, the
mean transport equations for these quantities reduce to

$\displaystyle 0$
$\displaystyle=\langle\mathcal{B}\rangle-\langle\varepsilon_{K}\rangle+\langle\Pi_{K}\rangle+\langle\mathcal{F}_{K}\rangle,$
(13) $\displaystyle 0$
$\displaystyle=-\langle\mathcal{B}\rangle-\langle\varepsilon_{P}\rangle+\langle\Pi_{P}\rangle.$
(14)

In the following discussion, we will first consider the case where the forcing
is confined to the large scales (such as in Lindborg (2006)), a flow we here
refer to as forced stably stratified turbulence (FSST), for which we introduce
the scale $\ell_{F}$ as the scale below which the forcing plays a sub-leading
role in in the flow. This behavior will then be compared to that of SSST to
understand the differences.

Focussing first on the TKE, in the limit $\ell/\eta\to\infty$,
$\langle\Pi_{K}\rangle\to 0$,
$\langle\varepsilon_{K}\rangle\to\langle\epsilon_{K}\rangle$,
$\langle\mathcal{F}_{K}\rangle\to\langle\mathcal{F}_{K}\rangle^{\infty}$,
$\langle\mathcal{B}\rangle\to\langle\mathcal{B}\rangle^{\infty}$ , so that we
have

$\displaystyle\langle\mathcal{F}_{K}\rangle^{\infty}$
$\displaystyle\sim-\langle\mathcal{B}\rangle^{\infty}+\langle\epsilon_{K}\rangle,$
(15)

reflecting a balance between the total injection of TKE by the forcing and the
total energy lost due to viscous dissipation and conversion to TPE. The
behavior of (13) and (14) as the scale $\ell$ decreases depends on the
dynamical scales of the system that were introduced in §2.2.

For FSST
$\langle\Pi_{K}\rangle\sim-\langle\mathcal{B}\rangle+\langle\epsilon_{K}\rangle$
when $\ell_{F}\gg\ell\gg\eta$, and for $\ell_{O}/\ell_{F}\to\infty$ (neutrally
buoyant), this would correspond to a TKE cascade
$\langle\Pi_{K}\rangle\sim\langle\epsilon_{K}\rangle$. However, for a stably
stratified flow where $\langle\mathcal{B}\rangle^{\infty}<0$, since buoyancy
effects reduce with decreasing scale,
$\langle\Pi_{K}\rangle\sim-\langle\mathcal{B}\rangle+\langle\epsilon_{K}\rangle$
implies that the TKE flux $\langle\Pi_{K}\rangle$ will actually reduce as
$\ell$ decreases, until it approaches a constant value in the regime
$\ell_{O}>\ell\gg\eta$, as discussed in Riley & Lindborg (2012); Kumar et al.
(2014). As a result of this, the TKE flux in stratified turbulence cannot be
in the form of a cascade (which would require a constant energy flux). This is
simply a reflection of the fact that TKE is being lost as it is passed down to
smaller scales due to conversion of TKE to TPE.

For SSST, the behavior of $\langle\Pi_{K}\rangle$ is quite different. In this
case, since $\langle\mathcal{F}_{K}\rangle$ operates down to the scale
$\ell_{C}$, then because $\ell_{C}<\ell_{O}$ the behavior
$\langle\Pi_{K}\rangle\sim-\langle\mathcal{B}\rangle+\langle\epsilon_{K}\rangle$
never emerges, and the regime that emerges instead for $Ri\ll 1$ is
$\langle\Pi_{K}\rangle\sim-\langle\mathcal{F}_{K}\rangle+\langle\epsilon_{K}\rangle$
for $\ell\gg\eta$. In this regime, the TKE flux is again not in the form of a
cascade since $\langle\mathcal{F}_{K}\rangle$ depends on $\ell$. Indeed, based
on the behavior of the co-spectrum (Katul et al., 2013) we might expect the
following behavior to emerge at high $Re$ in SSST

$\displaystyle\langle\Pi_{K}\rangle\sim-\gamma^{2}\langle\epsilon_{K}\rangle^{1/3}\ell^{4/3}+\langle\epsilon_{K}\rangle,\quad\text{for}\quad
L\gg\ell\gg\eta.$ (16)

Therefore, in this range $\langle\Pi_{K}\rangle$ increases as $\ell$
decreases, until it asymptotes to the constant flux cascade regime
$\langle\Pi_{K}\rangle\sim\langle\epsilon_{K}\rangle$ for
$\ell_{C}\gg\ell\gg\eta$. This is the opposite behavior to that discussed in
Kumar et al. (2014) for FSST where the forcing is confined to scales
$\gg\ell_{O}$ and where $\langle\Pi_{K}\rangle$ reduces as $\ell$ decreases
towards $\ell_{O}$.

In the above discussion, it has been assumed that
$\langle\varepsilon_{K}\rangle\sim\langle\epsilon_{K}\rangle$ for
$\ell\gg\eta$. While this is true for isotropic turbulence, it will not in
general apply in strongly stratified flows with $Fr\ll 1$. This is because
when $Fr\ll 1$, the flow structures are highly anisotropic, and so while the
horizontal lengthscale may be large enough for viscous effects to be
negligible, the vertical lengthscale of the structure may be small enough for
viscous effects to be important (Riley & Lindborg, 2012). In such a situation,
even if $\langle\Pi_{K}\rangle\sim\langle\varepsilon_{K}\rangle$ emerges in
some range of scales, this would not correspond to an inertial cascade since
$\langle\varepsilon_{K}\rangle$ would depend on $\ell$.

Concerning the TPE field, for flows with $\nu=\kappa$, then for $\ell\gg\eta$

$\displaystyle\langle\Pi_{P}\rangle\sim\langle\mathcal{B}\rangle+\langle\epsilon_{P}\rangle.$
(17)

Hence, in the regime $\ell>\ell_{O}$, as $\ell$ decreases,
$\langle\Pi_{P}\rangle$ increases, while for $\ell_{O}\gg\ell\gg\eta$, the
constant flux cascade regime
$\langle\Pi_{P}\rangle\sim\langle\epsilon_{P}\rangle$ emerges. Therefore,
there cannot be a TPE cascade in the strict sense at scales where density is
an active scalar, i.e. $\ell>\ell_{O}$, since the scale-dependency of
$\langle\mathcal{B}\rangle$ at these scales leads to a non-constant flux of
TPE at these scales. This observation seems to be in conflict with the well-
known Bolgiano-Obhukov (BO) scaling (Bolgiano Jr., 1959; Obukhov, 1959) that
has been proposed for stratified turbulent flows, in which
$\langle\Pi_{K}\rangle$ is supposed to decay with reducing $\ell$ due to
buoyancy as $\langle\Pi_{K}\rangle\propto\ell^{4/5}$ while
$\langle\Pi_{P}\rangle$ is constant
$\langle\Pi_{P}\rangle\sim\langle\epsilon_{P}\rangle$. This behavior, however,
seems problematic. For example, in FSST then for $\ell_{F}\gg\ell\gg\eta$

$\displaystyle\langle\Pi_{K}\rangle$
$\displaystyle\sim-\langle\mathcal{B}\rangle+\langle\epsilon_{K}\rangle,$ (18)
$\displaystyle\langle\Pi_{P}\rangle$
$\displaystyle\sim\langle\mathcal{B}\rangle+\langle\epsilon_{P}\rangle.$ (19)

It would seem that BO scaling could only emerge if
$\langle\epsilon_{P}\rangle\gg\langle\epsilon_{K}\rangle=O(|\langle\mathcal{B}\rangle|)$,
as noted in Alam et al. (2019), so that we could have
$\langle\Pi_{K}\rangle\sim-\langle\mathcal{B}\rangle+\langle\epsilon_{K}\rangle$
in which $\langle\Pi_{K}\rangle$ reduces as $\ell$ reduces while
$\langle\Pi_{P}\rangle\sim\langle\epsilon_{P}\rangle$. However, DNS of FSST
show $\langle\epsilon_{P}\rangle\lesssim\langle\epsilon_{K}\rangle$ over a
range of $Fr$ and $Re_{b}$ (Lindborg, 2006; de Bruyn Kops, 2015) which
includes the “moderately stratified” regime $Re\gg 1$ and $Fr\approx 1$ for
which it has been argued that BO should apply (Alam et al., 2019). In SSST, BO
scaling cannot emerge because there are by definition no scales at which the
dynamics are given by a balance between inertia and buoyancy forces, since the
mean-shear production is larger than the buoyancy at all scales at which it is
active, reflected in that SSST requires $\ell_{C}<\ell_{O}$.

###  2.4 Mechanisms governing the TKE and TPE fluxes

Since stable stratification acts to two-dimensionalize the flow, it had been
conjectured that such flows might feature an inverse energy cascade (Gage,
1979; Lilly, 1983). However, it has now been demonstrated that this does not
occur, but that there is a forward/downscale cascade of both kinetic and
potential energy (Riley & de Bruyn Kops, 2003; Waite & Bartello, 2006;
Lindborg, 2005, 2006; Brethouwer et al., 2007). The TKE cascade in
homogeneous, isotropic turbulence (HIT) is also downscale on average in three-
dimensions, however, the underlying mechanisms driving the energy transfers in
HIT and SSST might be quite different. In the context of HIT, it has long been
thought that the key mechanism driving the TKE cascade is that of vortex
stretching (VS) (Taylor, 1922, 1938; Tennekes & Lumley, 1972; Davidson, 2004;
Doan et al., 2018). However, recent studies have demonstrated quantitatively
that while VS plays an important role, the largest contribution to the TKE
cascade comes from the dynamical process of the self-amplification of the
strain-rate field (Carbone & Bragg, 2020; Johnson, 2020, 2021). An important
question is how this understanding applies to SSST, where effects such as
internal waves and mean-shear can play a role (whether directly or indirectly)
in how TKE is transferred between scales, as well as the question of the
mechanism driving the TPE transfer.

In Johnson (2020, 2021) a powerful, exact relationship was derived for
$\Pi_{K}$ that assumes only that the filtering kernel $\mathcal{G}_{\ell}$
used in constructing $\widetilde{\boldsymbol{u}}$ is Gaussian. The result is

$\displaystyle\Pi_{K}$
$\displaystyle=\Pi_{K}^{F,SSA}+\Pi_{K}^{F,VS}+\Pi_{K}^{SG,SSA}+\Pi_{K}^{SG,VS}+\Pi_{K}^{SG,C},$
(20) $\displaystyle\Pi_{K}^{F,SSA}$
$\displaystyle=-\ell^{2}\boldsymbol{\widetilde{s}}\boldsymbol{:}(\boldsymbol{\widetilde{s}}\boldsymbol{\cdot}\boldsymbol{\widetilde{s}}),$
(21) $\displaystyle\Pi_{K}^{F,VS}$
$\displaystyle=(1/4)\ell^{2}\boldsymbol{\widetilde{s}}\boldsymbol{:}(\boldsymbol{\widetilde{\omega}}\boldsymbol{\widetilde{\omega}}),$
(22)

where superscript ‘$F$’ denotes that the quantity only depends on the filtered
fields, while superscript ‘$SG$’ denotes that the quantity depends on the sub-
grid fields as well as the filtered fields. Explicit integral formalae for the
$SG$ terms can be found in Johnson (2020, 2021); we do not quote them here as
they will not be considered in detail in the present paper. Note that in
Johnson (2020, 2021), the ‘$F$’ contributions are referred to as the “local”
contributions while the ‘$SG$’ contributions are referred to as “non-local”
contributions. While the use of this terminology is appropriate in studies of
turbulence (L’vov & Falkovich, 1992; Aluie & Eyink, 2009), we prefer to avoid
it, because in turbulence, a local energy flux is usually typically taken to
imply an energy flux dominated by the interaction between scales of similar
size, which is not necessarily a property possessed by $\Pi_{K}^{F,SSA}$ and
$\Pi_{K}^{F,VS}$ which in principle may involve the interaction of any scales
in the flow of size $\geq\ell$.

The term $\Pi_{K}^{F,SSA}$ describes the contribution to the TKE flux arising
from the process of strain self-amplification (SSA) whereby the filtered
strain-rate field $\boldsymbol{\widetilde{s}}$ interacts with itself due to
nonlinearity to either amplify (if
$\boldsymbol{\widetilde{s}}\boldsymbol{:}(\boldsymbol{\widetilde{s}}\boldsymbol{\cdot}\boldsymbol{\widetilde{s}})<0$)
or else suppress (if
$\boldsymbol{\widetilde{s}}\boldsymbol{:}(\boldsymbol{\widetilde{s}}\boldsymbol{\cdot}\boldsymbol{\widetilde{s}})>0$)
the strain-rate magnitude $\|\boldsymbol{\widetilde{s}}\|$. The term
$\Pi_{K}^{F,VS}$ describes the contribution to the TKE flux arising from the
process of vortex stretching (VS) whereby the filtered strain-rate field
$\boldsymbol{\widetilde{s}}$ either amplifies (if
$\boldsymbol{\widetilde{s}}\boldsymbol{:}(\boldsymbol{\widetilde{\omega}}\boldsymbol{\widetilde{\omega}})>0$)
or else suppresses (if
$\boldsymbol{\widetilde{s}}\boldsymbol{:}(\boldsymbol{\widetilde{\omega}}\boldsymbol{\widetilde{\omega}})<0$)
the enstrophy $\|\boldsymbol{\widetilde{\omega}}\|^{2}$ through the process of
vortex stretching (or compression).

The average of the contribution from the filtered field
$\Pi_{K}^{F}=\Pi_{K}^{F,SSA}+\Pi_{K}^{F,VS}$ in (20) is similar to the
expression for the TKE flux in a two-point Karman-Hawarth equation derived
using filtering and an asymptotic expansion (Carbone & Bragg, 2020). Moreover,
due to the relation of Betchov (1956),
$\langle\Pi_{K}^{F,SSA}\rangle=3\langle\Pi_{K}^{F,VS}\rangle$ for
incompressible, homogeneous turbulence, so that the contribution of SSA to
$\langle\Pi_{K}^{F}\rangle$ is three times larger than that from VS. The sub-
grid contributions $\Pi_{K}^{SG,SSA}$ and $\Pi_{K}^{SG,VS}$ are similar to
$\Pi_{K}^{F,SSA}$ and $\Pi_{K}^{F,VS}$ except that in the $SG$ contributions,
the filtered strain-rate $\boldsymbol{\widetilde{s}}$ acts to amply the
strain-rate at sub-grid scales (for $\Pi_{K}^{SG,SSA}$) and vorticity at sub-
grid scales (for $\Pi_{K}^{SG,VS}$). DNS data in Johnson (2020, 2021) showed
that $\langle\Pi_{K}^{SG,SSA}\rangle\approx\langle\Pi_{K}^{SG,VS}\rangle$ in
the inertial range of isotropic turbulence, so that overall, the strain-self
amplification mechanism contributes more to the average energy cascade
$\langle\Pi_{K}\rangle$ than vortex stretching does, contrary to the
traditional explanation according to which vortex stretching is seen as the
key mechanism driving the energy cascade (Taylor, 1938; Tennekes & Lumley,
1972; Davidson, 2004; Doan et al., 2018).

The result in (20) also applies to SSST because it makes no assumption about
the flow dynamics. As such, in SSST, SSA and VS will still be the key
dynamical processes governing the TKE flux. Moreover, the relation of Betchov
(1956) still applies (since it only assumes incompressibility and homogeneity;
it does not assume isotropy). Therefore, at least with respect to the filtered
field contribution $\langle\Pi_{K}^{F}\rangle$, SSA still contributes three
times as much as VS to the total TKE flux. The relative contribution of
$\langle\Pi_{K}^{SG,SSA}\rangle$ and $\langle\Pi_{K}^{SG,VS}\rangle$ to
$\langle\Pi_{K}^{SG}\rangle$, however, may differ from that in HIT. Moreover,
the actual behavior of the SSA and VS processes in SSST may differ appreciably
from that in isotropic turbulence, since, for example, in SSST internal waves
can contribute to the behavior of $\boldsymbol{\widetilde{s}}$ and
$\boldsymbol{\widetilde{\omega}}$, and anisotropy in the flow can modify the
alignments between $\boldsymbol{\widetilde{s}}$ and
$\boldsymbol{\widetilde{\omega}}$ (see Sujovolsky & Mininni (2020)) which
affects the VS process. Some of these more involved questions will be the
subject of a forthcoming work. A key point to be explored here is the relative
contributions of $\Pi_{K}^{F}$ and $\Pi_{K}^{SG}$ to $\Pi_{K}$, and the
correlations between these terms which is important to understand for LES
modeling of SSST.

Following the same procedure outlined in Johnson (2020, 2021), a result
analogous to (20) can be derived for $\Pi_{P}$

$\displaystyle\Pi_{P}$
$\displaystyle=\Pi_{P}^{F}+\Pi_{P}^{SG,S}+\Pi_{P}^{SG,V},$ (23)
$\displaystyle\Pi_{P}^{F}$
$\displaystyle=-\ell^{2}\boldsymbol{\widetilde{s}}\boldsymbol{:}(\boldsymbol{\nabla}\widetilde{\phi}\boldsymbol{\nabla}\widetilde{\phi}),$
(24) $\displaystyle\Pi_{P}^{SG,S}$
$\displaystyle=-\boldsymbol{\nabla}\widetilde{\phi}\boldsymbol{\cdot}\int_{0}^{\ell^{2}}\tau_{\beta}\Big{(}\boldsymbol{\nabla}\widetilde{\phi}^{\sqrt{\alpha}},\boldsymbol{\widetilde{s}}^{\sqrt{\alpha}}\Big{)}\,d\alpha,$
(25) $\displaystyle\Pi_{P}^{SG,V}$
$\displaystyle=-\boldsymbol{\nabla}\widetilde{\phi}\boldsymbol{\cdot}\int_{0}^{\ell^{2}}\tau_{\beta}\Big{(}\boldsymbol{\nabla}\widetilde{\phi}^{\sqrt{\alpha}},\boldsymbol{\widetilde{r}}^{\sqrt{\alpha}}\Big{)}\,d\alpha,$
(26)

where $\widetilde{(\cdot)}^{\sqrt{\alpha}}$ denotes filtering at scale
$\sqrt{\alpha}$ (rather than scale $\ell$), and for arbitrary first order
$\boldsymbol{p}$ and second order $\boldsymbol{q}$ tensor fields,
$\tau_{\beta}$ is defined as
$\tau_{\beta}(\boldsymbol{p},\boldsymbol{q})\equiv\widetilde{\boldsymbol{p\cdot
q}}^{\beta}-\widetilde{\boldsymbol{p}}^{\beta}\boldsymbol{\cdot}\widetilde{\boldsymbol{q}}^{\beta}$
with $\beta\equiv\sqrt{\ell^{2}-\alpha}$, and
${\boldsymbol{r}}\equiv(\boldsymbol{\nabla}{\boldsymbol{u}}-[\boldsymbol{\nabla}{\boldsymbol{u}}]^{\top})/2$
is the rotation-rate tensor. The rotational motion in the flow makes no
explicit contribution to the filtered flux $\Pi_{P}^{F}$ because
$\boldsymbol{\widetilde{r}}\boldsymbol{:}(\boldsymbol{\nabla}\widetilde{\phi}\boldsymbol{\nabla}\widetilde{\phi})=0$,
although it does affect $\Pi_{P}^{F}$ implicitly since rotation in the fluid
affects the alignment between $\boldsymbol{\widetilde{s}}$ and
$\boldsymbol{\nabla}\widetilde{\phi}$.

The contribution from the filtered field $\Pi_{P}^{F}$ describes the flux of
TPE associated with the amplification (if
$\boldsymbol{\widetilde{s}}\boldsymbol{:}(\boldsymbol{\nabla}\widetilde{\phi}\boldsymbol{\nabla}\widetilde{\phi})<0$)
or suppression (if
$\boldsymbol{\widetilde{s}}\boldsymbol{:}(\boldsymbol{\nabla}\widetilde{\phi}\boldsymbol{\nabla}\widetilde{\phi})>0$)
of $\|\boldsymbol{\nabla}\widetilde{\phi}\|$ due to the filtered strain-rate
$\boldsymbol{\widetilde{s}}$. The contribution $\Pi_{P}^{F}$ is similar in
form to $\Pi_{K}^{F,VS}$, in that both depend upon the strain-rate field
$\boldsymbol{\widetilde{s}}$ acting to amplify or suppress another dynamical
field, rather than SSA in which the strain-rate field amplifies or suppresses
itself. However, the negative sign appearing in $\Pi_{P}^{F}$ but not in
$\Pi_{K}^{F,VS}$ leads to a significant difference in how the strain-rate
contributes to these energy fluxes. In particular, if we define
$\widetilde{\lambda}_{i}$ and $\widetilde{\boldsymbol{v}}_{i}$ as the $i^{th}$
ordered eigenvalues and eigenvectors of $\boldsymbol{\widetilde{s}}$, then
$\langle\Pi_{K}^{F,VS}\rangle=(1/4)\ell^{2}\langle\widetilde{\lambda}_{i}\|\boldsymbol{\widetilde{\omega}}\|^{2}\cos^{2}\theta_{\omega,i}\rangle$
and
$\langle\Pi_{P}^{F}\rangle=-\ell^{2}\langle\widetilde{\lambda}_{i}\|\boldsymbol{\nabla\widetilde{\phi}}\|^{2}\cos^{2}\theta_{\nabla\phi,i}\rangle$,
where $\theta_{\omega,i}$ denotes the angle between
$\boldsymbol{\widetilde{\omega}}$ and $\boldsymbol{v}_{i}$, and similarly for
$\theta_{\nabla\phi,i}$. In view of this, if $\langle\Pi_{K}^{F,VS}\rangle>0$
so that VS contributes to the average downscale TKE flux, then
$\langle\widetilde{\lambda}_{i}\|\boldsymbol{\widetilde{\omega}}\|^{2}\cos^{2}\theta_{\omega,i}\rangle$
must be dominated by contributions from the extensional eigendirections of
$\boldsymbol{\widetilde{s}}$. In isotropic turbulence and for $\ell/\eta\to
0$, $\boldsymbol{\widetilde{\omega}}$ preferentially aligns with
$\boldsymbol{v}_{2}$, but the contribution to
$\langle\widetilde{\lambda}_{i}\|\boldsymbol{\widetilde{\omega}}\|^{2}\cos^{2}\theta_{\omega,i}\rangle$
from $\widetilde{\lambda}_{1}$ is larger than that from
$\widetilde{\lambda}_{2}$ since $\widetilde{\lambda}_{1}$ tends to be much
larger than the positive values of $\widetilde{\lambda}_{2}$ (Tsinober, 2001;
Buaria et al., 2020). This is all the more true for $\ell/\eta$ in the
inertial range where the contribution from $\widetilde{\lambda}_{1}$ to
$\langle\widetilde{\lambda}_{i}\|\boldsymbol{\widetilde{\omega}}\|^{2}\cos^{2}\theta_{\omega,i}\rangle$
is much larger than that from $\widetilde{\lambda}_{2}$ (Carbone & Bragg,
2020). In an analogous way, if
$\langle\Pi_{P}^{F}\rangle=-\ell^{2}\langle\widetilde{\lambda}_{i}\|\boldsymbol{\nabla\widetilde{\phi}}\|^{2}\cos^{2}\theta_{\nabla\phi,i}\rangle>0$,
then
$\langle\widetilde{\lambda}_{i}\|\boldsymbol{\nabla\widetilde{\phi}}\|^{2}\cos^{2}\theta_{\nabla\phi,i}\rangle$
must be dominated by contributions from the compressive eigendirections of
$\boldsymbol{\widetilde{s}}$. Hence, while the VS contribution to the
downscale TKE flux is generated by the amplification of
$\boldsymbol{\widetilde{\omega}}$ due to the extensional straining motions in
the flow, the downscale TPE flux is generated by the amplification of
$\boldsymbol{\nabla}\widetilde{\phi}$ due to compressional straining motions
in the flow.

The sub-grid TPE flux associated with the strain-rate field $\Pi_{P}^{SG,S}$
and rotation-rate field $\Pi_{P}^{SG,V}$ do not have simple interpretations
(unlike $\Pi_{K}^{SG,SSA}$ and $\Pi_{K}^{SG,VS}$), but are related to the
amplification of the scalar gradients at multiple scales. It is interesting,
however, that the rotation-rate (and therefore vorticity) field makes an
explicit contribution to the TPE flux through the sub-grid term
$\Pi_{P}^{SG,V}$, even though it makes no explicit contribution to the
filtered flux $\Pi_{P}^{F}$.

The discussion above highlights that the physical mechanisms governing the TKE
and TPE fluxes in stratified turbulence are quite different, so that while
both are positive on average (e.g. Lindborg (2006)), it may not be reasonable
to model them in similar ways. We will explore this further by considering the
statistical correlation between $\Pi_{K}$ and $\Pi_{P}$ which will provide
insights into the extent to which they might be modeled using similar sub-grid
closures in the context of LES.

##  3 Direct Numerical Simulations

Table 1: Table of parameters in DNS. $Re_{b}$ | $Ri$ | $Fr$ | $L/\eta$ | $\ell_{O}/\eta$ | $\ell_{C}/\eta$
---|---|---|---|---|---
81 | 0.154 | 0.50 | 76.37 | 27 | 6.64

In the next section, data from DNS will be analyzed by computing various terms
in (6) and (8) to understand the processes and mechanisms controlling the
behavior of the TKE and TPE across scales in SSST. The DNS data used is from
the same data set presented in Portwood et al. (2019), which we here
summarize. In the DNS, the unfiltered versions of (1) and (2) are solved with
constant mean velocity gradient $\gamma$ and mean density gradient $\zeta$
using the Fourier pseudospectral scheme described in de Bruyn Kops (2015);
Almalkie & de Bruyn Kops (2012) but with the mean-shear term handled using an
integrating factor (Brucker et al., 2007; Chung & Matheou, 2012; Sekimoto et
al., 2016).

In order to generate a statistically stationary flow, a method similar to
Taylor et al. (2016) is adopted wherein the Richardson number of the flow is
adjusted via $g$ using a mass-spring damper control system with a target value
for the total TKE $(1/2)\|\boldsymbol{u}\|^{2}$. With this method, the
Richardson and Froude numbers are emergent in the flow, rather than imposed.
In the absence of stratification, the TKE in a homogeneous turbulent shear
flow grows with time (while the flows scales remain smaller than the domain
size). With stratification, mechanisms such as spontaneous shear instabilities
that form in the flow due to the layering provide for a downscale TKE flux. A
statistically stationary, homogeneous, sheared, stably stratified flow
therefore exists at a special equilibrium point where mechanisms generating a
downscale TKE flux prevent the growth of TKE that would occur in the flow in
the absence of stratification. Put another way, coupling of the fluctuating
momentum to the density fields provides for additional energy dissipation
mechanisms, which allows the sheared flow to attain a stead-state.

The flow is well-resolved, with isotropic grid spacing and maximum wavenumber
$k_{max}$ satisfying $k_{max}\eta\approx 2$. In order to allow for the
development of large anisotropic flow scales, a large domain of size
$\mathcal{L}_{x}/\mathcal{L}_{y}=2$, $\mathcal{L}_{x}/\mathcal{L}_{z}=4$,
$\mathcal{L}_{x}/L=40$ was used, where
$\mathcal{L}_{x},\mathcal{L}_{y},\mathcal{L}_{z}$ are the dimensions of the
domain in the $x,y,z,$ directions. In the $x$ direction 1792 grid points were
used. Table 1 summarizes the parameters in the DNS.

##  4 Results and Discussion

###  4.1 Mean-field behavior

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{tauij.pdf}
\put(100.0,-2.0){$\ell/\eta$} \put(98.0,90.0){\footnotesize$\langle
e_{K}\rangle/\langle E_{T}\rangle$}
\put(85.0,73.0){\footnotesize$\langle\tau_{xx}\rangle/[2\langle
E_{T}\rangle]$}
\put(85.0,58.0){\footnotesize$\langle\tau_{yy}\rangle/[2\langle
E_{T}\rangle]$}
\put(85.0,43.0){\footnotesize$\langle\tau_{zz}\rangle/[2\langle
E_{T}\rangle]$} \put(98.0,27.0){\footnotesize$\langle e_{P}\rangle/\langle
E_{T}\rangle$} \end{overpic} (a)  \begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{TKE_PE_dissipation.pdf}
\put(100.0,-2.0){$\ell/\eta$}
\put(98.0,42.0){\footnotesize$\langle\varepsilon_{K}\rangle/\langle\epsilon_{T}\rangle$}
\put(98.0,27.0){\footnotesize$\langle\varepsilon_{P}\rangle/\langle\epsilon_{T}\rangle$}
\end{overpic} (b)

Figure 1: (a) Plot of mean small-scale TKE $\langle e_{K}\rangle$, TPE
$\langle e_{P}\rangle$, and diagonal components of
$\langle\boldsymbol{\tau}\rangle/2$, normalized by total energy $\langle
E_{T}\rangle\equiv\lim_{\ell/\eta\to\infty}[\langle e_{K}\rangle+\langle
e_{P}\rangle]$, as a function of filter scale $\ell$. (b) Plot of mean small-
scale TKE $\langle\varepsilon_{K}\rangle$ and TPE
$\langle\varepsilon_{P}\rangle$ dissipation rates, normalized by the total
turbulent energy dissipation rate
$\langle\epsilon_{T}\rangle\equiv\lim_{\ell/\eta\to\infty}[\langle\varepsilon_{K}\rangle+\langle\varepsilon_{P}\rangle]$.

We begin by considering the behavior of $\langle e_{K}\rangle$, $\langle
e_{P}\rangle$ and the diagonal components of
$\langle\boldsymbol{\tau}\rangle/2$ (which correspond to the TKE associated
with different components of $\boldsymbol{u}$) as a function of filter length
$\ell$ in order to understand how energy is partitioned in the flow at
different scales. The results in figure 1(a) show that at all scales, $\langle
e_{K}\rangle>\langle e_{P}\rangle$, with $\langle e_{K}\rangle/\langle
e_{P}\rangle=O(10)$ at the large scales of the flow. The flow is therefore far
from a state of equipartition of large-scale energy among the TKE and TPE
fields, unlike the behavior that is thought to emerge for strongly stratified
flows with $Fr\ll 1$ (Billant & Chomaz, 2001). Below the Corrsin scale
$\ell_{C}=?$, the difference between $\langle e_{K}\rangle$ and $\langle
e_{P}\rangle$ reduces, but remains significant even for $\ell/\eta\ll 1$. When
considering the components $\langle\boldsymbol{\tau}\rangle/2$, it is seen
that although the total TKE is much larger than the TPE, the TKE associated
with particular components of the velocity field is comparable to the TPE. In
particular, $\langle\tau_{yy}\rangle/2\approx\langle e_{P}\rangle$ and
$\langle\tau_{zz}\rangle/2\approx(2/3)\langle e_{P}\rangle$ at the large
scales. As such, even though the energy contained in the TPE field is small
compared to that in the total TKE field, it is of the same order as that
contained in the spanwise ($y$) and vertical ($z$) directions of the flow, and
therefore plays an important energetic role in the system. In terms of
scaling, $\langle e_{K}\rangle$ and $\langle e_{P}\rangle$ both show the
expected behavior $\langle e_{K}\rangle\propto\ell^{2}$ and $\langle
e_{P}\rangle\propto\ell^{2}$ at $\ell/\eta\leq O(1)$ where viscous effects
lead to smooth velocity and density fields. Outside of this a well defined
scaling regime does not emerge, due to the limited Reynolds number of the flow
(as well as the fact that $\ell_{C}/\eta$ is too small for an inertially
dominated regime to emerge).

In figure 1(b) we consider the mean small-scale turbulent kinetic
$\langle\varepsilon_{K}\rangle$ and potential $\langle\varepsilon_{P}\rangle$
energy dissipation rates. For $\ell/\eta\to\infty$ these satisfy
$\langle\varepsilon_{K}\rangle\to\langle\epsilon_{K}\rangle$ and
$\langle\varepsilon_{P}\rangle\to\langle\epsilon_{P}\rangle$, and the results
show that these are in the ratio
$\langle\epsilon_{K}\rangle/\langle\epsilon_{P}\rangle\approx 5$. Both
$\langle\varepsilon_{K}\rangle$ and $\langle\varepsilon_{P}\rangle$ are only
approximately independent of $\ell$ down to $\ell/\eta=O(10)$, consistent with
the usual observation that the Kolmogorov scale underestimates the scale at
which viscous forces become important (Pope, 2000).

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{TKE_budget.pdf}
\put(80.0,-3.0){$\ell/\eta$}
\put(167.0,121.0){$-\langle\mathcal{B}\rangle/\langle\epsilon_{K}\rangle$}
\put(167.0,105.0){$\langle\Pi_{K}\rangle/\langle\epsilon_{K}\rangle$}
\put(167.0,89.0){$\langle\Pi_{K}^{F}\rangle/\langle\epsilon_{K}\rangle$}
\put(167.0,74.0){$-\langle\varepsilon_{K}\rangle/\langle\epsilon_{K}\rangle$}
\put(167.0,59.0){$\langle\mathcal{F}_{K}\rangle/\langle\epsilon_{K}\rangle$}
\end{overpic} Figure 2: Plot of terms in the average small-scale TKE budget
equation (13).

We now turn to consider the contributions to the equations governing $\langle
e_{K}\rangle$ and $\langle e_{P}\rangle$. The results in figure 2 show that
the flow exhibits the $\ell/\eta\to\infty$ asymptotic behavior
$\langle\mathcal{F}_{K}\rangle^{\infty}\sim-\langle\mathcal{B}\rangle^{\infty}+\langle\epsilon_{K}\rangle$
for $\ell/\eta\gtrsim O(100)$. As $\ell/\eta$ is decreased below $O(100)$,
$\langle\Pi_{K}\rangle$ begins to increase and becomes positive, indicating a
downscale flux of TKE, as has previously been observed for stratified
turbulent flows (Riley & de Bruyn Kops, 2003; Waite & Bartello, 2006;
Lindborg, 2005, 2006; Brethouwer et al., 2007). However, it does not become
significant until $\ell/\eta\lesssim O(10)$. Moreover, when it reaches its
peak value, it does not exhibit the behavior
$\langle\Pi_{K}\rangle/\langle\epsilon_{K}\rangle\sim 1$ which would be
expected for non-stratified isotropic turbulence. This is because
$\langle\Pi_{K}\rangle$ cannot grow significantly until $\ell$ is small enough
such that
$\langle\mathcal{F}_{K}\rangle\ll\langle\mathcal{F}_{K}\rangle^{\infty}$,
which only occurs for $\ell\leq O(\ell_{C})$. Once this regime is obtained,
then since $\ell_{C}<\ell_{O}$ the TKE balance becomes
$\langle\Pi_{K}\rangle\sim\langle\varepsilon_{K}\rangle$, a behavior that can
be observed in figure 2 to be approached at small scales. However, in the
present flow, this behavior does not give rise to an inertial TKE cascade
since at the scales where
$\langle\Pi_{K}\rangle\sim\langle\varepsilon_{K}\rangle$,
$\langle\varepsilon_{K}\rangle$ is a decreasing function of $\ell$. In order
to observe an inertial TKE cascade with
$\langle\Pi_{K}\rangle\sim\langle\epsilon_{K}\rangle$ would require
considering a flow possessing a range of scales
$\ell_{O}>\ell_{C}\gg\ell\gg\eta$. The results in figure 2 also show that the
contribution to the TKE flux coming from the filtered field dynamics, i.e.
$\langle\Pi_{K}^{F}\rangle$, makes a significant contribution to
$\langle\Pi_{K}\rangle$ at scales where $\langle\Pi_{K}\rangle$ plays a
significant role in the small-scale TKE budget equation, but also that the
contribution involving the sub-grid fields $\langle\Pi_{K}^{SG}\rangle$ is
also significant, just as for isotropic turbulence (Johnson, 2020, 2021). At
sufficiently small $\ell$,
$\langle\Pi_{K}\rangle\approx\langle\Pi_{K}^{F}\rangle$, consistent with the
exact limiting behavior $\lim_{\ell/\eta\to 0}\Pi_{K}^{F}\to\Pi_{K}$.

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{PE_budget.pdf}
\put(80.0,-3.0){$\ell/\eta$}
\put(167.0,121.0){$-\langle\mathcal{B}\rangle/\langle\epsilon_{P}\rangle$}
\put(167.0,105.0){$\langle\Pi_{P}\rangle/\langle\epsilon_{P}\rangle$}
\put(167.0,89.0){$\langle\Pi_{P}^{F}\rangle/\langle\epsilon_{P}\rangle$}
\put(167.0,74.0){$-\langle\varepsilon_{P}\rangle/\langle\epsilon_{P}\rangle$}
\end{overpic} Figure 3: Plot of terms in the average small-scale TPE budget
equation (14).

Concerning the TPE behavior, the results in figure 3 exhibit the
$\ell/\eta\to\infty$ asymptotic behavior
$\langle\mathcal{B}\rangle^{\infty}\sim-\langle\epsilon_{P}\rangle$ for
$\ell/\eta\gtrsim O(100)$. The behavior of $\langle\Pi_{P}\rangle$ is very
similar to that of $\langle\Pi_{K}\rangle$ in that $\langle\Pi_{P}\rangle$
does not become significant until $\ell/\eta\lesssim O(10)$, and does not
exhibit the behavior $\langle\Pi_{P}\rangle/\langle\epsilon_{P}\rangle\sim 1$
that would be expected for a passive scalar advected in isotropic turbulence.
Whereas the growth of $\langle\Pi_{K}\rangle$ is constrained in SSST to scales
$\ell\leq O(\ell_{C})$, as discussed previously, the growth of
$\langle\Pi_{P}\rangle$ is constrained in SSST to scales $\ell\leq
O(\ell_{O})$. In our DNS these two ranges are similar, and so the range of
scales over which $\langle\Pi_{K}\rangle$ and $\langle\Pi_{P}\rangle$ are
active is similar. However, in an SSST flow with $Ri\lll 1$ and $Re\ggg 1$
this could lead to an interesting scenario where
$\langle\Pi_{P}\rangle/\langle\epsilon_{P}\rangle\sim 1$ while
$\langle\Pi_{K}\rangle/\langle\epsilon_{K}\rangle\approx 0$ for
$\ell_{O}\gg\ell\gg\ell_{C}$. The results in figure 3 also show that like the
TKE results, the contribution to the TPE flux coming from the filtered field
dynamics, i.e. $\langle\Pi_{P}^{F}\rangle$, makes a significant contribution
to $\langle\Pi_{P}\rangle$ at scales where $\langle\Pi_{P}\rangle$ plays a
significant role in the small-scale TPE budget equation, but that the
contribution involving the sub-grid fields $\langle\Pi_{P}^{SG}\rangle$ is
also significant. At sufficiently small $\ell$,
$\langle\Pi_{P}\rangle\approx\langle\Pi_{P}^{F}\rangle$, consistent with the
exact limiting behavior $\lim_{\ell/\eta\to 0}\Pi_{P}^{F}\to\Pi_{P}$.

###  4.2 Fluctuations about the mean-field

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
213.39566pt,scale={0.32},clip,tics=20]{PDF_eK.pdf}
\put(83.0,0.0){$e_{K}/\langle e_{K}\rangle$}
\put(-3.0,60.0){\rotatebox{90.0}{PDF}} \put(108.0,120.0){$\ell/\eta=0.25$}
\put(108.0,109.0){$\ell/\eta=6$} \put(108.0,97.0){$\ell/\eta=16$}
\put(108.0,85.0){$\ell/\eta=60$} \end{overpic} (a)  \begin{overpic}[trim=0.0pt
170.71652pt 0.0pt 213.39566pt,scale={0.32},clip,tics=20]{PDF_eP.pdf}
\put(83.0,0.0){$e_{P}/\langle e_{P}\rangle$}
\put(-3.0,60.0){\rotatebox{90.0}{PDF}} \end{overpic} (b)

\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
213.39566pt,scale={0.32},clip,tics=20]{PDF_PiK_PiP.pdf} \put(99.0,0.0){$\xi$}
\put(-3.0,60.0){\rotatebox{90.0}{PDF}} \end{overpic} (c)
\begin{overpic}[trim=0.0pt 170.71652pt 0.0pt
213.39566pt,scale={0.32},clip,tics=20]{PDF_B.pdf}
\put(89.0,0.0){$\mathcal{B}/\langle\mathcal{B}\rangle$}
\put(-3.0,60.0){\rotatebox{90.0}{PDF}} \end{overpic} (d)

Figure 4: Plots of PDF of (a) $e_{K}/\langle e_{K}\rangle$, (b) $e_{P}/\langle
e_{P}\rangle$ and (d) $\mathcal{B}/\langle\mathcal{B}\rangle$ for different
filter lengths $\ell$. In plot (c), the solid lines correspond to the PDFs of
$\Pi_{K}/\langle\Pi_{K}\rangle$ and the dashed lines correspond to the PDFs of
$\Pi_{P}/\langle\Pi_{P}\rangle$. Different colors/symbols correspond to
different $\ell/\eta$ as indicated by the legend in plot (a).

In order to understand the energetics of the flow beyond its mean-field
behavior, we will now consider the Probability Density Functions (PDFs) of
various quantities. We begin in figure 4(a),(b) with the PDFs of the
normalized small-scale energies $e_{K}/\langle e_{K}\rangle$ and
$e_{P}/\langle e_{P}\rangle$ for different filtering lengths $\ell/\eta$. At
the small scales, the PDFs are highly non-Gaussian, such that the TKE and TPE
fields exhibit frequent large fluctuations about their mean-field behaviour.
As the filter length increases, the PDFs approach a Gaussian shape, associated
with the dominance of the linear mean-shear and buoyancy forces at the large
scales of the flow. The results also show that the TPE field exhibits stronger
fluctuations from the mean-field behavior than the TKE field. This is
consistent with the fact that scalar fields are known to be more intermittent
in turbulent flows because they lack the pressure gradient in their dynamics
that regulates large fluctuations in the velocity gradient field.

In figure 4(c) we show the PDFs of $\Pi_{K}/\langle\Pi_{K}\rangle$ and
$\Pi_{P}/\langle\Pi_{P}\rangle$ for different filtering lengths $\ell/\eta$.
At larger scales, the PDFs have a large variance, showing that at these
scales, the fluxes of TKE and TPE can significantly exceed their mean values.
The PDFs are also close to being symmetric, so that the probability of upscale
or downscale fluxes TKE are similar, and the same for the TPE. As one moves to
smaller scales where $\langle\Pi_{K}\rangle$ and $\langle\Pi_{P}\rangle$ both
increase, the probability of large fluctuations about the mean-field behavior
reduce. However, there is still a significant probability of observing regions
of the flow where the local values of $\Pi_{K}$ and $\Pi_{P}$ significantly
exceed their mean values in magnitude. As $\ell$ is reduced, the mode of the
PDF remains close to zero, while the mean and skewness both become positive,
with the probability of downscale fluxes significantly exceeding that of
upscale fluxes. The PDFs for TKE and TPE are similar in shape, and at scales
where $\langle\Pi_{K}\rangle$ and $\langle\Pi_{P}\rangle$ play an important
role in the TKE and TPE scale-wise budgets, the fluctuations of $\Pi_{K}$
about $\langle\Pi_{K}\rangle$ are larger than those of $\Pi_{P}$ about
$\langle\Pi_{P}\rangle$.

\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_PE_flux_1.pdf}
\put(80.0,8.0){$\Pi_{P}/\langle\Pi_{P}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{K}/\langle\Pi_{K}\rangle$}}
\end{overpic} (a)  \begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_PE_flux_2.pdf}
\put(80.0,8.0){$\Pi_{P}/\langle\Pi_{P}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{K}/\langle\Pi_{K}\rangle$}}
\end{overpic} (b)

\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_PE_flux_3.pdf}
\put(80.0,8.0){$\Pi_{P}/\langle\Pi_{P}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{K}/\langle\Pi_{K}\rangle$}}
\end{overpic} (c)  \begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_PE_flux_4.pdf}
\put(80.0,8.0){$\Pi_{P}/\langle\Pi_{P}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{K}/\langle\Pi_{K}\rangle$}}
\end{overpic} (d)

Figure 5: Contour plot of the logarithm of the joint PDF of
$\Pi_{P}/\langle\Pi_{P}\rangle$ and $\Pi_{K}/\langle\Pi_{K}\rangle$ for (a)
$\ell/\eta=0.25$, (b) $\ell/\eta=6$, (c) $\ell/\eta=16$, (d) $\ell/\eta=60$.
Colors correspond to the logarithm of the PDF.

In figure 4(d) we show the PDFs of $\mathcal{B}/\langle\mathcal{B}\rangle$ for
different filtering lengths $\ell/\eta$. At larger scales where
$-\langle\mathcal{B}\rangle$ is the dominant source term in the mean TPE
budget, the fluctuations about $\langle\mathcal{B}\rangle$ are not that
strong, and remarkably, the probability to observe
$\mathcal{B}/\langle\mathcal{B}\rangle<0$, which would correspond to
conversion of TPE to TKE, is zero. For smaller scales, the probability to
observe $\mathcal{B}/\langle\mathcal{B}\rangle<0$ becomes finite and increases
with decreasing $\ell$. The significant fluctuations of $\mathcal{B}$ about
its mean value indicate that the effects of buoyancy could be felt at scales
considerably smaller than the mean-field description suggests.

In figure 5 we show the joint PDFs of $\Pi_{P}/\langle\Pi_{P}\rangle$ and
$\Pi_{K}/\langle\Pi_{K}\rangle$ for different filtering lengths $\ell/\eta$.
The results reveal that there is a weak positive correlation between $\Pi_{P}$
and $\Pi_{K}$ at larger scales, which becomes stronger as $\ell/\eta$ is
increased. Nevertheless, even in the range of scales where the TKE and TPE
mean fluxes are strongest, the positive correlation between $\Pi_{P}$ and
$\Pi_{K}$ is not very strong. Indeed, the probability of observing a downscale
TKE flux with an upscale TPE flux and vice-versa, is quite significant. This
relatively weak correlation can be understood in view of the discussion in
§2.4 regarding the mechanisms governing the TKE and TPE fluxes. It is common
to model TKE and TPE fluxes using similar closures but with a coefficient that
accounts for their differences only through the Prandtl number effects.
However, the weak correlation between $\Pi_{K}$ and $\Pi_{P}$ revealed in
figure 5 shows that such an approach is not suitable, and that the closures
need to reflect to some extent the different physical processes dominating
these fluxes, which leads to the relatively weak correlation between $\Pi_{K}$
and $\Pi_{P}$.

\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_FSG_1.pdf}
\put(80.0,8.0){$\Pi_{K}^{F}/\langle\Pi_{K}^{F}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{K}/\langle\Pi_{K}\rangle$}}
\end{overpic} (a)  \begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_FSG_2.pdf}
\put(80.0,8.0){$\Pi_{K}^{F}/\langle\Pi_{K}^{F}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{K}/\langle\Pi_{K}\rangle$}}
\end{overpic} (b)

\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_FSG_3.pdf}
\put(80.0,8.0){$\Pi_{K}^{F}/\langle\Pi_{K}^{F}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{K}/\langle\Pi_{K}\rangle$}}
\end{overpic} (c)  \begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_FSG_4.pdf}
\put(80.0,8.0){$\Pi_{K}^{F}/\langle\Pi_{K}^{F}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{K}/\langle\Pi_{K}\rangle$}}
\end{overpic} (d)

Figure 6: Contour plot of the logarithm of the joint PDF of
$\Pi_{K}/\langle\Pi_{K}\rangle$ and $\Pi_{K}^{F}/\langle\Pi_{K}^{F}\rangle$
for (a) $\ell/\eta=0.25$, (b) $\ell/\eta=6$, (c) $\ell/\eta=16$, (d)
$\ell/\eta=60$. Colors correspond to the logarithm of the PDF.

In the context of LES, the contributions to $\Pi_{K}$ and $\Pi_{P}$ from the
filtered field dynamics, namely
$\Pi_{K}^{F}\equiv\Pi_{K}^{F,SSA}+\Pi_{K}^{F,VS}$ and $\Pi_{P}^{F}$ (see
§2.4), are closed since they depend only on the filtered field quantities, not
the sub-grid fields. It is therefore of particular interest to understand how
these terms contribute to the total TKE and TPE fluxes in the flow. To explore
this, in figure 6 we consider the joint PDF of $\Pi_{K}$ and $\Pi_{K}^{F}$ for
different filtering lengths $\ell/\eta$. The filtered fluxes satisfy
$\lim_{\ell/\eta\to 0}\Pi_{K}^{F}\to\Pi_{K}$, and the results for
$\ell/\eta=0.25$ are approaching this regime, showing a very small probability
of events deviating from the state $\Pi_{K}=\Pi_{K}^{F}$. As $\ell/\eta$ is
increased, there remains a significant positive correlation between the
variables $\Pi_{K}$ and $\Pi_{K}^{F}$, indicating that $\Pi_{K}^{F}$ makes an
important contribution to the total flux $\Pi_{K}$, consistent with the mean-
field results in figure 2.

\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_PE_FSG_1.pdf}
\put(80.0,8.0){$\Pi_{P}^{F}/\langle\Pi_{P}^{F}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{P}/\langle\Pi_{P}\rangle$}}
\end{overpic} (a)  \begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_PE_FSG_2.pdf}
\put(80.0,8.0){$\Pi_{P}^{F}/\langle\Pi_{P}^{F}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{P}/\langle\Pi_{P}\rangle$}}
\end{overpic} (b)

\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_PE_FSG_3.pdf}
\put(80.0,8.0){$\Pi_{P}^{F}/\langle\Pi_{P}^{F}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{P}/\langle\Pi_{P}\rangle$}}
\end{overpic} (c)  \begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_PE_FSG_4.pdf}
\put(80.0,8.0){$\Pi_{P}^{F}/\langle\Pi_{P}^{F}\rangle$}
\put(5.0,65.0){\rotatebox{90.0}{$\Pi_{P}/\langle\Pi_{P}\rangle$}}
\end{overpic} (d)

Figure 7: Contour plot of the logarithm of the joint PDF of
$\Pi_{P}/\langle\Pi_{P}\rangle$ and $\Pi_{P}^{F}/\langle\Pi_{P}^{F}\rangle$
for (a) $\ell/\eta=0.25$, (b) $\ell/\eta=6$, (c) $\ell/\eta=16$, (d)
$\ell/\eta=60$. Colors correspond to the logarithm of the PDF.

However, there is a significant probability of events where either
$\Pi_{K}^{F}$ and $\Pi_{K}$ have significantly different magnitudes, or else
even have opposite signs. In terms of the physical mechanisms discussed in
§2.4, this means that at a given scale $\ell$, if the strain and vorticity
fields are being amplified at that scale by nonlinearity such that
$\Pi_{K}^{F}>0$, this may nevertheless not lead to a downscale flux of TKE
associated with $\Pi_{K}>0$ since at scales smaller than $\ell$ the strain and
vorticity fields may be experiencing suppression (rather than amplification)
due to nonlinearity yielding $\Pi_{K}^{SG}<0$, and if this is strong enough
then $\Pi_{K}=\Pi_{K}^{F}+\Pi_{K}^{SG}<0$. The spread of the PDF about the
line $\Pi_{K}=\Pi_{K}^{F}$ also implies that the relative contribution of
$\Pi_{K}^{F}$ to $\Pi_{K}$ is similar both during events where
$\Pi_{K}\sim\langle\Pi_{K}\rangle$ and large fluctuations where
$|\Pi_{K}|\ll\langle\Pi_{K}\rangle$ or $|\Pi_{K}|\gg\langle\Pi_{K}\rangle$.
For $\ell/\eta=60$, the results show that $\Pi_{K}$ and $\Pi_{K}^{F}$ are
almost uncorrelated, and at this scale,
$\langle\Pi_{K}\rangle/\langle\Pi_{K}^{F}\rangle\approx 43$, so that the
filtered field makes a small, uncorrelated contribution to $\Pi_{K}$, which is
itself very small at $\ell/\eta=60$. The results for the joint PDF of
$\Pi_{P}$ and $\Pi_{P}^{F}$ are shown in figure 7. These are very similar to
the TKE flux results, with $\Pi_{P}^{F}\sim\Pi_{P}$ observed at the smallest
scales, and a significant correlation between the two at intermediate scales.
The correlation between $\Pi_{P}^{F}$ and $\Pi_{P}$ is, however, slightly
stronger than that between $\Pi_{K}$ and $\Pi_{K}^{F}$.

\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_B_1.pdf}
\put(5.0,70.0){\rotatebox{90.0}{$\mathcal{B}/\langle\mathcal{B}\rangle$}}
\put(80.0,8.0){$\Pi_{K}/\langle\Pi_{K}\rangle$} \end{overpic} (a)
\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_B_2.pdf}
\put(5.0,70.0){\rotatebox{90.0}{$\mathcal{B}/\langle\mathcal{B}\rangle$}}
\put(80.0,8.0){$\Pi_{K}/\langle\Pi_{K}\rangle$} \end{overpic} (b)

\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_B_3.pdf}
\put(5.0,70.0){\rotatebox{90.0}{$\mathcal{B}/\langle\mathcal{B}\rangle$}}
\put(80.0,8.0){$\Pi_{K}/\langle\Pi_{K}\rangle$} \end{overpic} (c)
\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_TKE_B_4.pdf}
\put(5.0,70.0){\rotatebox{90.0}{$\mathcal{B}/\langle\mathcal{B}\rangle$}}
\put(80.0,8.0){$\Pi_{K}/\langle\Pi_{K}\rangle$} \end{overpic} (d)

Figure 8: Contour plot of the logarithm of the joint PDF of
$\mathcal{B}/\langle\mathcal{B}\rangle$ and $\Pi_{K}/\langle\Pi_{K}\rangle$
for (a) $\ell/\eta=0.25$, (b) $\ell/\eta=6$, (c) $\ell/\eta=16$, (d)
$\ell/\eta=60$. Colors correspond to the logarithm of the PDF.

In figure 8 we consider the joint PDF of $\mathcal{B}$ and $\Pi_{K}$ for
different filter lengths $\ell$. The results show that as $\ell$ is decreased,
the PDF reorients from being extended along the $\Pi_{K}$ axis, to being
extended along the $\mathcal{B}$ axis. This corresponds to the transition from
larger scales where buoyancy plays a more dominant role in the TKE energetics,
to smaller scales where the nonlinear energy flux plays a dominant role. As
also observed in figure 4(d), at larger scales the probability to observe
$\mathcal{B}/\langle\mathcal{B}\rangle<0$ is very low, suggesting that
convective motion is very rare at these scales. As $\ell$ decreases, however,
the probability to observe $\mathcal{B}/\langle\mathcal{B}\rangle<0$ increases
significantly. At $\ell/\eta=16$,
$-\langle\mathcal{B}\rangle\approx\langle\Pi_{K}\rangle$, and the results in
figure 8 show that at this scale, convective motion
$\mathcal{B}/\langle\mathcal{B}\rangle<0$ can occur, but the PDF is strongly
skewed towards stably stratified regions that have
$\mathcal{B}/\langle\mathcal{B}\rangle>0$. Moreover, for this scale the PDF is
stretched significantly along the $\Pi_{K}$ axis, showing that fluctuations of
$\Pi_{K}$ are considerably stronger than those of $\mathcal{B}$ about their
respective mean values which are approximately equal in magnitude at this
scale. The results for the joint PDF of $\mathcal{B}$ and $\Pi_{P}$ are
similar to those of the joint PDF of $\mathcal{B}$ and $\Pi_{K}$, and
therefore for brevity we do not show them.

\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_QR_1.pdf}
\put(80.0,8.0){$R/\langle\|\widetilde{\boldsymbol{s}}\|^{2}\rangle^{3/2}$}
\put(5.0,65.0){\rotatebox{90.0}{$Q/\langle\|\widetilde{\boldsymbol{s}}\|^{2}\rangle$}}
\end{overpic} (a)  \begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_QR_2.pdf}
\put(80.0,8.0){$R/\langle\|\widetilde{\boldsymbol{s}}\|^{2}\rangle^{3/2}$}
\put(5.0,65.0){\rotatebox{90.0}{$Q/\langle\|\widetilde{\boldsymbol{s}}\|^{2}\rangle$}}
\end{overpic} (b)

\begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_QR_3.pdf}
\put(80.0,8.0){$R/\langle\|\widetilde{\boldsymbol{s}}\|^{2}\rangle^{3/2}$}
\put(5.0,65.0){\rotatebox{90.0}{$Q/\langle\|\widetilde{\boldsymbol{s}}\|^{2}\rangle$}}
\end{overpic} (c)  \begin{overpic}[trim=0.0pt 142.26378pt 0.0pt
199.16928pt,scale={0.32},clip,tics=20]{JPDF_QR_4.pdf}
\put(80.0,8.0){$R/\langle\|\widetilde{\boldsymbol{s}}\|^{2}\rangle^{3/2}$}
\put(5.0,65.0){\rotatebox{90.0}{$Q/\langle\|\widetilde{\boldsymbol{s}}\|^{2}\rangle$}}
\end{overpic} (d)

Figure 9: Contour plot of the logarithm of the joint PDF of $Q$ and $R$ for
(a) $\ell/\eta=0.25$, (b) $\ell/\eta=6$, (c) $\ell/\eta=16$, (d)
$\ell/\eta=60$. Colors correspond to the logarithm of the PDF.

Finally, given the importance of the velocity gradient dynamics for the TKE
flux, in figure 9 we plot the results for the joint PDF of the filtered
velocity gradient invariants
$Q\equiv{\boldsymbol{\nabla\widetilde{u}}\boldsymbol{:}\boldsymbol{\nabla\widetilde{u}}}$
and
$R\equiv{(\boldsymbol{\nabla\widetilde{u}}\boldsymbol{\cdot}\boldsymbol{\nabla\widetilde{u}})\boldsymbol{:}\boldsymbol{\nabla\widetilde{u}}}$
at different filter scales $\ell/\eta$. In Danish & Meneveau (2018) this joint
PDF was considered for isotropic turbulence, and they found that the classic
sheared-drop shape of the isoprobability lines is preserved as $\ell/\eta$ is
varied, while other statistical characterizations of the velocity gradients
also reveal behavior that is qualitatively similar as the filter scale is
increased (Tom et al., 2020). As a result of this, the change in shape of the
PDF contours shown in figure 9 as $\ell/\eta$ is increased is strictly due to
the competing influence of buoyancy and mean-shear with nonlinearity in the
flow. For $\ell/\eta=0.25$, the shape of the PDF is identical to that observed
in isotropic turbulence (Meneveau, 2011), with the contours extending down the
right Viellefosse tail, and a strong preference for quadrants $Q>0,R<0$ and
$Q<0,R>0$. As $\ell/\eta$ is increased, however, the PDF becomes more
symmetric in $R$. This is due to the fact that as $\ell/\eta$ is increased,
the nonlinear term that generates the asymmetry in $R$ (which is associated
with the preference for vortex stretching over compression, and strain-rate
self-amplification over self-suppression (Tsinober, 2001)) becomes subleading
compared to buoyancy and mean-shear, and these terms, being linear, generate
time-reversible dynamics with no preference for nonlinear amplification or
suppression of the strain-rate and vorticity fields. Figure 9 (d) still shows,
however, that for the outer isoprobability lines, an asymmetry still exists,
with a preference for the regions $Q>0,R<0$ and $Q<0,R>0$. This corresponds to
events in the flow where locally the nonlinear term dominates in the flow over
buoyancy, allowing for strong strain-rate and/or vorticity amplification in
that localized region.

##  5 Conclusions

In this paper, we have analyzed the scale-dependent TKE and TPE in sheared,
stably stratified turbulence (SSST) using a filtering approach, where the flow
has constant mean velocity gradient and mean density gradient. Equations for
the scale-dependent TKE and TPE are to explore the competing effects in the
flow as well as the physical mechanisms governing the TKE and TPE fluxes
between scales. Various quantities in these equations were then evaluated
using data from DNS of SSST, with attention given to both the mean field
behavior of the flow, as well as fluctuations about this mean field state.

In terms of the mean field properties, while the mean TKE $\langle
e_{K}\rangle$ is larger than the mean TPE $\langle e_{P}\rangle$ by an order
of magnitude at the large scales, the difference between them reduces with
reducing scale $\ell$. The TPE nevertheless plays an import energetic role in
the system, and is larger than the vertical component of TKE at all scales in
the flow. The mean small-scale TKE dissipation rate
$\langle\varepsilon_{K}\rangle$ is also significantly larger than the
potential dissipation rate $\langle\varepsilon_{P}\rangle$. The mean TKE and
TPE fluxes between scales, $\langle\Pi_{K}\rangle$ and
$\langle\Pi_{P}\rangle$, respectively, do not reveal a cascade regime in the
flow due to the impact of the mean-shear down to relatively small-scales in
the flow. The contributions to these fluxes from the filtered fields are also
shown the be significant at scales where the energy fluxes play significant
roles in the energy budget equations, however, the contribution from the sub-
grid fields are also import, just as has been observed for isotropic
turbulence.

To understand the flow energetics beyond the mean fields, PDFs of various
quantities (normalized by their mean values) have been studied and discussed.
The PDFs of small-scale TKE and TPE are highly non-Gaussian at the smallest
scales, indicating large fluctuations of the TKE and TPE about their mean
field values. The TPE shows stronger fluctuations from its mean field value
than the TKE, which is consistent with the known result that scalar fields are
more intermittent the velocity fields in turbulence, since scalar fields lack
a pressure gradient term in their equation which can act to suppress large
fluctuations. As the filter scale is increased, these PDFs approach a Gaussian
distribution at scales where the linear forces in the flow (mean-shear and
buoyancy) dominate the dynamics. At larger scales in the flow, buoyancy seems
always convert TKE to TPE, with the PDF of the buoyancy term $\mathcal{B}$
showing almost zero probability of locally convective events where
$\mathcal{B}>0$. The probability of locally convective regions increases as
the filter scale is decreased, however, the probability remains small at
scales downward to the Ozimidov scale. The TKE and TPE fluxes between scales
are on average positive (downscale flux), with their instantaneous values
positively correlated, and increasingly so with decreasing filter scale.
However, the correlation is quite weak and this is because the physical
mechanisms governing the TKE and TPE fluxes are quite different, as discussed
in our theoretical analysis in §2.4. Indeed, the joint PDFs of the TKE and TPE
fluxes reveals a significant probability of events where the TKE and TPE
fluxes may have very different magnitudes, and even opposite signs. Finally,
the PDF of the principal invariants of the filtered velocity gradients (the
so-called $Q,R$ invariants (Tsinober, 2001)) reveals the classical sheared-
drop shape for the smallest filter scales. However, the PDF becomes
increasingly symmetric about $R=0$ as the filter scale increases. This
contrasts with the behavior observed in isotropic turbulence (Danish &
Meneveau, 2018) and is due to the increasingly dominant effects of the linear
mean shear and buoyancy effects over the nonlinear inertial effects at these
scales. This symmetry that emerges means that vortex stretching and
compression become equally probably, as well as strain self-amplification and
self-suppression.

##  6 Acknowledgements

A.D. Bragg and X. Zhang acknowledge support through a National Science
Foundation (NSF) CAREER award # 2042346. This work used the Extreme Science
and Engineering Discovery Environment (XSEDE), which is supported by NSF grant
number ACI-1548562 (Towns et al., 2014). Specifically, the Comet and Expanse
clusters were used under allocation CTS170009. The DNS database was generated
with funding from the U.S. Office of Naval Research via grant
N00014-15-1-2248. High performance computing resources for the DNS were
provided through the U.S. Department of Defense High Performance Computing
Modernization Program by the Army Engineer Research and Development Center and
the Army Research Laboratory under Frontier Project FP-CFD-FY14-007.

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344.

##  7 Declaration of Interests

The authors report no conflict of interest.

## References

  * Alam et al. (2019) Alam, Shadab, Guha, Anirban & Verma, Mahendra K. 2019 Revisiting bolgiano–obukhov scaling for moderately stably stratified turbulence. Journal of Fluid Mechanics 875, 961–973.
  * Almalkie & de Bruyn Kops (2012) Almalkie, Saba & de Bruyn Kops, Stephen M. 2012 Kinetic energy dynamics in forced, homogeneous, and axisymmetric stably stratified turbulence. Journal of Turbulence 13, N29, arXiv: https://doi.org/10.1080/14685248.2012.702909.
  * Aluie & Eyink (2009) Aluie, Hussein & Eyink, Gregory L. 2009 Localness of energy cascade in hydrodynamic turbulence. ii. sharp spectral filter. Physics of Fluids 21 (11), 115108, arXiv: https://doi.org/10.1063/1.3266948.
  * Ayet et al. (2020) Ayet, A., Katul, G. G., Bragg, A. D. & Redelsperger, J. L. 2020 Scalewise return to isotropy in stratified boundary layer flows. Journal of Geophysical Research: Atmospheres 125 (16), e2020JD032732, e2020JD032732 10.1029/2020JD032732, arXiv: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020JD032732.
  * Betchov (1956) Betchov, R. 1956 An inequality concerning the production of vorticty in isotropic turbulence. J. Fluid Mech. 1, 497–504.
  * Billant & Chomaz (2001) Billant, Paul & Chomaz, Jean-Marc 2001 Self-similarity of strongly stratified inviscid flows. Physics of Fluids 13 (6), 1645–1651, arXiv: https://doi.org/10.1063/1.1369125.
  * Bolgiano Jr. (1959) Bolgiano Jr., R. 1959 Turbulent spectra in a stably stratified atmosphere. Journal of Geophysical Research (1896-1977) 64 (12), 2226–2229, arXiv: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/JZ064i012p02226.
  * Brethouwer et al. (2007) Brethouwer, G., Billant, P., Lindborg, E. & Chomaz, J.-M. 2007 Scaling analysis and simulation of strongly stratified turbulent flows. Journal of Fluid Mechanics 585, 343?368.
  * Brucker et al. (2007) Brucker, K. A., Isaza, J. C., Vaithianathan, T. & Collins, L. R. 2007 Efficient algorithm for simulating homogeneous turbulent shear flow without remeshing. J. Comp. Phys. 225, 20–32.
  * de Bruyn Kops (2015) de Bruyn Kops, S. M. 2015 Classical turbulence scaling and intermittency in stably stratified Boussinesq turbulence. J. Fluid Mech. 775, 436–463.
  * Buaria et al. (2020) Buaria, Dhawal, Bodenschatz, Eberhard & Pumir, Alain 2020 Vortex stretching and enstrophy production in high reynolds number turbulence. Phys. Rev. Fluids 5, 104602.
  * Carbone & Bragg (2020) Carbone, M. & Bragg, A. D. 2020 Is vortex stretching the main cause of the turbulent energy cascade? Journal of Fluid Mechanics 883, R2.
  * Chung & Matheou (2012) Chung, D. & Matheou, G. 2012 Direct numerical simulation of stationary homogeneous stratified sheared turbulence. Journal of Fluid Mechanics 696, 434–467.
  * Danish & Meneveau (2018) Danish, Mohammad & Meneveau, Charles 2018 Multiscale analysis of the invariants of the velocity gradient tensor in isotropic turbulence. Phys. Rev. Fluids 3, 044604.
  * Davidson (2004) Davidson, P. A. 2004 Turbulence: an introduction for scientists and engineers. Oxford.
  * Doan et al. (2018) Doan, N. A. K., Swaminathan, N., Davidson, P. A. & Tanahashi, M. 2018 Scale locality of the energy cascade using real space quantities. Phys. Rev. Fluids 3, 084601.
  * Ferrari & Wunsch (2009) Ferrari, Raffaele & Wunsch, Carl 2009 Ocean circulation kinetic energy: Reservoirs, sources, and sinks. Annual Review of Fluid Mechanics 41 (1), 253–282.
  * Gage (1979) Gage, K. S. 1979 Evidence far a k?5/3 law inertial range in mesoscale two-dimensional turbulence. Journal of the Atmospheric Sciences 36 (10), 1950–1954.
  * Germano (1992) Germano, M. 1992 Turbulence: the filtering approach. Journal of Fluid Mechanics 238, 325–336.
  * Gregg et al. (2018) Gregg, M.C., D’Asaro, E.A., Riley, J.J. & Kunze, E. 2018 Mixing efficiency in the ocean. Annual Review of Marine Science 10 (1), 443–473, pMID: 28934598.
  * Jayne (2009) Jayne, Steven R. 2009 The Impact of Abyssal Mixing Parameterizations in an Ocean General Circulation Model. Journal of Physical Oceanography 39 (7), 1756–1775.
  * Johnson (2020) Johnson, Perry L. 2020 Energy transfer from large to small scales in turbulence by multiscale nonlinear strain and vorticity interactions. Phys. Rev. Lett. 124, 104501.
  * Johnson (2021) Johnson, Perry L. 2021 On the role of vorticity stretching and strain self-amplification in the turbulence energy cascade. Journal of Fluid Mechanics 922, A3.
  * Katul et al. (2013) Katul, Gabriel G., Porporato, Amilcare, Manes, Costantino & Meneveau, Charles 2013 Co-spectrum and mean velocity in turbulent boundary layers. Physics of Fluids 25 (9), 091702, arXiv: https://doi.org/10.1063/1.4821997.
  * Khani & Waite (2015) Khani, Sina & Waite, Michael L. 2015 Large eddy simulations of stratified turbulence: the dynamic smagorinsky model. Journal of Fluid Mechanics 773, 327–344.
  * Kolmogorov (1941) Kolmogorov, A. N. 1941 The local structure of turbulence in an incompressible viscous fluid for very large Reynolds numbers. Dokl. Akad. Nauk. SSSR 30, 299–303.
  * Kumar et al. (2014) Kumar, Abhishek, Chatterjee, Anando G. & Verma, Mahendra K. 2014 Energy spectrum of buoyancy-driven turbulence. Phys. Rev. E 90, 023016.
  * Lee et al. (1990) Lee, Moon Joo, Kim, John & Moin, Parviz 1990 Structure of turbulence at high shear rate. Journal of Fluid Mechanics 216, 561?583.
  * Lilly (1983) Lilly, D. K. 1983 Stratified turbulence and the mesoscale variability of the atmosphere. Journal of the Atmospheric Sciences 40 (3), 749–761.
  * Lindborg (2005) Lindborg, Erik 2005 The effect of rotation on the mesoscale energy cascade in the free atmosphere. Geophysical Research Letters 32 (1).
  * Lindborg (2006) Lindborg, Erik 2006 The energy cascade in a strongly stratified fluid. Journal of Fluid Mechanics 550, 207?242.
  * L’vov & Falkovich (1992) L’vov, Victor & Falkovich, Gregory 1992 Counterbalanced interaction locality of developed hydrodynamic turbulence. Phys. Rev. A 46, 4762–4772.
  * Meneveau (2011) Meneveau, Charles 2011 Lagrangian dynamics and models of the velocity gradient tensor in turbulent flows. Annual Review of Fluid Mechanics 43 (1), 219–245.
  * Obukhov (1959) Obukhov, A 1959 Effect of archimedean forces on the structure of the temperature field in a turbulent flow. In Dokl. Akad. Nauk SSSR, , vol. 125, pp. 1246–1248.
  * Peltier & Caulfield (2003) Peltier, W. R. & Caulfield, C. P. 2003 Mixing efficiency in stratified shear flows. Annual Review of Fluid Mechanics 35 (1), 135–167.
  * Pope (2000) Pope, S. B. 2000 Turbulent Flows. New York: Cambridge University Press.
  * Portwood et al. (2019) Portwood, G. D., de Bruyn Kops, S. M. & Caulfield, C. P. 2019 Asymptotic dynamics of high dynamic range stratified turbulence. Phys. Rev. Lett. 122, 194504.
  * Riley & de Bruyn Kops (2003) Riley, James J. & de Bruyn Kops, Stephen M. 2003 Dynamics of turbulence strongly influenced by buoyancy. Physics of Fluids 15 (7), 2047–2059.
  * Riley & Lelong (2000) Riley, James J. & Lelong, Marie-Pascale 2000 Fluid motions in the presence of strong stable stratification. Annual Review of Fluid Mechanics 32 (1), 613–657.
  * Riley & Lindborg (2012) Riley, James J. & Lindborg, Erik 2012 Recent Progress in Stratified Turbulence, p. 269–317. Cambridge University Press.
  * Sekimoto et al. (2016) Sekimoto, Atsushi, Dong, Siwei & Jiménez, Javier 2016 Direct numerical simulation of statistically stationary and homogeneous shear turbulence and its relation to other shear flows. Physics of Fluids 28 (3), 035101, arXiv: https://doi.org/10.1063/1.4942496.
  * Sujovolsky & Mininni (2020) Sujovolsky, N. E. & Mininni, P. D. 2020 From waves to convection and back again: The phase space of stably stratified turbulence. Phys. Rev. Fluids 5, 064802.
  * Taylor (1922) Taylor, G. I. 1922 Diffusion by continuous movements. Proc. Lond. Math. Soc. 20, 196–212.
  * Taylor (1938) Taylor, G. I. 1938 Production and dissipation of vorticity in a turbulent fluid. Proc. Roy. Soc. Lond. A Mat. 164 (916), 15–23.
  * Taylor et al. (2016) Taylor, J. R., Deusebio, E., Caulfield, C. P. & Kerswell, R. R. 2016 A new method for isolating turbulent states in transitional stratified plane couette flow. Journal of Fluid Mechanics 808, R1.
  * Tennekes & Lumley (1972) Tennekes, H. & Lumley, J. L. 1972 A first course in turbulence. Cambridge: MIT Press.
  * Tom et al. (2020) Tom, Josin, Carbone, Maurizio & Bragg, Andrew D 2020 Exploring the turbulent velocity gradients at different scales from the perspective of the strain-rate eigenframe, arXiv: 2005.04300.
  * Towns et al. (2014) Towns, J., Cockerill, T., Dahan, M., Foster, I., Gaither, K., Grimshaw, A., Hazlewood, V., Lathrop, S., Lifka, D., Peterson, G. D., Roskies, R., Scott, J. R. & Wilkins-Diehr, N. 2014 Xsede: Accelerating scientific discovery. Computing in Science & Engineering 16 (5), 62–74.
  * Tsinober (2001) Tsinober, Arkady 2001 An informal introduction to turbulence. Kluwer Academic Publishers.
  * Vallis (2006) Vallis, G. K. 2006 Atmospheric and Oceanic Fluid Dynamics. Cambridge, U.K.: Cambridge University Press.
  * Vreman et al. (1994) Vreman, Bert, Geurts, Bernard & Kuerten, Hans 1994 Realizability conditions for the turbulent stress tensor in large-eddy simulation. Journal of Fluid Mechanics 278, 351–362.
  * Waite & Bartello (2006) Waite, Michael L. & Bartello, Peter 2006 The transition from geostrophic to stratified turbulence. Journal of Fluid Mechanics 568, 89?108.
  * Wyngaard (2010) Wyngaard, J. C. 2010 Turbulence in the Atmosphere. Cambridge: Cambridge University Press.
  * Zorzetto et al. (2018) Zorzetto, E., Bragg, A. D. & Katul, G. 2018 Extremes, intermittency, and time directionality of atmospheric turbulence at the crossover from production to inertial scales. Phys. Rev. Fluids 3, 094604.