# Effective Drift Velocity from Turbulent Transport by Vorticity

Hussein Aluie<EMAIL_ADDRESS>Department of Mechanical Engineering,
University of Rochester, Rochester, New York  Laboratory for Laser Energetics,
University of Rochester, Rochester, New York  Shikhar Rai  Department of
Mechanical Engineering, University of Rochester, Rochester, New York  Hao Yin
Department of Mechanical Engineering, University of Rochester, Rochester, New
York  Aarne Lees  Department of Mechanical Engineering, University of
Rochester, Rochester, New York  Laboratory for Laser Energetics, University of
Rochester, Rochester, New York  Dongxiao Zhao  Department of Mechanical
Engineering, University of Rochester, Rochester, New York  Stephen M. Griffies
NOAA/Geophysical Fluid Dynamics Laboratory, Princeton, New Jersey  Princeton
University Program in Atmospheric and Oceanic Sciences, Princeton, New Jersey
Alistair Adcroft  NOAA/Geophysical Fluid Dynamics Laboratory, Princeton, New
Jersey  Princeton University Program in Atmospheric and Oceanic Sciences,
Princeton, New Jersey  Jessica K. Shang  Department of Mechanical Engineering,
University of Rochester, Rochester, New York  Laboratory for Laser Energetics,
University of Rochester, Rochester, New York

###### Abstract

We highlight the differing roles of vorticity and strain in the transport of
coarse-grained scalars at length-scales larger than $\ell$ by smaller scale
(subscale) turbulence. We use the first term in a multiscale gradient
expansion due to Eyink Eyink06a , which exhibits excellent correlation with
the exact subscale physics when the partitioning length $\ell$ is any scale
smaller than that of the spectral peak. We show that unlike subscale strain,
which acts as an anisotropic diffusion/anti-diffusion tensor, subscale
vorticity’s contribution is solely a conservative advection of coarse-grained
quantities by an eddy-induced non-divergent velocity, ${\bf v}_{*}$, that is
proportional to the curl of vorticity. Therefore, material (Lagrangian)
advection of coarse-grained quantities is accomplished not by the coarse-
grained flow velocity, $\overline{\bf u}_{\ell}$, but by the effective
velocity, $\overline{\bf u}_{\ell}+{\bf v}_{*}$, the physics of which may
improve commonly used LES models.

##  I Introduction

Basic considerations from fluid dynamics (e.g. kundu2015fluid ) indicate that
the distance between particles in a laminar flow is determined by the strain.
Vorticity merely imparts a rotation on their separation vector ${\bf r}$
without affecting its magnitude. This behavior can be seen by considering the
velocity, ${\bf u}$, difference between particles $P$ and $Q$ at positions
${\bf x}$ and ${\bf x}+{\bf r}$, respectively,

$\displaystyle{\bf u}_{Q}-{\bf u}_{P}=\delta{\bf u}={\bf u}({\bf x}+{\bf
r})-{\bf u}({\bf x})={\bf
r}{\mbox{\boldmath$\cdot$}}{\mbox{\boldmath$\nabla$}}{\bf u}{\big{|}}_{{\bf
x}}+\dots~{},$ (1)

where a Taylor-series expansion is justified for short distances $|{\bf r}|$
over which the flow is sufficiently smooth. In the Lagrangian frame of $P$ at
${\bf x}$, the separation from $Q$ evolves as

$\displaystyle\frac{D{\bf r}}{Dt}=\delta{\bf u}={\bf
r}{\mbox{\boldmath$\cdot$}}{\bf S}+\underbrace{~{}{\bf
r}{\mbox{\boldmath$\cdot$}}{\mbox{\boldmath$\Omega$}}~{}}_{\frac{1}{2}{\mbox{\boldmath$\omega$}}{\mbox{\boldmath$\times$}}{\bf
r}},$ (2)

where the velocity gradient tensor, ${\mbox{\boldmath$\nabla$}}{\bf u}={\bf
S}+{\mbox{\boldmath$\Omega$}}$, has been decomposed into the symmetric strain
rate tensor ${\bf S}=[{\mbox{\boldmath$\nabla$}}{\bf
u}+({\mbox{\boldmath$\nabla$}}{\bf u})^{T}]/2$ and the antisymmetric vorticity
tensor ${\mbox{\boldmath$\Omega$}}=[{\mbox{\boldmath$\nabla$}}{\bf
u}-({\mbox{\boldmath$\nabla$}}{\bf
u})^{T}]/2=-\frac{1}{2}\epsilon_{ijk}\omega_{k}$. Here,
${\mbox{\boldmath$\omega$}}={\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\times$}}{\bf
u}$ is vorticity and $\epsilon_{ijk}$ is the Levi-Civita symbol. Taking an
inner product of eq. (2) with ${\bf r}$,

$\displaystyle\frac{1}{2}\frac{D|{\bf r}|^{2}}{dt}={\bf
r}{\mbox{\boldmath$\cdot$}}{\bf S}{\mbox{\boldmath$\cdot$}}{\bf r}~{},$ (3)

shows that the distance is determined by the strain. Vorticity in eq. (2) only
acts to rotate ${\bf r}$ without changing its magnitude.

These considerations hinge on the critical assumption that the flow is
sufficiently smooth over separations ${\bf r}$, which is patently invalid in a
turbulent flow for ${\bf r}$ at inertial scales Pope00 . However, a version of
this story survives thanks to the property of scale-locality, which justifies
an expansion in scale as we shall discuss below. The main result of this paper
is eq. (43), which is an expression for the eddy-induced advection velocity
${\bf v}_{*}$ at length-scales larger than $\ell$, which may be the size of
grid cells in a simulation. The non-divergent velocity, ${\bf v}_{*}$, arises
from vortical motions associated with subscale nonlinear interactions. Fig. 1
captures the essential insight behind eq. (43).

Figure 1: Schematic of how eddy-induced advection by a non-divergent velocity,
${\bf
v}_{*}\sim{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\times$}}\overline{{\mbox{\boldmath$\omega$}}}_{\ell}$,
(blue arrow) from eq. (43) can arise due to subscale vorticity (red). In this
example, two symmetric but counter-rotating eddies (red) are of a size smaller
than that of the coarse-graining box (dashed). At length-scales larger than
that of the box, $\ell$, which may represent a grid cell in a simulation or
the resolution limit of an observation, these eddies exert an eddy-induced
advection with velocity ${\bf v}_{*}$ proportional to the curl of vorticity,
depicted by the blue arrow at location ${\bf x}$.

In section II, we discuss subscale transport, Eyink’s expansion in length-
scale, and the leading order term. In section III, we present empirical
support from (i) 3D direct numerical simulation (DNS) of a compressible
turbulent flow and (ii) two-layer stacked shallow water simulation of
geophysical fluid turbulence. Section IV contains the main result, which
expresses transport by subscale vorticity as a conservative advection by an
effective drift velocity ${\bf v}_{*}$ proportional to the curl of vorticity.
Throughout this paper we attempt to make the presentation accessible to
readers from various backgrounds.

##  II Multi-Scale Dynamics

To analyze the dynamics of different scales in a flow, we use the coarse-
graining approach, which has proven to be a natural and versatile framework to
understand and model scale interactions (e.g. MeneveauKatz00 ; Eyink05 ). The
approach is standard in partial differential equations and distribution theory
(e.g. Strichartz03 ; Evans10 ). It became common in Large Eddy Simulation
(LES) modeling of turbulence thanks to the foundational works of Leonard
Leonard74  and Germano Germano92 .

For any field ${\bf a}({\bf x})$, a coarse-grained or (low-pass) filtered
version of this field, which contains spatial variations at scales $>\ell$, is
defined in $n$-dimensions as

$\overline{\bf a}_{\ell}({\bf x})=\int\mathrm{d}^{n}{\bf r}~{}G_{\ell}({\bf
x}-{\bf r})\,{\bf a}({\bf r}).$ (4)

Here, $G_{\ell}$ is a convolution kernel derived from a “mother” kernel
$G({\bf s})$ (borrowing the term from wavelet analysis daubechies1992ten ).
$G({\bf s})$ is normalized, $\int\mathrm{d}^{n}{\bf s}\,G({\bf s})=1$. It is
an even function such that, $\int{\bf s}\,G({\bf s})\,\mathrm{d}^{n}{\bf
s}=0$, which ensures that local averaging is symmetric. $G({\bf s})$ also has
its main support (or variance) over a region of order unity in diameter,
$\int\mathrm{d}^{n}\,{\bf s}\,|{\bf s}|^{2}\,G({\bf s})=O(1)$. The dilated
version of the kernel,

$G_{\ell}({\bf r})=\ell^{-n}G({\bf r}/\ell),$ (5)

inherits all those properties, except that its main support is over a region
of diameter $\ell$. An example is the Top-Hat kernel in non-dimensional
coordinates ${\bf s}$,

$\displaystyle G({\bf s})=\begin{cases}1,&\text{if $|s_{i}|<1/2$ for
$i=1,\dots,n$},\\\ 0,&\text{otherwise},\\\ \end{cases}$ (6)

and its dilated version in dimensional coordinates ${\bf r}$,

$\displaystyle G_{\ell}({\bf r})=\begin{cases}\ell^{-n},&\text{if
$|r_{i}|<\ell/2$ for $i=1,\dots,n$},\\\ 0,&\text{otherwise}.\\\ \end{cases}$
(7)

The scale decomposition in eq. (4) is essentially a partitioning of scales for
a field into large scales ($\gtrsim\ell$), captured by $\overline{\bf
a}_{\ell}$, and small scales ($\lesssim\ell$), captured by the residual

${\bf a}^{\prime}_{\ell}={\bf a}-\overline{\bf a}_{\ell}.$ (8)

In the remainder of this paper, we shall omit subscript $\ell$ from variables
if there is no risk for confusion.

###  II.1 Coarse-grained Transport

Consider the transport equation of the concentration (per unit volume) $C({\bf
x},t)$:

$\partial_{t}C+{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}(C\,{\bf
u})={\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}[\rho\,\kappa\,{\mbox{\boldmath$\nabla$}}\phi].$
(9)

Here, $\rho$ is mass density, ${\bf u}$ is the advecting flow velocity, which
is possibly compressible. The scalar field $C=\rho\,\phi$, where $\phi$ is
concentration per unit mass, can be either passive or active, with diffusivity
$\kappa$, which can be due to the microphysics or an effective (turbulent)
diffusivity used in a simulation. Analysis of the incompressible momentum
transport follows similar reasoning as below, with $C$ replaced by ${\bf u}$
when $\rho=({\rm const.})$, so we shall focus on $C$ for simplicity.

We shall call eq. (9) the “bare” dynamics Aluie17 ; eyink2018review  as it
governs flow in an experiment or natural system. In a coarse-resolution
simulation, such as in implicit or explicit LES grinstein2007implicit ;
MeneveauKatz00 , or when taking spatially under-resolved measurements in an
experiment, ${\bf u}$ and $C$ are not resolved down to the microphysical
scales. The equation that governs the coarse-grained scalar
$\overline{C}_{\ell}$ is obtained by filtering eq. (9) to get (e.g.
higgins2004heat )

$\displaystyle\partial_{t}\overline{C}_{\ell}+{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}(\overline{{\bf
u}}_{\ell}\,\overline{C}_{\ell})=-{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}\overline{\tau}_{\ell}+{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}\overline{(\rho\,\kappa{\mbox{\boldmath$\nabla$}}\phi)}_{\ell}~{},$
(10a) where $\displaystyle\overline{\tau}_{\ell}=\overline{\tau}_{\ell}({\bf
u},C)\equiv\overline{({\bf u}\,C)}_{\ell}-\overline{{\bf
u}}_{\ell}\,\overline{C}_{\ell}$ (10b)

is the subscale flux of $\overline{C}_{\ell}$ Germano92 . Eq.(10) is exact,
without any approximation, and describes the scalar transport at length-scales
larger than $\ell$. The last term,
${\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}\overline{(\rho\,\kappa{\mbox{\boldmath$\nabla$}}\phi)}_{\ell}$,
was shown mathematically in Aluie13  to have an upper bound proportional to
$\alpha_{\textrm{rms}}\,\ell^{-2}\,\phi_{\textrm{rms}}$ at every location in
the flow and every time, where $\alpha=\rho\,\kappa$. Therefore, it is
guaranteed to be negligible for sufficiently large $\ell$ or sufficiently
small $\alpha$, which was demonstrated numerically in ZhaoAluie18 . Therefore,
${\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}\overline{(\rho\,\kappa{\mbox{\boldmath$\nabla$}}\phi)}_{\ell}$
can be dropped from eq. (10) (but not from the bare eq. (9)).

An evident role of the subscale flux $\overline{\tau}_{\ell}(C,{\bf u})$ in
eq. (10) is accounting for the influence of unresolved scales
$C^{\prime}_{\ell}$ and ${\bf u}^{\prime}_{\ell}$ on the transport of
$\overline{C}_{\ell}$. Another important role of
$\overline{\tau}_{\ell}(C,{\bf u})$ in eq. (10) that may not be as well-
appreciated is balancing $\overline{{\bf u}}_{\ell}\,\overline{C}_{\ell}$,
which can interact with scales $<\ell$ that are approximately within the band
$[\ell/2,\ell]$. This is analogous to the term
$\rho\kappa{\mbox{\boldmath$\nabla$}}\phi$ balancing ${\bf u}\,C$ at
microphysical scales $\sim\ell_{\kappa}$, below which molecular diffusion
dominates advection. To illustrate, assume that at a time $t$ the scalar and
velocity fields are fully resolved, $\overline{C}_{\ell}=C$ and
$\overline{{\bf u}}_{\ell}={\bf u}$. Moreover, assume that scale $\ell$ is
sufficiently large, $\ell\gg\ell_{\kappa}$, such that microphysical
diffusivity is negligible and does not damp nonlinear interactions. That small
scales are absent, $C^{\prime}_{\ell}=0$ and ${\bf u}^{\prime}_{\ell}=0$, does
not necessarily imply that the subscale flux $\overline{\tau}_{\ell}(C,{\bf
u})=\overline{({\bf u}\,C)}_{\ell}-{\bf u}\,C$ is zero unless there is
sufficient scale separation between $\ell$ and the smallest scales of $C$ and
${\bf u}$. Therefore, the general transport of a coarse-grained scalar
$\overline{C}_{\ell}$ is not described by the typical (or bare) transport
equation (9) even when $C$ and ${\bf u}$ are resolved ($C^{\prime}_{\ell}=0$
and ${\bf u}^{\prime}_{\ell}=0$), because coarse-graining also alters the
representation of their nonlinear interactions, which may be at scales smaller
than $\ell$.

From definition (10b), it is possible to decompose
$\overline{\tau}_{\ell}(C,{\bf u})$ into contributions that are solely due to
resolved fields $\overline{C}_{\ell}$ and $\overline{{\bf u}}_{\ell}$ (e.g.
the Leonard term Leonard74 ; germano1986proposal ) and other contributions
that involve unresolved fields $C^{\prime}_{\ell}$ and ${\bf
u}^{\prime}_{\ell}$. In the discussion below, we use the term ‘subscale’ to
refer to any contribution from $\overline{\tau}_{\ell}(C,{\bf u})$ because it
arises from coarse-graining the dynamics at scale $\ell$. We shall see that
the resolved fields make a dominant contribution to the subscale flux
$\overline{\tau}_{\ell}(C,{\bf u})$ thanks to the property of ultraviolet
scale-locality.

###  II.2 Multiscale Expansion

The subscale term $\overline{\tau}_{\ell}({\bf u},C)$ does not lend itself to
a straightforward intuitive interpretation. To that end, we utilize the
multiscale gradient expansion of Eyink Eyink06a . This expansion is
accomplished by decomposing a field

${\bf a}({\bf x})=\sum_{p=0}^{\infty}{\bf a}^{[p]}({\bf x})$ (11)

into a sum of band-passed fields Eyink06a

$\displaystyle{\bf a}^{[p]}$ $\displaystyle=\overline{\bf
a}_{\ell_{p}}-\overline{\bf a}_{\ell_{p-1}},~{}~{}~{}p=1,2,\dots,$ (12a)
$\displaystyle{\bf a}^{[0]}$ $\displaystyle=\overline{\bf a}_{\ell_{0}}$ (12b)

for a sequence of scales $\ell_{p}=\lambda^{-p}\ell_{0}$ with the non-
dimensional number $\lambda>1$, e.g. $\lambda=2$. $\ell_{0}$ can be a
characteristic large scale such as the domain size. A similar decomposition
can be done for the subscale flux,

$\displaystyle\overline{\tau}_{\ell}({\bf
u},C)=\sum_{p=0}^{\infty}\sum_{q=0}^{\infty}\tau_{\ell}^{[p,q]}({\bf u},C),$
(13a) with $\displaystyle\tau_{\ell}^{[p,q]}({\bf u},C)=\overline{({\bf
u}^{[p]}C^{[q]})}_{\ell}-\overline{({\bf
u}^{[p]})}_{\ell}\overline{(C^{[q]})}_{\ell}~{}.$ (13b)

The term $\tau_{\ell}^{[p,q]}({\bf u},C)$ represents the subscale contribution
from scalar and velocity scales within bands $[q]$ and $[p]$, respectively.
Eyink Eyink05 ; Eyink06a  proved that the series (13) converges absolutely and
at a geometric rate under the condition that the scaling exponents of
increments of $\delta{\bf u}(r)\sim r^{\sigma^{u}}$ and $\delta C(r)\sim
r^{\sigma^{C}}$ satisfy $\sigma^{u}>0$ and $\sigma^{C}>0$ 111More precisely,
the condition is that $\sigma^{u}_{p}>0$ and $\sigma^{C}_{p}>0$ if increments
scale as $\langle|\delta{\bf u}(r)|^{p}\rangle^{1/p}\sim
u_{\text{rms}}A_{p}(r/\ell_{0})^{\sigma^{u}_{p}}$ and $\langle|\delta
C(r)|^{p}\rangle^{1/p}\sim C_{\text{rms}}B_{p}(r/\ell_{0})^{\sigma^{C}_{p}}$
for some dimensionless constants $A_{p}$ and $B_{p}$, where
$\langle\dots\rangle$ is a space average. The condition has to hold for at
least $p=2$, with higher values of $p$ representing more stringent conditions.
Note that a power-law scaling of increments is not required; the condition is
still satisfied if increments decay faster than any power-law with $r\to 0$,
for example exponentially..

Exponents $\sigma^{u}$ and $\sigma^{C}$ reflect the smoothness of their
respective fields (see Fig. 1 in Aluie17  and the associated discussion), with
larger values of $\sigma$ implying smoother (or, heuristically, less
turbulent) fields. For example, the Kolmogorov-Obukhov theory
kolmogorov1941local ; obukhov1949structure ; corrsin1951spectrum  predicts
$\sigma^{u}=\sigma^{C}=1/3$ for velocity and passive scalar increments over
scales within the inertial–convective range of a high Reynolds number
turbulent flow that is statistically homogeneous. In contrast, for low
Reynolds number flows, even if they are chaotic, we can expect $\sigma^{u}\geq
1$. Ignoring intermittency corrections, the scaling exponent $\sigma$ of
increments is related to the power-law scaling of a spectrum, $E(k)$, of
either ${\bf u}$ or $C$, via $E(k)\sim k^{-2\sigma-1}$. Therefore, the
condition $\sigma>0$ corresponds to a wavenumber scaling of spectra that
decays faster than $k^{-1}$ in wavenumber, which is a fairly weak condition
and is expected to hold in most flows. The conditions are the same as those
required for $\overline{\tau}_{\ell}({\bf u},C)$ to be ultraviolet scale-local
Eyink05 ; eyink2009localness ; aluie2009localness , which means that length-
scales $\delta\ll\ell$ make a negligible contribution to
$\overline{\tau}_{\ell}({\bf u},C)$,

$|\overline{\tau}_{\ell}(C^{\prime}_{\delta},{\bf
u}^{\prime}_{\delta})|\ll|\overline{\tau}_{\ell}(C,{\bf u})|.$ (14)

Under these conditions of convergence, the leading order term in eq. (13) is
the dominant contribution to the series, which justifies the approximation

$\overline{\tau}_{\ell}({\bf u},C)\approx\tau_{\ell}^{[0,0]}({\bf
u},C)=\overline{\tau}_{\ell}(\overline{{\bf u}}_{\ell},\overline{C}_{\ell}).$
(15)

This approximation corresponds to the well-known similarity (or Bardina) model
Bardinaetal80 ; MeneveauKatz00 . To readers who are perhaps less familiar with
LES or coarse-graining, the term

$\overline{\tau}_{\ell}(\overline{{\bf
u}}_{\ell},\overline{C}_{\ell})=\overline{(\overline{{\bf
u}}_{\ell}\overline{C}_{\ell})}_{\ell}-\overline{(\overline{{\bf
u}}_{\ell})}_{\ell}\overline{(\overline{C}_{\ell})}_{\ell}$ (16)

may seem insignificant, especially considering that if one were to replace
operation $\overline{(\dots)}_{\ell}$ with a Reynolds (or ensemble) average,
expression (16) would be identically zero. However, operation
$\overline{(\dots)}_{\ell}$ is a decomposition in length-scale, which is
inherently different from Reynolds averaging. To aid in the conceptual
understanding of expression (16), assume that $\overline{(\dots)}_{\ell}$ is a
sharp-spectral projection in wavenumber space,
$\overline{(\dots)}_{\ell}=(\dots)^{<K}$. Applying $(\dots)^{<K}$ removes all
wavenumbers larger than $K=\ell^{-1}$ (i.e. small length-scales) while keeping
wavenumbers smaller than $K$ intact. Then

$\displaystyle\overline{\tau}_{\ell}(\overline{{\bf
u}}_{\ell},\overline{C}_{\ell})$ $\displaystyle=({\bf
u}^{<K}C^{<K})^{<K}-({\bf u}^{<K})^{<K}(C^{<K})^{<K}$ (17a)
$\displaystyle=({\bf u}^{<K}C^{<K})^{<K}-{\bf u}^{<K}C^{<K},$ (17b)

using the sharp projection property, $((\dots)^{<K})^{<K}=(\dots)^{<K}$, to
arrive at the last expression. The term ${\bf u}^{<K}C^{<K}$ in eq.(17b),
being a quadratic product, contains wavenumbers $<2K$. Of these, the
contribution $({\bf u}^{<K}C^{<K})^{<K}$ from wavenumbers $<K$ is subtracted
in eq.(17b), such that $\overline{\tau}_{\ell}(\overline{{\bf
u}}_{\ell},\overline{C}_{\ell})$ represents wavenumbers within the band
$[K,2K]$ or, equivalently, between length-scales $\ell$ and $\ell/2$. As
discussed in previous work on scale-locality Eyink05 ; Aluie11 , this band of
scales is expected to make the dominant contribution to the subscale flux
$\overline{\tau}_{\ell}({\bf u},C)$ if the spectrum decays faster than
$k^{-1}$ in wavenumber over the range of scales $<\ell$. This is essentially
the condition under which eq. (15) was derived. Physically, a spectral decay
faster than $k^{-1}$ implies that dyadic wavenumber band $[k,2k]$ has more
energy than band $[2k,4k]$ (e.g. (frazier1991littlewood, ; AluieEyink09, )),
which justifies retaining the smallest wavenumbers via the leading order term
in expansion (13) aluie2009localness . A spectral decay faster than $k^{-1}$
is a fairly weak condition and is expected in most flows (turbulent or
laminar), but may fail, for example, if $\ell$ approaches scales of the
spectral peak or larger. We emphasize that filtering with a sharp-spectral
kernel in eq. (17) is merely for conceptual understanding. Unlike the Top-Hat
or Gaussian, a sharp-spectral kernel is not positive semi-definite in physical
space, which violates physical realizability conditions Vremanetal94 ,
rendering subscale energy or coarse-grained mass density negative in physical
space Eyink05 ; Aluie13 .

###  II.3 Relation to Increments and Gradients

Using the usual definition of an increment,

$\delta f({\bf x};{\bf r})=f(\mathbf{x}+\mathbf{r})-f(\mathbf{x}),$ (18)

the subscale flux (10b) can be rewritten exactly in terms of $\delta C$ and
$\delta{\bf u}$ Constantinetal94 ; Eyink05 ; Aluie11 :

$\overline{\tau}(C,{\bf u})=\left\langle\delta{C}\,\delta{{\bf
u}}\right\rangle_{\ell}-\left\langle\delta
C\right\rangle_{\ell}\left\langle\delta{{\bf u}}\right\rangle_{\ell}.$ (19)

Equation (19) is exact and follows directly from eq. (10b) without requiring
any assumption about scale-locality or the scaling exponents, where

$\left\langle\delta{f}\right\rangle_{\ell}({\bf
x})\equiv\int\mathrm{d}^{n}{\bf r}\,G_{\ell}(-{\bf r})\delta{f}({\bf x};{\bf
r})$ (20)

is a local average around ${\bf x}$ over all separations ${\bf r}$ weighted by
the kernel $G_{\ell}$. A spatially localized (or compact) kernel effectively
limits the average in eq. (20) to separations $|{\bf r}|\leq\ell/2$. Note that
$G_{\ell}(-{\bf r})=G_{\ell}({\bf r})$ for even kernels, such as a Gaussian or
a Top-Hat. Relation (19) is a key step and reveals the connection between the
subscale physics and increments (spatial variations) over distances smaller
than $\ell$. This connection is underscored by noting that expression (20)
equals (minus) the high-pass field containing scales $<\ell$,

$\int\mathrm{d}^{n}{\bf r}\,G_{\ell}(-{\bf r})\,\delta{f}({\bf x};{\bf
r})=-[f({\bf x})-\overline{f}_{\ell}({\bf x})]=-f^{\prime}_{\ell}({\bf
x})~{}.$ (21)

We may also connect these considerations with the argument presented in the
introduction, where increments reflect gradients Eyink06a , including strain
and vorticity of the subscale flow.

Approximating the subscale flux by the leading order term in eq.(15) allows us
to use identity (19) for the filtered quantities,

$\overline{\tau}(\overline{C}_{\ell},\overline{{\bf
u}}_{\ell})=\left\langle\delta\overline{C}_{\ell}\,\delta\overline{{\bf
u}}_{\ell}\right\rangle_{\ell}-\left\langle\delta\overline{C}_{\ell}\right\rangle_{\ell}\left\langle\delta\overline{{\bf
u}}_{\ell}\right\rangle_{\ell}.$ (22)

Since a filtered field $\overline{f}_{\ell}({\bf x})$ is smooth, we can Taylor
expand its increments around ${\bf x}$,

$\delta\overline{f}_{\ell}({\bf x};{\bf r})=\overline{f}_{\ell}({\bf x}+{\bf
r})-\overline{f}_{\ell}({\bf
x})\approx\mathbf{r}{\mbox{\boldmath$\cdot$}}{\mbox{\boldmath$\nabla$}}\overline{f}_{\ell}({\bf
x})+\dots,$ (23)

where we neglect higher order terms. Substituting the first term in the Taylor
expansion of each of $\delta\overline{C}$ and $\delta\overline{{\bf u}}$ into
eq. (22) gives in $n$-dimensions

$\displaystyle\overline{\tau}_{\ell}(\overline{C},\overline{u}_{i})$
$\displaystyle=\left[\partial_{k}\overline{C}_{\ell}\right]\left[\partial_{m}\overline{(u_{i})}_{\ell}\right]\left[\left\langle
r_{k}\,r_{m}\right\rangle_{\ell}-\left\langle
r_{k}\right\rangle_{\ell}\left\langle r_{m}\right\rangle_{\ell}\right]$ (24a)
$\displaystyle=\left[\partial_{k}\overline{C}_{\ell}\right]\left[\partial_{m}\overline{(u_{i})}_{\ell}\right]\left[\frac{1}{n}\delta_{km}\,\ell^{2}\int\mathrm{d}^{n}\mathbf{s}\,G(\mathbf{s})\left|\mathbf{s}\right|^{2}\right]$
(24b)
$\displaystyle=\frac{1}{n}\,M_{2}\,\ell^{2}\,\partial_{k}\overline{(u_{i})}_{\ell}\,\partial_{k}\overline{C}_{\ell}~{}.$
(24c)

In deriving the second line, we used the symmetry of the kernel such that
$\left\langle r_{k}\right\rangle_{\ell}=0$. In the final expression, the
mother kernel’s second moment,
$M_{2}=\int\mathrm{d}^{n}\mathbf{s}\,G(\mathbf{s})\left|\mathbf{s}\right|^{2}$,
depends solely on the shape of kernel $G$ and, in particular, is independent
of scale $\ell$. For the Top-Hat in eq. (6), $M_{2}=n/12$ in $n$-dimensions.

###  II.4 Summary and Interpretation

From the derivation above, we see that the final expression in eq. (24c)
represents the leading order contribution to $\overline{\tau}_{\ell}(C,{\bf
u})$. Therefore, to the extent expression (24c) faithfully approximates
$\overline{\tau}_{\ell}(C,{\bf u})$, it also represents the subscale physics.
How can expression (24c), which involves only coarse-grained (or resolved)
quantities represent subscale interactions? It is analogous to inferring the
value of a function $f(x_{0}+h)$ a distance $h$ away from $x_{0}$ based solely
on its local properties at $x_{0}$ via a Taylor series expansion, the accuracy
of which depends on smoothness properties of $f$. Here, the expansion is done
in length-scale Eyink06a  rather than in space to arrive at eq. (24c). The
analogue to smoothness is ultraviolet scale-locality, which is guaranteed
under weak spectral scaling conditions Eyink05 ; eyink2009localness ;
aluie2009localness ; Aluie11 .

Expression (24c) corresponds to the nonlinear (or Clark) model
clark1979evaluation , which is a closure of the subscale flux using a product
of gradients in LES modeling and has been shown to be an excellent (_a
priori_) approximation of the subscales in several previous studies (e.g.
Liuetal94 ; porte2001priori ; higgins2004heat ). Some studies have also
offered viable ideas for utilizing it in a posteriori LES Bardinaetal80 ;
kosovic1997subgrid ; bouchet2003parameterization . Below, we provide
additional _a priori_ support from simulations of 3D compressible turbulence
and two-layer stacked shallow water simulations of geophysical fluid
turbulence.

In the discussion above, we followed Eyink’s multiscale gradient expansion
Eyink06a  to derive expression (24c), which is more general than the standard
derivations of the nonlinear model found in the LES literature (e.g.
clark1979evaluation ; Bardinaetal80 ; Liuetal94 ; leonard1997large ;
BorueOrszag98 ; Pope00 ; bouchet2003parameterization ). Since Eyink’s
multiscale gradient expansion is convergent, it provides a rigorous foundation
to justify neglecting higher order terms. Another advantage of using Eyink’s
expansion is in bringing to the fore the underlying physics, including the
assumptions under which the series (13) converges and, therefore, the
conditions under which the approximation may fail. Moreover, being a
systematic expansion, it allows for future improvements by retaining higher
order terms or by using a closed form of those terms following the recent work
of Johnson johnson2020energy ; johnson2021role .

Approximation (24c) provides us with physical insight into the mechanisms by
which the subscales transport $\overline{C}_{\ell}$. Below, we show that
subscale vorticity’s contribution is an eddy-induced drift velocity. Prior to
doing so, we discuss how well eq. (24c) approximates
$\overline{\tau}_{\ell}({\bf u},C)$ from simulations.

##  III Empirical Support

Approximating $\overline{\tau}_{\ell}({\bf u},C)$ by expression (24c) has
significant support in the LES literature from _a priori_ tests of the
nonlinear model. It has been shown to exhibit excellent agreement when $C$ is
taken as the momentum (e.g. Liuetal94 ; BorueOrszag98 ), as the temperature
field in the atmospheric boundary layer porte2001priori ; higgins2004heat ,
and as a passive scalar kang2001passive ; chumakov2008priori .

In this section, we present additional supporting evidence for active scalars,
when $C$ represents (i) density in a 3D compressible turbulent flow and (ii)
layer thickness in a two-layer stacked shallow water simulation of geophysical
fluid turbulence. The correlation is generally better than $0.9$ at inertial
scales smaller than the scale at which the spectrum peaks.

###  III.1 Compressible Turbulence

We carry out a direct numerical simulation (DNS) of forced compressible
turbulence in a periodic box of non-dimensional size $2\pi$ with grid
resolution $1024^{3}$. From the simulation data, we test how well the
expression in eq. (24c) approximates $\overline{\tau}_{\ell}({\bf u},C)$ with
$C$ replaced by the density field, $\rho$, which is an active scalar.

Figure 2: Left to right: visualization of x, y and z components, respectively,
of momentum, $\rho\,{\bf u}$, in physical space from the $1024^{3}$ DNS. These
images are shown at an instant of time after the flow has reached steady
state. Shocks can be seen as discontinuities. The spectrum of this flow is
shown in Fig. 15 in the Appendix.

The DNS follows the configuration of lees2019baropycnal  and solves the fully
compressible Navier-Stokes equations using our DiNuSUR code ZhaoAluie18 ;
lees2019baropycnal :

$\displaystyle\partial_{t}\rho$ $\displaystyle+\partial_{j}(\rho\,u_{j})=0,$
(25) $\displaystyle\partial_{t}(\rho\,u_{i})$
$\displaystyle+\partial_{j}(\rho\,u_{i}\,u_{j})=-\partial_{i}P+\partial_{j}\sigma_{ij}+\rho\,F_{i}~{},$
(26) $\displaystyle\partial_{t}(\rho\,E)$
$\displaystyle+\partial_{j}(\rho\,E\,u_{j})=-\partial_{j}(P\,u_{j})+\partial_{j}[2\,\mu\,u_{i}\,(S_{ij}-\frac{1}{3}S_{kk}\,\delta_{ij})]-\partial_{j}q_{j}+\rho\,u_{i}\,F_{i}-\mathcal{RL}~{}.$
(27)

Here, $E=|{\bf u}|^{2}/2+e$ is total energy per unit mass, where $e$ is
specific internal energy, $P$ is thermodynamic pressure, $\mu$ is dynamic
viscosity, ${\bf q}=-\kappa{\mbox{\boldmath$\nabla$}}T$ is the heat flux with
a thermal conductivity $\kappa$ and temperature $T$. Both dynamic viscosity
and thermal conductivity are spatially variable, where $\mu({\bf
x})=\mu_{0}(T({\bf x})/T_{0})^{0.76}$. Thermal conductivity is set to satisfy
a Prandtl number $Pr=c_{p}\,\mu/\kappa=0.7$, where
$c_{p}=R\,\gamma/(\gamma-1)$ is the specific heat with specific gas constant
$R$ and $\gamma=5/3$ for a monoatomic gas. We use the ideal gas equation of
state, $P=\rho\,R\,T$. We stir the flow using an external acceleration field
$F_{i}$, and $\mathcal{RL}=\rho\,u_{i}\,F_{i}$ represents radiation losses
from internal energy (see JagannathanDonzis16 ; lees2019baropycnal ).
$S_{ij}=(\partial_{j}u_{i}+\partial_{i}u_{j})/2$ is the symmetric rate of
strain tensor and $\sigma_{ij}$ is the deviatoric (traceless) viscous stress,

$\sigma_{ij}=2\,\mu\,(S_{ij}-\frac{1}{3}S_{kk}\,\delta_{ij}).$ (28)

We solve the above equations using the pseudo-spectral method with $2/3$rd
dealiasing. We advance in time using the $4^{th}$-order Runge-Kutta scheme
with a variable time step.

The acceleration, ${\bf F}$, is similar to that in Federrathetal08 . In
Fourier space, the acceleration is defined as

$\widehat{F}_{i}(\mathbf{k})=\widehat{f}_{j}(\mathbf{k})\,P_{ij}^{\zeta}(\mathbf{k}),$
(29)

where $\mathbf{\widehat{f}}$ is constructed from independent Ornstein-
Uhlenbeck stochastic processes EswaranPope88 . The projection operator

$P_{ij}^{\zeta}(\mathbf{k})=\zeta\,\delta_{ij}+(1-2\,\zeta)\,\frac{k_{i}\,k_{j}}{|\mathbf{k}|^{2}}$
(30)

allows for controlling the ratio of solenoidal
(${\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}\mathbf{F}=0$) to
dilatational
(${\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\times$}}\mathbf{F}=0$)
components of the forcing using the parameter $\zeta$. When $\zeta=0$, the
forcing is purely dilatational and when $\zeta=1$, the forcing is purely
solenoidal. For the simulation we use here, shown in Fig. 2, we set
$\zeta=0.01$ to simulate turbulence with high compressibility, such that
dilatation, ${\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}{\bf u}\neq
0$, is significant.

Our simulated flow has significant compressibility, which can be seen from
Fig. 2 and inferred from the following three time-averaged quantities. (1) The
ratio $({\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}{\bf
u})_{rms}/({\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\times$}}{\bf
u})_{rms}=5.61$ in our flow, which measures compressibility at small scales.
(2) We also Helmholtz decompose the velocity field into the dilatational,
${\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\times$}}{\bf u}^{d}=0$, and
solenoidal, ${\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}{\bf
u}^{s}=0$, components, ${\bf u}={\bf u}^{d}+{\bf u}^{s}$. The ratio of
dilatational kinetic energy,
$K^{d}=\langle\rho\,u^{d}_{i}\,u^{d}_{i}/2\rangle$, to solenoidal kinetic
energy, $K^{s}=\langle\rho\,u^{s}_{i}\,u^{s}_{i}/2\rangle$, yields a measure
of compressibility at large scales, which in our flow is $K^{d}/K^{s}=4.80$.
(3) We also calculate the turbulence Mach number, $M_{t}=\langle|{\bf
u}|^{2}\rangle^{1/2}/c=0.23$, where $c=2.18$ is the mean non-dimensional sound
speed. Spectra of density and velocity are plotted in Fig. 15 in the Appendix.

We use the Top-Hat kernel in eq. (7) for the coarse-graining. Fig. 3 shows a
2D slice from the 3D flow at the instant of time shown in Fig. 2, comparing
the exact $\overline{\tau}_{\ell}(u_{x},\rho)$ at $\ell=0.19635$ (domain size
is $2\pi$) with its approximation (24c). The other two components are shown in
the Appendix Fig. 16. All indicate an excellent agreement. Fig. 4 (left
column) shows a more quantitative evaluation of all three components of the
subscale transport approximation, where we can see that the approximation is
remarkably accurate almost everywhere but underestimates
$\overline{\tau}_{\ell}(u_{x},\rho)$ at locations of strong shocks. This
underestimation indicates that higher order terms in the expansion (13) are
needed for better agreement at those locations since these rare but strong
discontinuities have a significant contribution from scales $<\ell$ and are an
extreme test case for the approximation. Indeed, at a location $x_{0}$ of a
shock, the increment scales as $\delta u(x_{0};r)=u(x_{0}+r)-u(x_{0})\sim
r^{\sigma}$, with $\sigma\approx 0$ for distances $r$ larger than the viscous
shock width (see Fig. 1 in Aluie17  and the associated discussion). This value
is at the margin of validity of the $\sigma>0$ condition 222Formally, the
condition $\sigma>0$ for Eyink’s expansion to converge is not based on the
scaling of an increment $\delta u(x;r)$ at a single point $x$. Rather, it is
based on the scaling of structure functions, which are spatial averages of
increments. for Eyink’s expansion (13) to converge.

Fig. 4 (right column) shows the joint probability density function (PDF)
between $\overline{\tau}_{\ell}({\bf u},\rho)$ and its approximation, with
correlations exceeding $0.9$. In Figs. 11-13 in the Appendix, we show similar
comparisons at other length-scales. The correlation coefficient between two
datasets, $A$ and $B$, is defined as

$r(A,B)=\frac{\langle A\,B\rangle-\langle A\rangle\langle B\rangle}{[(\langle
A^{2}\rangle-\langle A\rangle^{2})(\langle B^{2}\rangle-\langle
B\rangle^{2})]^{\frac{1}{2}}}~{},$ (31)

where $\langle\dots\rangle$ is a space-time average over the entire domain and
time after which the flow has become statistically stationary. The joint PDF
in Fig. 4 uses the so-called ‘z-scores’ of $\tau_{m}$ and
$\overline{\tau}_{\ell}$, denoted by $\tau^{*}_{m}$ and $\tau^{*}$ (all other
figures show the actual $\overline{\tau}_{\ell}$ and $\tau_{m}$ values). The
z-score of A is defined as

$A^{*}=\frac{A-\langle A\rangle}{\langle(A-\langle
A\rangle)^{2}\rangle^{1/2}},$ (32)

which helps focus on the spatio-temporal variation patterns after subtracting
the mean and normalizing by the standard deviation. For completeness, we plot
the joint PDFs of the actual $\tau_{m}$ and $\overline{\tau}_{\ell}$ in Fig.
14 in the Appendix. Using z-scores in the joint PDF only affects the slope of
the linear fit between $\tau_{m}$ and $\overline{\tau}_{\ell}$. The
correlation coefficient is the same whether using the actual values or their
z-scores,

$r(A^{*},B^{*})=r(A,B).$ (33)

Therefore, z-scores show how well the approximation $\tau_{m}$ is able to
capture the exact subscale physics $\overline{\tau}_{\ell}$ up to a
proportionality constant.

Fig. 5 plots the correlation coefficient as a function of $\ell$, which has
values exceeding $0.90$ for scales $\ell$ smaller than that of the peak and
decreasing to $\approx 0.85$ for scales $\ell=\pi$ (half the domain size).
These strong correlations justify our usage of expression (24c) as a proxy for
$\overline{\tau}_{\ell}({\bf u},\rho)$ to gain insight into its underlying
mechanics.

Figure 3: A 2D slice at $x=0$ from a snapshot of the 3D compressible
turbulence DNS, comparing the (left) exact
$\overline{\tau}_{\ell}(u_{x},\rho)$ at $\ell=0.19635$ with (right) its
approximation from eq. (24c),
$\tau_{m}=\frac{1}{3}M_{2}\,\ell^{2}\,\partial_{k}\overline{\rho}\,\partial_{k}\overline{u_{x}}$.
Figure 4: Comparison of $\overline{\tau}_{\ell}(\rho,{\bf u})$ with its
approximation
$\tau_{m}=\frac{1}{3}M_{2}\,\ell^{2}\,\partial_{k}\overline{\rho}\,\partial_{k}\overline{\bf
u}$ at $\ell=0.19635$. (A) x-component of $\tau$ and $\tau_{m}$ along a
diagonal transect through the domain from a single snapshot. (B) Joint PDF
between $\tau^{*}$ and $\tau_{m}^{*}$, where superscript ‘$*$’ indicates
z-scores, which emphasizes trends (see eq. (32) and associated discussion in
the text). Color bar indicates the number of grid-points on a logarithmic
scale, proportional to the probability. Black line is the best linear fit, and
the correlation-squared, $r^{2}$, quantifies the fraction of variability in
$\tau$ that is captured by $\tau_{m}$. (C), (D) Same as (A), (B) but for the
y-component. (E), (F) Same as (A), (B) but for the z-component. Figure 5:
Correlation coefficient $r$ between $\overline{\tau}_{\ell}(\rho,{\bf u})$ and
its approximation
$\tau_{m}=\frac{1}{3}M_{2}\,\ell^{2}\,\partial_{k}\overline{\rho}\,\partial_{k}\overline{\bf
u}$ as a function of scale $\ell$. (A) x-component. (B) y-component. (C)
z-component.

###  III.2 Geophysical fluid turbulence

We also ran simulations of the two-layer stacked shallow water equations
(34),(35) using the ocean general circulation model MOM6 adcroft2019gfdl ;
griffies2020primer . We simulated an ocean current hallberg2013using  that is
baroclinically unstable phillips1954energy ; vallisatmospheric , resulting in
mesoscale eddy generation and geostrophic turbulence, which we visualize in
Fig. 7.

Details of the simulation setup can be found in hallberg2013using . The
equations solved in MOM6 are

$\frac{\partial\boldsymbol{u}_{m}}{\partial
t}+(f+\omega)\,\hat{\bf{z}}\times\boldsymbol{u}_{m}=-{\mbox{\boldmath$\nabla$}}(M_{m}+|\boldsymbol{u}_{m}|^{2}/2)-{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}\boldsymbol{T}-\delta_{m2}\,C_{D}\,|\boldsymbol{u}_{2}|\,\boldsymbol{u}_{2}~{},$
(34) $\frac{\partial h_{m}}{\partial
t}+{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}(h_{m}{\bf
u}_{m})=(3-2m)\left[\gamma\left(\overline{\eta_{3/2}}^{x}-\eta_{3/2,ref}\right)-{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}\left(K_{h}{\mbox{\boldmath$\nabla$}}\eta_{3/2}\right)\right]~{}.$
(35)

The subscript $m=1,2$ labels the top and bottom layers, respectively, and
repeated indices are not summed in the thickness equation (35). $\eta_{1/2}$
is the vertical position of the top layer interface (surface height),
$\eta_{3/2}$ is the height of the interface between the two layers,
$\boldsymbol{u}_{m}$ is the horizontal velocity, $\omega$ is vertical
component of relative vorticity, $f=f_{0}+\beta y$ is the Coriolis parameter
for Earth’s rotation (i.e., the planetary vorticity), where $f_{0}$ is the
value of $f$ at southernmost boundary ($y=0$) and $\beta$ is the linear rate
of change of $f$ along meridional direction. $\boldsymbol{T}$ is the
horizontal stress tensor parameterized by Smagorinsky biharmonic viscosity
griffies2000biharmonic , $C_{D}$ is the dimensionless bottom drag coefficient
and $M_{1}=g\,\eta_{1/2}$ and
$M_{2}=M_{1}+\left(g\frac{\Delta\rho}{\rho_{0}}\right)\eta_{3/2}$ are the
pressure of the layers 1 and 2 normalized by $\rho_{0}$, the density of the
top layer. $\Delta\rho$ is the density difference between the two layers,
$h_{m}$ is the layer thickness, $K_{h}$ is the interface height diffusivity
coefficient, and $\gamma$ is a rate that is proportional to the mass flux
between the layers hallberg2013using . $K_{h}$ allows for the dissipation of
available potential energy gent1990isopycnal  and $\gamma$ forces the
interface height between the two layers to reference zonal mean state
$\eta_{3/2,ref}$. Here, $\overline{(\cdot)}^{x}$ represents a zonal average
(i.e., in the east-west direction). The flow is able to attain a statistically
steady state in the absence of direct wind forcing in eq. (34) by nudging the
layer interface $\eta_{3/2}$ to the reference profile $\eta_{3/2,ref}$ in eq.
(35), which drives the flow by providing a source of potential energy
hallberg2013using , mimicking the sloping of isopycnals by Ekman pumping in
the real ocean vallisatmospheric .

Figure 6: Initial condition for two-layer stacked shallow water simulation
following hallberg2013using . The top layer has zonal (east-west) velocity
with a Gaussian profile peaking at $\approx 0.35~{}$m/s (positive to the east)
while the bottom layer is initially quiescent. Bottom topography is flat. The
surface has initial height $\eta_{1/2}=0~{}$m. The interface height
$\eta_{3/2}$ between the two layers is initialized to the reference
$\eta_{3/2,ref}$ shown to satisfy geostrophic balance. $\eta_{3/2}$ is allowed
to evolve freely while being relaxed back to $\eta_{3/2,ref}$ at the rate
$\gamma$.

The flat bottom domain is $1200~{}$km $\times$ $1600~{}$km in horizontal
extent and is discretized on a Cartesian grid with a grid-cell size of 5 km.
The total depth is $2~{}$km and the bottom surface is flat. It is zonally
periodic and subject to free-slip walls at the northern and southern
boundaries. The initial conditions have a zonal jet $200~{}$km wide with a
Gaussian profile peaking at $\approx 0.35~{}$m/s for top layer and bottom
layer is quiescent, similar to initial conditions of hallberg2013using . The
interface between the two layers is such that it satisfies geostrophic balance
and is as shown in figure 6. The Coriolis parameter at the southernmost
boundary is $f_{0}=6.49\times 10^{-5}s^{-1}$ and $\beta=\partial_{y}f=2\times
10^{-11}m^{-1}s^{-1}$. The value of $\gamma=1.1574\times 10^{-6}~{}s^{-1}$ and
$K_{h}=8000~{}m^{2}s^{-1}$.

Figure 7: Visualization of the relative vorticity normalized by the Coriolis
paramater $f$ at the domain’s mid-latitude from our simulation of an eastward
flowing ocean current that is baroclinically unstable, resulting in eddy
generation and turbulence. The simulation solves the two-layer shallow water
equations, with the top and bottom panels showing the respective top and
bottom layers. Time increases from left to right: panels A1 and A2 are from
day 60, B1 and B2 from day 93, C1 and C3 from day 180. Our analysis uses model
data after day 169, when the flow has become statistically steady.

The simulation is run for 180 days, sufficiently long for the geostrophic
turbulence to develop as a result of baroclinic instabilities. For our
results, we analyze 100 snapshots taken every 2.5 hours starting from day 169.
Spectra of ${\bf u}_{1}$ and $h_{1}$ are shown in Appendix Fig. 21. For our
purpose, we only analyze the fields from the top layer since the flow in the
bottom layer is much weaker. Hereafter, we drop the layer subcript since all
results pertain to the top layer.

We compare $\overline{\tau}_{\ell}(h,{\bf u})$ with its approximation
$\frac{1}{2}\ell^{2}\,M_{2}\,\partial_{k}\overline{h}\,\partial_{k}\overline{{\bf
u}}$ in eq. (24c). We use a graded Top-Hat kernel following our previous work
analyzing oceanic flow aluie2018mapping ; Raietal2021 ,

$G_{\ell}(r)=0.5-0.5\tanh((r-\ell/2)/10),$ (36)

where $\ell$ denotes the kernel width. This kernel is essentially the same as
the Top-Hat in eq.(7), but with smoothed edges to avoid discretization noise
from the nonuniform grids aluie2018mapping  commonly used in general
circulation models. Since the model domain is flat, $r$ is a simple Euclidean
distance as in eq. (7), without complications due to Earth’s surface curvature
aluie2019convolutions .

Figure 8: Comparing the (left) exact $\overline{\tau}_{\ell}({\bf u},h)$ at
$\ell=100~{}$km with (right) its approximation from eq. (24c),
$\tau_{m}=\frac{1}{2}M_{2}\,\ell^{2}\,\partial_{k}\overline{h}\,\partial_{k}\overline{{\bf
u}}$, from a single snapshot. Top panels are for the zonal (eastward)
component. Bottom panels are for the meridional (northward) component. Grey
bands at northern and southern boundaries are excluded from our analysis to
avoid details of the coarse-graining next to boundaries.

Fig. 8 compares the $\overline{\tau}_{\ell}({\bf u},h)$ at $\ell=100~{}$km
with its approximation (24c). It indicates an excellent agreement. Fig. 9
(left column) shows a more quantitative evaluation of both components of the
subscale transport approximation, where we can see that the approximation is
remarkably accurate almost everywhere. Fig. 9 (right column) shows the joint
PDFs between $\overline{\tau}_{\ell}({\bf u},\rho)$ and its approximation for
both components (eastward and northward directions). They show correlations
exceeding $0.95$. In Figs. 18-20 in the Appendix, we show similar comparisons
at other length-scales of 50 km, 200 km, and 400 km. The joint PDF in Fig. 9
uses the z-scores of $\overline{\tau}_{\ell}$ and $\tau_{m}$, denoted by
$\tau^{*}$ and $\tau^{*}_{m}$ and (all other figures show the actual
$\overline{\tau}_{\ell}$ and $\tau_{m}$ values). As mentioned in eq. (33), the
correlation coefficient between $\tau^{*}$ and $\tau^{*}_{m}$ is the same as
that between $\overline{\tau}_{\ell}$ and $\tau_{m}$. We use z-scores to focus
on how well the approximation $\tau_{m}$ is able to capture the exact subscale
physics $\overline{\tau}_{\ell}$ up to a proportionality constant. For
completeness, we plot the joint PDFs of the actual $\tau_{m}$ and
$\overline{\tau}_{\ell}$ in Fig. 25 in the Appendix.

Noting that the spectral peak (Fig. 21) is at $\approx 300~{}$km, Fig. 10
plots the correlation coefficients between the exact
$\overline{\tau}_{\ell}({\bf u},h)$ and its approximation (24c) as a function
of $\ell$, which has values exceeding $0.90$ for scales $\ell$ smaller than
the that of the spectral peak and decreasing to $\approx 0.75$ for scales
$\ell=400~{}$km. These high correlation values justify our usage of expression
(24c) as a proxy for $\overline{\tau}_{\ell}({\bf u},h)$ to gain insight into
its underlying mechanics at least for scales $\ell$ smaller than that of the
spectral peak.

Figure 9:  Comparison of $\overline{\tau}_{\ell}(h,{\bf u})$ with its
approximation
$\tau_{m}=\frac{1}{2}M_{2}\,\ell^{2}\,\partial_{k}\overline{h}\,\partial_{k}\overline{\bf
u}$ at $\ell=100~{}$km. (A) Zonal (eastward) component of $\tau$ and
$\tau_{m}$ along the transect $y=802.5~{}$km from a snapshot on day 171 of the
simulation. (B) Joint PDF for all space-time points between $\tau^{*}$ and
$\tau_{m}^{*}$, where superscript ‘$*$’ indicates z-scores, which emphasizes
trends (see eq. (32) and associated discussion in the text). Color bar
indicates the number of grid-points on a logarithmic scale, proportional to
the probability. Red line is the best linear fit, and the correlation-squared,
$r^{2}$, quantifies the fraction of variability in $\tau$ that is captured by
$\tau_{m}$. Figure 10: Correlation coefficient $r$ between
$\overline{\tau}_{\ell}(h,{\bf u})$ and its approximation
$\tau_{m}=\frac{1}{2}M_{2}\,\ell^{2}\,\partial_{k}\overline{h}\,\partial_{k}\overline{\bf
u}$ as a function of scale $\ell$. (A) Zonal (eastward) direction. (B)
Meridional (northward) direction.

##  IV Transport by Subscale Vorticity and Strain

Insight from eq. (24c) into how the subscale physics contributes to large-
scale transport can be gleaned by using a fundamental result from Moffatt1983
, MiddletonLoder1989 , and griffies1998gent , who showed in different physical
systems the importance of scalar transport by an anti-symmetric skew diffusion
tensor, or equivalently an eddy-induced advective flux. We elaborate on this
idea in this section, which contains the main result of the paper in eq. (43).

Replacing $\overline{\tau}_{\ell}({\bf u},C)$ in eq. (10) by expression (24c),
we have

$\partial_{t}\overline{C}_{\ell}+{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}(\overline{{\bf
u}}_{\ell}\,\overline{C}_{\ell})={\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}\left(-A\,\ell^{2}\,\partial_{j}\overline{{\bf
u}}_{\ell}\,\partial_{j}\overline{C}_{\ell}\right),$ (37)

with the constant $A=M_{2}/n$. Therefore, transport of $\overline{C}_{\ell}$
by the subscale flux,

$R(C)={\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}[\underbrace{\left(-A~{}\ell^{2}\partial_{j}\overline{{\bf
u}}_{\ell}\right)}_{K_{ij}}\,\partial_{j}\overline{C}_{\ell}],$ (38)

resembles that of generalized transport tensor

$K_{ij}=-A\,\ell^{2}\,\partial_{j}\overline{u}_{i}.$ (39)

Note that $K_{ij}$ is a function of the length scale $\ell$ through both the
$\ell^{2}$ factor as well as the coarse-grained velocity $\overline{{\bf
u}}_{\ell}$. For notational brevity we drop the $\ell$ when writing individual
components of $\overline{{\bf u}}_{\ell}$.

Physical properties of the subscale transport tensor $K_{ij}$ can be
illuminated by separating it into a symmetric part $K^{\mbox{\tiny sym}}_{ij}$
representing subscale strain, and anti-symmetric part (also called skew-
symmetric) $K^{\mbox{\tiny skew}}_{ij}$ representing subscale vorticity. The
tensor is rewritten as

$K_{ij}=-A\,\ell^{2}\,\partial_{j}\overline{u}_{i}=-A~{}\ell^{2}\,[\underbrace{\frac{1}{2}\left(\partial_{j}\overline{u}_{i}+\partial_{i}\overline{u}_{j}\right)}_{\overline{S}_{ij}}+\underbrace{\frac{1}{2}\left(\partial_{j}\overline{u}_{i}-\partial_{i}\overline{u}_{j}\right)}_{\overline{\Omega}_{ij}}]=K^{\mbox{\tiny
sym}}_{ij}+K^{\mbox{\tiny skew}}_{ij}$ (40)

with vorticity and strain underlying kinematically distinct flow processes
(e.g., griffies1998gent, ). Note that the symmetric part satisfies
$K^{\mbox{\tiny sym}}_{ij}=K^{\mbox{\tiny sym}}_{ji}$ while the anti-symmetric
part satisfies $K^{\mbox{\tiny skew}}_{ij}=-K^{\mbox{\tiny skew}}_{ji}$.

Transport due to the subscale strain,

$\displaystyle R^{\mbox{\tiny sym}}(C)=\partial_{i}[K^{\mbox{\tiny
sym}}_{ij}\,\partial_{j}\overline{C}_{\ell}],$ (41)

can be thought of in the local eigenframe reference of $K^{\mbox{\tiny
sym}}_{ij}$ as anisotropic diffusion/anti-diffusion of $\overline{C}_{\ell}$
in directions that are not necessarily parallel to
${\mbox{\boldmath$\nabla$}}\overline{C}_{\ell}$. Such transport333The tensor
$K^{\mbox{\tiny sym}}$ would have to be positive semi-definite to qualify as a
diffusion tensor in the strict sense. was analyzed in detail in
higgins2004heat . Positive eigenvalues are associated with a stretching flow
direction (diffusion) while negative eigenvalues are associated with a
contracting flow direction (anti-diffusion). If the underlying flow is non-
solenoidal, ${\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}{\bf u}\neq
0$, such as in our simulations presented here, then $K^{\mbox{\tiny
sym}}_{ij}$ can also account for isotropic diffusion or coalescence due to the
subscale physics.

To understand the role of subscale vorticity, represented by the anti-
symmetric (or skew-symmetric) tensor $K^{\mbox{\tiny skew}}_{ij}$, we follow
the analysis of Moffatt1983 ; MiddletonLoder1989 ; griffies1998gent  and
rewrite

$\displaystyle R^{\mbox{\tiny skew}}(C)=\partial_{i}[K^{\mbox{\tiny
skew}}_{ij}\,\partial_{j}\overline{C}_{\ell}]=[\partial_{i}K^{\mbox{\tiny
skew}}_{ij}]\partial_{j}\overline{C}_{\ell}=\partial_{j}[\underbrace{\partial_{i}K^{\mbox{\tiny
skew}}_{ij}}_{-({\bf v}_{*})_{j}}\overline{C}_{\ell}]$ (42)

where the second equality follows from the chain rule and noting that the
product of $K^{\mbox{\tiny skew}}_{ij}$ (antisymmetric tensor) with
$\partial_{i}\partial_{j}\overline{C}_{\ell}$ (symmetric tensor) vanishes.
Similarly, the third equality follows from
$\partial_{i}\partial_{j}K^{\mbox{\tiny skew}}_{ij}=0$. The interpretation of
${\bf v}_{*}$ is an eddy-induced velocity arising from the stirring induced by
subscale vortical motions. It is straightforward to show that

${\bf v}_{*}=-{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}{\bf
K}^{\mbox{\tiny
skew}}=\frac{1}{2}\,A\,\ell^{2}\,{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\times$}}({\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\times$}}\overline{{\bf
u}}_{\ell}),$ (43)

which is the main result of this paper. An example of how it can arise in
swirling flow is sketched in Fig. 1. Note that this eddy-induced velocity is
solenoidal,

${\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}{\bf v}_{*}=0,$ (44)

due to the antisymmetry of $K^{\mbox{\tiny skew}}_{ij}$.

###  IV.1 Summary and Interpretation

The above result is a manifestation of the equivalence between transport by
(i) flux due to an anti-symmetric diffusion tensor and (ii) an advective flux
(Moffatt1983, ; MiddletonLoder1989, ; griffies1998gent, ). It suggests that
coarse-grained simulations such as LES may need to solve

$\displaystyle\partial_{t}\overline{C}_{\ell}+{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}\left[(\overline{{\bf
u}}_{\ell}+{\bf
v}_{*})\overline{C}_{\ell}\right]=-{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\cdot$}}\left({\bf
J}(C)\right)$ (45)

to represent the unresolved (subgrid) vorticity physics self-consistently.

${\bf J}(C)$ in eq. (45) can represent traditional subgrid models such as
turbulent (Smagorinsky) diffusion smagorinsky1963general , ${\bf
J}(C)=-\alpha_{turb}{\mbox{\boldmath$\nabla$}}\overline{C}_{\ell}$. The
evaluation of ${\bf
v}_{*}\sim{\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\times$}}({\mbox{\boldmath$\nabla$}}{\mbox{\boldmath$\times$}}\overline{{\bf
u}}_{\ell})$ presents no additional cost since the resolved velocity
$\overline{\bf u}_{\ell}$ is known. Figs. 17-26 in the Appendix show the joint
PDFs of ${\bf v}_{*}$ and $\overline{{\bf u}}_{\ell}$ at different length-
scales $\ell$, which demonstrate that ${\bf v}_{*}$ and the coarse velocity
$\overline{{\bf u}}_{\ell}$ are often of a comparable magnitude.

Our interpretation of ${\bf v}_{*}$ as an eddy-induced velocity representing
subscale vortical motions can be justified by appealing to a relation by
Johnson444Johnson’s relation (46) is exact only when filtering with a Gaussian
kernel. johnson2020energy ; johnson2021role ,

$\overline{\tau}_{\ell}(u_{i},C)=2\int_{0}^{\ell}d\delta\,\delta\,\overline{\left(\partial_{k}\overline{(u_{i})}_{\delta}\,\partial_{k}\overline{C}_{\delta}\right)}_{\sqrt{\ell^{2}-\delta^{2}}}~{}~{}.$
(46)

Eq. (46) highlights that subscale flux $\overline{\tau}_{\ell}({\bf u},C)$ is
due to the cumulative contribution of subscale velocity and scalar gradients
at all scales $<\ell$. Since subscale velocity gradients can be decomposed
into strain and vorticity, and since expression (24c) is the dominant
contribution, we are justified in interpreting ${\bf v}_{*}$ as an eddy-
induced velocity representing subscale vortical motions.

##  V Conclusion

We have shown that unlike subscale strain, which acts as an anisotropic
diffusion/anti-diffusion tensor, subscale vorticity’s contribution at leading
order is solely a conservative advection of coarse-grained scalars by an eddy-
induced velocity ${\bf v}_{*}$ proportional to the curl of vorticity. Our
analysis relied on the leading order term in Eyink’s multiscale gradient
expansion, which coincides with the nonlinear model and is known to be an
excellent _a priori_ representation of the subscale physics. Here, we provided
additional evidence for this excellent agreement from a three-dimensional
compressible turbulence simulation and from a two-layer shallow water model of
geophysical turbulence. Since the convergence of Eyink’s expansion and,
therefore, the dominance of the leading order term relies on ultraviolet
scale-locality, our results and conclusions may not hold at length-scales
larger than those of the spectral peak.

While the focus of this paper was on the transport of scalars, a similar
analysis may also apply to the transport of momentum. Indeed, this is
straightforward for incompressible flows where approximation (24c) has been
shown to work very well for momentum transportclark1979evaluation ; Liuetal94
; BorueOrszag98 . Johnson’s recent work on the energy cascade johnson2021role
seems to indicate that approximation (24c) may be deficient in its magnitude
(by a factor $\approx\frac{1}{2}$) compared to the true momentum flux. Yet,
since the two exhibit excellent spatio-temporal correlation as was shown in
Liuetal94 ; BorueOrszag98 , a simple rescaling of expression (24c) by a space-
time independent coefficient (possibly $\approx 2$) would suffice to match the
true momentum flux. Once approximation (24c) (possibly $\times 2$) is accepted
as representative of the subscale flux, the same argument we presented here
follows through for momentum after replacing $C$ with ${\bf u}$ in eq. (24c).
In flows with significant density variations such as compressible turbulence,
the subscale flux of momentum has a somewhat different functional form than in
eq. (10b) due to the density-weighted filtering required to satisfy the
inviscid criterionAluie13 ; livescu2020turbulence . Yet, we conjecture that a
result similar to what we derived here is still possible for subscale momentum
flux in variable-density flows, which was shown to be scale-local in Aluie11 .

## Acknowledgement

We thank Brandon Reichl at GFDL and the reviewers for valuable feedback and
suggestions. This research was funded by US DOE grant DE-SC0020229. Partial
support from US NSF grants PHY-2020249, OCE-2123496, US NASA grant
80NSSC18K0772, and US NNSA grant DE-NA0003856 is acknowledged. HA was also
supported by US DOE grants DE-SC0014318, DE-SC0019329, NSF grant PHY-2206380,
and US NNSA grant DE-NA0003914. JKS was also supported by US DOE grant DE-
SC0019329, US NNSA grant DE-NA0003914, and US NSF grant CBET-2143702.
Computing time was provided by NERSC under Contract No. DE-AC02-05CH11231 and
NASA’s HEC Program through NCCS at Goddard Space Flight Center. MOM6 is
publicly available and can be accessed at https://github.com/mom-ocean/MOM6.
The coarse-graining analysis was performed using codes similar to FlowSieve,
which is publicly available at https://github.com/husseinaluie/FlowSieve.
Other data used here may be obtained from the corresponding author upon
reasonable request.

## References

  * [1] Gregory L Eyink.  Multi-scale gradient expansion of the turbulent stress tensor.  Journal of Fluid Mechanics, 549(-1):159, February 2006.
  * [2] Pijush K Kundu, Ira M Cohen, and David R Dowling.  Fluid mechanics.  Academic press, 2015.
  * [3] S. B. Pope.  Turbulent flows.  Cambridge University Press, New York, 2000.
  * [4] C. Meneveau and J. Katz.  Scale-Invariance and Turbulence Models for Large-Eddy Simulation.  Ann. Rev. Fluid Mech., 32:1–32, 2000.
  * [5] G. L. Eyink.  Locality of turbulent cascades.  Physica D, 207:91–116, 2005.
  * [6] R S Strichartz.  A guide to distribution theory and Fourier transforms.  World Scientific Publishing Company, 2003.
  * [7] Lawrence C Evans.  Partial Differential Equations.  Amer Mathematical Society, April 2010.
  * [8] A. Leonard.  Energy Cascade in Large-Eddy Simulations of Turbulent Fluid Flows.  Adv. Geophys., 18:A237, 1974.
  * [9] M Germano.  Turbulence: the filtering approach.  Journal of Fluid Mechanics, 238:325–336, 1992.
  * [10] Ingrid Daubechies.  Ten lectures on wavelets.  SIAM, 1992.
  * [11] Hussein Aluie.  Coarse-grained incompressible magnetohydrodynamics: analyzing the turbulent cascades.  New Journal of Physics, January 2017.
  * [12] Gregory L Eyink.  Review of the onsager” ideal turbulence” theory.  arXiv preprint arXiv:1803.02223, 2018.
  * [13] Fernando F Grinstein, Len G Margolin, and William J Rider.  Implicit large eddy simulation, volume 10.  Cambridge university press Cambridge, 2007.
  * [14] Chad W Higgins, Marc B Parlange, and Charles Meneveau.  The heat flux and the temperature gradient in the lower atmosphere.  Geophysical research letters, 31(22), 2004.
  * [15] H. Aluie.  Scale decomposition in compressible turbulence.  Physica D: Nonlinear Phenomena, 247(1):54–65, March 2013.
  * [16] Dongxiao Zhao and Hussein Aluie.  Inviscid criterion for decomposing scales .  Physical Review Fluids, 3:054603, May 2018.
  * [17] Massimo Germano.  A proposal for a redefinition of the turbulent stresses in the filtered navier–stokes equations.  The Physics of fluids, 29(7):2323–2324, 1986.
  * [18] Andrey Nikolaevich Kolmogorov.  The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers.  Dokl. Akad. Nauk SSSR, 30(4):299–303, 1941.
  * [19] Alexander Mikhailovich Obukhov.  The structure of the temperature field in a turbulent flow.  Izv. Akad. Nauk. SSSR, Ser. Geogr. and Geophys., 13:58, 1949.
  * [20] Stanley Corrsin.  On the spectrum of isotropic temperature fluctuations in an isotropic turbulence.  Journal of Applied Physics, 22(4):469–473, 1951.
  * [21] Gregory L Eyink and Hussein Aluie.  Localness of energy cascade in hydrodynamic turbulence. i. smooth coarse graining.  Physics of Fluids, 21(11):115107, 2009.
  * [22] Hussein Aluie and Gregory L Eyink.  Localness of energy cascade in hydrodynamic turbulence. ii. sharp spectral filter.  Physics of Fluids, 21(11):115108, 2009.
  * [23] J Bardina, J H Ferziger, and W C Reynolds.  Improved subgrid-scale models for large-eddy simulation.  American Institute of Aeronautics and Astronautics, July 1980.
  * [24] H. Aluie.  Compressible Turbulence: The Cascade and its Locality.  Phys. Rev. Lett., 106(17):174502, April 2011.
  * [25] Michael Frazier, Michael W Frazier, Björn Jawerth, and Guido Weiss.  Littlewood-Paley theory and the study of function spaces.  Number 79 in Conference Board of the Mathematical Sciences regional conference series in mathematics. American Mathematical Society, Providence, Rhode Island, 1991.
  * [26] H. Aluie and G. Eyink.  Localness of energy cascade in hydrodynamic turbulence. II. Sharp spectral filter.  Phys. Fluids, 21(11):115108, November 2009.
  * [27] B. Vreman, B. Geurts, and H. Kuerten.  Realizability conditions for the turbulent stress tensor in large-eddy simulation.  J. Fluid Mech., 278:351–362, 1994.
  * [28] P. Constantin, E. Weinan, and E. S. Titi.  Onsager’s conjecture on the energy conservation for solutions of Euler’s equation.  Communications in Mathematical Physics, 165:207–209, October 1994\.
  * [29] Robert A Clark, Joel H Ferziger, and William Craig Reynolds.  Evaluation of subgrid-scale models using an accurately simulated turbulent flow.  Journal of fluid mechanics, 91(1):1–16, 1979.
  * [30] S Liu, C Meneveau, and J Katz.  On the properties of similarity subgrid-scale models as deduced from measurements in a turbulent jet.  Journal of Fluid Mechanics, 275:83–119, 1994.
  * [31] Fernando Porté-Agel, Marc B Parlange, Charles Meneveau, and William E Eichinger.  A priori field study of the subgrid-scale heat fluxes and dissipation in the atmospheric surface layer.  Journal of the atmospheric sciences, 58(18):2673–2698, 2001.
  * [32] Branko Kosović.  Subgrid-scale modelling for the large-eddy simulation of high-reynolds-number boundary layers.  Journal of Fluid Mechanics, 336:151–182, 1997.
  * [33] Freddy Bouchet.  Parameterization of two-dimensional turbulence using an anisotropic maximum entropy production principle.  arXiv preprint cond-mat/0305205, 2003.
  * [34] A Leonard.  Large-eddy simulation of chaotic convection and beyond.  In 35th Aerospace Sciences Meeting and Exhibit, page 204, 1997.
  * [35] V Borue and S A Orszag.  Local energy flux and subgrid-scale statistics in three-dimensional turbulence.  Journal of Fluid Mechanics, 366(1), 1998.
  * [36] Perry L Johnson.  Energy transfer from large to small scales in turbulence by multiscale nonlinear strain and vorticity interactions.  Physical review letters, 124(10):104501, 2020.
  * [37] Perry L Johnson.  On the role of vorticity stretching and strain self-amplification in the turbulence energy cascade.  Journal of Fluid Mechanics, 922, 2021.
  * [38] Hyung Suk Kang and Charles Meneveau.  Passive scalar anisotropy in a heated turbulent wake: new observations and implications for large-eddy simulations.  Journal of Fluid Mechanics, 442:161–170, 2001.
  * [39] Sergei G Chumakov.  A priori study of subgrid-scale flux of a passive scalar in isotropic homogeneous turbulence.  Physical Review E, 78(3):036313, 2008.
  * [40] Aarne Lees and Hussein Aluie.  Baropycnal work: a mechanism for energy transfer across scales.  Fluids, 4(2):92, 2019.
  * [41] Shriram Jagannathan and Diego A Donzis.  Reynolds and Mach number scaling in solenoidally-forced compressible turbulence using high-resolution direct numerical simulations.  Journal of Fluid Mechanics, 789:669–707, February 2016.
  * [42] Christoph Federrath, Ralf S. Klessen, and Wolfram Schmidt.  The density probability distribution in compressible isothermal turbulence: Solenoidal versus compressive forcing.  The Astrophysical Journal Letters, 688(2):L79, 2008.
  * [43] V Eswaran and SB Pope.  An examination of forcing in direct numerical simulations of turbulence.  Computers & Fluids, 16(3):257–278, 1988.
  * [44] Alistair Adcroft, Whit Anderson, V Balaji, Chris Blanton, Mitchell Bushuk, Carolina O Dufour, John P Dunne, Stephen M Griffies, Robert Hallberg, Matthew J Harrison, et al.  The gfdl global ocean and sea ice model om4. 0: Model description and simulation features.  Journal of Advances in Modeling Earth Systems, 11(10):3167–3211, 2019.
  * [45] Stephen M Griffies, Alistair Adcroft, and Robert W Hallberg.  A primer on the vertical lagrangian-remap method in ocean models based on finite volume generalized vertical coordinates.  Journal of Advances in Modeling Earth Systems, 12(10):e2019MS001954, 2020.
  * [46] Robert Hallberg.  Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects.  Ocean Modelling, 72:92–103, 2013.
  * [47] Norman A Phillips.  Energy transformations and meridional circulations associated with simple baroclinic waves in a two-level, quasi-geostrophic model.  Tellus, 6(3):274–286, 1954.
  * [48] GK Vallis.  Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation, 2006.
  * [49] Stephen M Griffies and Robert W Hallberg.  Biharmonic friction with a smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models.  Monthly Weather Review, 128(8):2935–2946, 2000.
  * [50] Peter R Gent and James C Mcwilliams.  Isopycnal mixing in ocean circulation models.  Journal of Physical Oceanography, 20(1):150–155, 1990.
  * [51] Hussein Aluie, Matthew Hecht, and Geoffrey K Vallis.  Mapping the energy cascade in the north atlantic ocean: The coarse-graining approach.  Journal of Physical Oceanography, 48(2):225–244, 2018.
  * [52] Shikhar Rai, Matthew Hecht, Matthew Maltrud, and Hussein Aluie.  Scale of oceanic eddy killing by wind from global satellite observations.  Science Advances, 7(28):eabf4920, 2021.
  * [53] Hussein Aluie.  Convolutions on the sphere: commutation with differential operators.  GEM-International Journal on Geomathematics, 10(1):9, 2019.
  * [54] H.K. Moffatt.  Transport effects associated with turbulence with particular attention to the influence of helicity.  Reports on Progress in Physics, 46:621–664, 1983.
  * [55] J. F. Middleton and J. W. Loder.  Skew fluxes in polarized wave fields.  Journal of Physical Oceanography, 19:68–76, 1989.
  * [56] Stephen M Griffies.  The gent–mcwilliams skew flux.  Journal of Physical Oceanography, 28(5):831–841, 1998.
  * [57] Joseph Smagorinsky.  General circulation experiments with the primitive equations: I. the basic experiment.  Monthly weather review, 91(3):99–164, 1963.
  * [58] Daniel Livescu.  Turbulence with large thermal and compositional density variations.  Annual Review of Fluid Mechanics, 52:309–341, 2020.
  * [59] Mahmoud Sadek and Hussein Aluie.  Extracting the spectrum of a flow by spatial filtering.  Physical Review Fluids, 3(12):124610, December 2018.
  * [60] Michele Buzzicotti, Benjamin A Storer, Stephen M Griffies, and Hussein Aluie.  A coarse-grained decomposition of surface geostrophic kinetic energy in the global ocean.  Earth and Space Science Open Archive, page 58, 2021.

## Appendix

###  V.1 Additional supporting material from the compressible turbulence
simulation

Figure 11:  Same as in Fig. 4 for $\ell=0.3927$. Figure 12:  Same as Fig. 4
for $\ell=0.7854$. Figure 13:  Same as Fig. 4 for $\ell=1.5708$. Figure 14:
Similar to the joint PDFs of the z-scores in Fig. 4, here we show the joint
PDFs of the actual $\overline{\tau}_{\ell}(\rho,{\bf u})$ and its
approximation
$\tau_{m}=\frac{1}{3}M_{2}\ell^{2}\partial_{k}\overline{\rho}\,\partial_{k}\overline{\bf
u}$. Top row (A, B, C) shows the joint PDFs between $\tau$ and $\tau_{m}$ at
at $\ell=0.19635$. Second row (D, E, F) is for $\ell=0.3927$. Third row (D, H,
I) is for $\ell=0.7854$. Bottom row (J, K, L) is for $\ell=1.5708$. Left,
middle, and right columns show, the x, y and z components, respectively.

Figure 15:  Spectra of (left) velocity, $E^{u}(k)$, and (right) density,
$E^{\rho}(k)$. Both decay with wavenumber $k$ sufficiently faster than the
$k^{-1}$ required for the multiscale gradient expansion to converge. Here
$E^{u}(k)=\sum_{\big{|}|{\bf k}|-k\big{|}<0.5}\frac{1}{2}|\hat{\bf u}({\bf
k})|^{2}$ and $E^{\rho}(k)=\sum_{\big{|}|{\bf
k}|-k\big{|}<0.5}|\hat{\rho}({\bf k})|^{2}$, where $\hat{\bf u}({\bf k})$ and
$\hat{\rho}({\bf k})$ Fourier coefficients.

Figure 16: Same as in Fig. 3, comparing the (left column) exact
$\overline{\tau}_{\ell}(\rho,{\bf u})$ with (right column) its approximation
for the (top row) y-components and (bottom row) z-components. Figure 17: Joint
PDFs of ${\bf v}_{*}$ and $\overline{\bf u}$. Color bar indicates the number
of grid-points on a logarithmic scale, proportional to the probability. Slope
of best linear fit and correlation-squared, $r^{2}$, are shown in the top left
of each panel. The best linear fits are omitted in the the plots due to small
slopes. The rms values are shown in the bottom right of each panel. Top row
(A, B, C) shows the joint PDFs between ${\bf v}_{*}$ and $\overline{\bf u}$ at
at $\ell=0.19635$. Second row (D, E, F) is for $\ell=0.3927$. Third row (D, H,
I) is for $\ell=0.7854$. Bottom row (J, K, L) is for $\ell=1.5708$. Left,
middle, and right columns show, the x, y and z components, respectively. This
figure highlights that the magnitudes of ${\bf v}_{*}$ and the coarse velocity
$\overline{{\bf u}}_{\ell}$ are often comparable.

###  V.2 Additional supporting material from the shallow water simulation

Figure 18:  Same as in Figure 9 for $\ell=50~{}km$. Figure 19:  Same as in
Figure 9 for $\ell=200~{}km$. Figure 20:  Same as in Figure 9 for
$\ell=400~{}km$. Figure 21: Spectra of (left) ${\bf u}$ and (right) $h$
computed from the two-layer stacked shallow water simulation. These are
filtering spectra, calculated by filtering at successive length-scales as
described in [59, 60]. Figure 22: Same as in Fig. 8 for $\ell=50~{}$km. Figure
23:  Same as in Figure 8 for $\ell=200~{}km$. Figure 24:  Same as in Figure 8
for $\ell=400~{}km$. Figure 25: Similar to the joint PDFs of the z-scores in
Fig. 9, here we show the joint PDFs of the actual
$\overline{\tau}_{\ell}(h,{\bf u})$ and its approximation
$\tau_{m}=\frac{1}{2}M_{2}\ell^{2}\partial_{k}\overline{h}\,\partial_{k}\overline{\bf
u}$, without normalization. Top and bottom panels show zonal (eastward) and
meridional (northward) components, respectively. A1 and A2 are for
$\ell=50~{}$km, B1 and B2 for $100~{}$km, C1 and C2 for $200~{}$km, D1 and D2
for $400~{}$km. Figure 26: Joint Probability Density of ${\bf v}_{*}$ and
$\overline{{\bf u}}$. The colorbars are in log scale. The top panels show
zonal(eastward) componet and the bottom panels show meridional(northward)
component. A1 and A2 are for $\ell$ = 50 km, B1 and B2 are for $\ell$ = 100
km, C1 and C2 are for $\ell$ = 200 km, D1 and D2 are for $\ell$ = 400 km. The
red lines are linear regression fits of which the slope and goodness of
fit($r^{2}$) are in the text in each panel. The rms values are shown in the
bottom right of each panel. This figure highlights that ${\bf v}_{*}$ is often
comparable in magnitude relative to the coarse velocity $\overline{{\bf
u}}_{\ell}$.