# Surface Quasigeostrophic Turbulence in Variable Stratification

###### Abstract

Numerical and observational evidence indicates that, in regions where mixed-
layer instability is active, the surface geostrophic velocity is largely
induced by surface buoyancy anomalies. Yet, in these regions, the observed
surface kinetic energy spectrum is steeper than predicted by uniformly
stratified surface quasigeostrophic theory. By generalizing surface
quasigeostrophic theory to account for variable stratification, we show that
surface buoyancy anomalies can generate a variety of dynamical regimes
depending on the stratification’s vertical structure. Buoyancy anomalies
generate longer range velocity fields over decreasing stratification and
shorter range velocity fields over increasing stratification. As a result, the
surface kinetic energy spectrum is steeper over decreasing stratification than
over increasing stratification. An exception occurs if the near surface
stratification is much larger than the deep ocean stratification. In this
case, we find an extremely local turbulent regime with surface buoyancy
homogenization and a steep surface kinetic energy spectrum, similar to
equivalent barotropic turbulence. By applying the variable stratification
theory to the wintertime North Atlantic, and assuming that mixed-layer
instability acts as a narrowband small-scale surface buoyancy forcing, we
obtain a predicted surface kinetic energy spectrum between $k^{-4/3}$ and
$k^{-7/3}$, which is consistent with the observed wintertime $k^{-2}$
spectrum. We conclude by suggesting a method of measuring the buoyancy
frequency’s vertical structure using satellite observations.

##  1 Introduction

###  1.1 Geostrophic flow induced by surface buoyancy

Geostrophic flow in the upper ocean is induced by either interior potential
vorticity anomalies, $q$, or surface buoyancy anomalies, $b|_{z=0}$. At first,
it was assumed that the surface geostrophic flow observed by satellite
altimeters is due to interior potential vorticity (Stammer, 1997; Wunsch,
1997). It was later realized, however, that under certain conditions, upper
ocean geostrophic flow can be inferred using the surface buoyancy anomaly
alone (Lapeyre and Klein, 2006; LaCasce and Mahadevan, 2006). Subsequently,
Lapeyre (2009) used a numerical ocean model to show that the surface buoyancy
induced geostrophic flow dominates the $q$-induced geostrophic flow over a
large fraction of the North Atlantic in January. Lapeyre then concluded that
the geostrophic velocity inferred from satellite altimeters in the North
Atlantic must usually be due to surface buoyancy anomalies instead of interior
potential vorticity.

Similar conclusions have been reached in later numerical studies using the
effective surface quasigeostrophic (eSQG, Lapeyre and Klein, 2006) method. The
eSQG method aims to reconstruct three-dimensional velocity fields in the upper
ocean: it assumes that surface buoyancy anomalies generate an exponentially
decaying streamfunction with a vertical attenuation determined by the buoyancy
frequency, as in the uniformly stratified surface quasigeostrophic model (Held
et al., 1995). Because the upper ocean does not typically have uniform
stratification, an ”effective” buoyancy frequency is used, which is also
intended to account for interior potential vorticity anomalies (Lapeyre and
Klein, 2006). In practice, however, this effective buoyancy frequency is
chosen to be the vertical average of the buoyancy frequency in the upper
ocean. Qiu et al. (2016) derived the surface streamfunction from sea surface
height in a $\nicefrac{{1}}{{30}}^{\circ}$ model of the Kuroshio Extension
region in the North Pacific and used the eSQG method to reconstruct the three-
dimensional vorticity field. They found correlations of 0.7-0.9 in the upper
1000 m between the reconstructed and model vorticity throughout the year. This
result was also found to hold in a $\nicefrac{{1}}{{48}}^{\circ}$ model with
tidal forcing (Qiu et al., 2020).

A clearer test of whether the surface flow is induced by surface buoyancy is
to reconstruct the geostrophic flow directly using the sea surface buoyancy or
temperature (Isern‐Fontanet et al., 2006). This approach was taken by
Isern‐Fontanet et al. (2008) in the context of a
$\nicefrac{{1}}{{10}}^{\circ}$ numerical simulation of the North Atlantic.
When the geostrophic velocity is reconstructed using sea surface temperature,
correlations between the reconstructed velocity and the model velocity
exceeded 0.7 over most of the North Atlantic in January. Subsequently,
Miracca-Lage et al. (2022) used a reanalysis product with a grid spacing of 10
km to reconstruct the geostrophic velocity using both sea surface buoyancy and
temperature over certain regions in the South Atlantic. The correlations
between the reconstructed streamfunctions and the model streamfunction had a
seasonal dependence, with correlations of 0.7-0.8 in winter and 0.2-0.4 in
summer.

Observations also support the conclusion that a significant portion of the
surface geostrophic flow may be due to surface buoyancy anomalies over a
substantial fraction of the World Ocean. González-Haro and Isern-Fontanet
(2014) reconstructed the surface streamfunction using
$\nicefrac{{1}}{{3}}^{\circ}$ satellite altimeter data (for sea surface
height) and $\nicefrac{{1}}{{4}}^{\circ}$ microwave radiometer data (for sea
surface temperature). If the surface geostrophic velocity is due to sea
surface temperature alone, then the streamfunction constructed from sea
surface temperature should be perfectly correlated with the streamfunction
constructed from sea surface height. The spatial correlations between the two
streamfunctions was found to be seasonal. For the wintertime northern
hemisphere, high correlations (exceeding 0.7-0.8) are observed near the Gulf
Stream and Kuroshio whereas lower correlations (0.5-0.6) are seen in the
eastern recirculating branch of North Atlantic and North Pacific gyres [a
similar pattern was found by Isern‐Fontanet et al. (2008) and Lapeyre (2009)].
In summer, correlations over the North Atlantic and North Pacific plummet to
0.2-0.5, again with lower correlations in the recirculating branch of the
gyres. In contrast to the strong seasonality observed in the northern
hemisphere, correlation over the Southern Ocean typically remain larger than
0.8 throughout the year.

Another finding is that more of the surface geostrophic flow is due to surface
buoyancy anomalies in regions with high eddy kinetic energy, strong thermal
gradients, and deep mixed layers (Isern‐Fontanet et al., 2008; González-Haro
and Isern-Fontanet, 2014; Miracca-Lage et al., 2022). These are the same
conditions under which we expect mixed-layer baroclinic instability to be
active (Boccaletti et al., 2007; Mensa et al., 2013; Sasaki et al., 2014;
Callies et al., 2015). Indeed, one model of mixed-layer instability consists
of surface buoyancy anomalies interacting with interior potential vorticity
anomalies at the base of the mixed-layer (Callies et al., 2016). The dominance
of the surface buoyancy induced velocity in regions of mixed-layer instability
suggests that, to a first approximation, we can think of mixed-layer
instability as energizing the surface buoyancy induced part of the flow
through vertical buoyancy fluxes and the concomitant release of kinetic energy
at smaller scales.

###  1.2 Surface quasigeostrophy in uniform stratification

The dominance of the surface buoyancy induced velocity suggests that a useful
model for upper ocean geostrophic dynamics is the surface quasigeostrophic
model (Held et al., 1995; Lapeyre, 2017), which describes the dynamics induced
by surface buoyancy anomalies over uniform stratification. The primary
difference between surface quasigeostrophic dynamics and two-dimensional
barotropic dynamics (Kraichnan, 1967) is that surface quasigeostrophic eddies
have a shorter interaction range than their two-dimensional barotropic
counterparts. One consequence of this shorter interaction range is a flatter
kinetic energy spectrum (Pierrehumbert et al., 1994). Letting $k$ be the
horizontal wavenumber, then two-dimensional barotropic turbulence theory
predicts a kinetic energy spectrum of $k^{-5/3}$ upscale of small-scale
forcing and a $k^{-3}$ spectrum downscale of large-scale forcing (Kraichnan,
1967). If both types of forcing are present, then we expect a spectrum between
$k^{-5/3}$ and $k^{-3}$, with the realized spectrum depending on the relative
magnitude of small-scale to large-scale forcing (Lilly, 1989; Maltrud and
Vallis, 1991). In contrast, the corresponding spectra for surface
quasigeostrophic turbulence are $k^{-1}$ (upscale of small-scale forcing) and
$k^{-5/3}$ (downscale of large-scale forcing) (Blumen, 1978), both of which
are flatter than the corresponding two-dimensional barotropic spectra.111The
uniformly stratified geostrophic turbulence theory of Charney (1971) provides
spectral predictions similar to the two-dimensional barotropic theory (See
Callies and Ferrari, 2013).

The above considerations lead to the first discrepancy between the surface
quasigeostrophic model and ocean observations. As we have seen, we expect
wintertime surface geostrophic velocities near major extratropical currents to
be primarily due to surface buoyancy anomalies. Therefore, the predictions of
surface quasigeostrophic theory should hold. If we assume that mesoscale
baroclinic instability acts as a large-scale forcing and that mixed-layer
baroclinic instability acts as a small-scale forcing to the upper ocean (we
assume a narrowband forcing in both cases, although this may not be the case,
see Khatri et al., 2021), then we expect a surface kinetic energy spectrum
between $k^{-1}$ and $k^{-5/3}$. However, both observations and numerical
simulations of the Gulf Stream and Kuroshio find kinetic energy spectra close
to $k^{-2}$ in winter (Sasaki et al., 2014; Callies et al., 2015; Vergara et
al., 2019), which is steeper than predicted.

The second discrepancy relates to the surface transfer function implied by
uniformly stratified surface quasigeostrophic theory. The surface transfer
function, $\mathcal{F}(\bm{k})$, is defined as (Isern-Fontanet et al., 2014)

$\mathcal{F}(\bm{k})=\frac{\hat{\psi}_{\bm{k}}}{\hat{b}_{\bm{k}}},$ (1)

where $\hat{\psi}_{\bm{k}}$ and $\hat{b}_{\bm{k}}$ are the Fourier amplitudes
of the geostrophic streamfunction, $\psi$, and the buoyancy, $b$, at the
ocean’s surface, and $\bm{k}$ is the horizontal wavevector. Uniformly
stratified surface quasigeostrophic theory predicts an isotropic transfer
function $\mathcal{F}(k)\sim k^{-1}$ (Held et al., 1995). Using a
$\nicefrac{{1}}{{12}}^{\circ}$ ocean model and focusing on the western coast
of Australia, González-Haro et al. (2020) confirmed that the transfer function
between sea surface temperature and sea surface height is indeed isotropic but
found that the transfer function is generally steeper than $k^{-1}$. In
another study using a $\nicefrac{{1}}{{16}}^{\circ}$ model of the
Mediterranean Sea, Isern-Fontanet et al. (2014) found that the transfer
function below 100 km has seasonal dependence closely related to mixed-layer
depth, fluctuating between $k^{-1}$ and $k^{-2}$.

In the remainder of this article, we account for these discrepancies by
generalizing the uniformly stratified surface quasigeostrophic model (Held et
al., 1995) to account for variable stratification (section 2). Generally, we
find that the surface kinetic energy spectrum in surface quasigeostrophic
turbulence depends on the stratification’s vertical structure (section 3); we
recover the Blumen (1978) spectral predictions only in the limit of uniform
stratification. Stratification controls the kinetic energy spectrum by
modifying the interaction range of surface quasigeostrophic eddies, and we
illustrate this dependence by examining the turbulence under various idealized
stratification profiles (section 4). We then apply the theory to the North
Atlantic in both winter and summer, and find that the surface transfer
function is seasonal, with a $\mathcal{F}(k)\sim k^{-3/2}$ dependence in
winter and a $\mathcal{F}(k)\sim k^{-1/2}$ dependence in summer. Moreover, in
the wintertime North Atlantic, the theory predicts a surface kinetic energy
spectrum between $k^{-4/3}$ and $k^{-7/3}$, which is consistent with both
observations and numerical simulations (section 5). Finally, in section 6, we
discuss the validity of theory at other times and locations.

##  2 The inversion function

###  2.1 Physical space equations

Consider an ocean of depth $H$ with zero interior potential vorticity $(q=0)$
so that the geostrophic streamfunction satisfies

$\nabla^{2}\psi+\frac{\partial}{\partial
z}\left(\frac{1}{\sigma^{2}}\frac{\partial\psi}{\partial
z}\right)=0\quad\text{for }z\in(-H,0).$ (2)

In this equation, $\nabla^{2}$ is the horizontal Laplacian, $\psi$ is the
geostrophic streamfunction, and

$\sigma(z)=N(z)/f,$ (3)

where $N(z)$ is the depth-dependent buoyancy frequency and $f$ is the constant
local value of the Coriolis frequency. We refer to $\sigma(z)$ as the
_stratification_ for the remainder of this article. The horizontal geostrophic
velocity is obtained from $\bm{u}=\hat{\bm{z}}\times\bm{\nabla}\bm{\psi}$
where $\hat{\bm{z}}$ is the vertical unit vector.

The upper surface potential vorticity is given by (Bretherton, 1966)

$\theta=-\frac{1}{\sigma_{0}^{2}}\frac{\partial\psi}{\partial
z}\bigg{|}_{z=0},$ (4)

where $\sigma_{0}=\sigma(0)$. The surface potential vorticity is related to
the surface buoyancy anomaly through

$b|_{z=0}=-f\,\sigma_{0}^{2}\,\theta.$ (5)

The time-evolution equation at the upper boundary is given by

$\frac{\partial\theta}{\partial
t}+\mathrm{J}\left(\psi,\theta\right)=F-D\quad\text{at }z=0,$ (6)

where
$\mathrm{J}\left(\theta,\psi\right)=\partial_{x}\theta\,\partial_{y}\psi-\partial_{y}\theta\,\partial_{x}\psi$
represents the advection of $\theta$ by the horizontal geostrophic velocity
$\bm{u}$, $F$ is the buoyancy forcing at the upper boundary, and $D$ is the
dissipation.

We assume a bottom boundary condition of

$\psi\rightarrow 0\text{ as }z\rightarrow-\infty,$ (7)

which is equivalent to assuming the bottom boundary, $z=-H$, is irrelevant to
the dynamics. In section 5, we find that this assumption is valid in the mid-
latitude North Atlantic open ocean at horizontal scales smaller than $\approx
500$ km. We consider alternative boundary conditions in appendix A.

###  2.2 Fourier space equations

Assuming a doubly periodic domain in the horizontal prompts us to consider the
Fourier expansion of $\psi$,

$\psi(\bm{r},z,t)=\sum_{\bm{k}}\hat{\psi}_{\bm{k}}(t)\,\Psi_{k}(z)\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{r}},$
(8)

where $\bm{r}=(x,y)$ is the horizontal position vector, $\bm{k}=(k_{x},k_{y})$
is the horizontal wavevector, and $k=\left\lvert\bm{k}\right\rvert$ is the
horizontal wavenumber. The wavenumber dependent non-dimensional vertical
structure, $\Psi_{k}(z)$, is determined by the boundary-value problem222To
derive the vertical structure equation (9), substitute the Fourier
representation (8) into the vanishing potential vorticity condition (2),
multiply by $\mathrm{e}^{-i\bm{l}\cdot\bm{r}}$, take an area average, and use
the identity
$\frac{1}{A}\,\int_{A}\mathrm{e}^{\mathrm{i}\left(\bm{k}-\bm{l}\right)\cdot\bm{r}}\,\mathrm{d}{\bm{r}}=\delta_{\bm{k},\bm{l}}$
where $\delta_{\bm{k},\bm{l}}$ is the Kronecker delta, and $A$ is the
horizontal area of the periodic domain.

$-\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{1}{\sigma^{2}}\frac{\mathrm{d}\Psi_{k}}{\mathrm{d}z}\right)+k^{2}\,\Psi_{k}=0,$
(9)

with upper boundary condition

$\Psi_{k}(0)=1$ (10)

and bottom boundary condition

$\Psi_{k}\rightarrow 0\text{ as }z\rightarrow-\infty.$ (11)

The upper boundary condition (10) is a normalization for the vertical
structure, $\Psi_{k}(z)$, chosen so that

$\psi(\bm{r},z=0,t)=\sum_{\bm{k}}\hat{\psi}_{\bm{k}}(t)\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{r}}.$
(12)

Consequently, the surface potential vorticity (4) is given by

$\theta(\bm{r},t)=\sum_{\bm{k}}\hat{\theta}_{\bm{k}}(t)\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{r}},$
(13)

where

$\hat{\theta}_{\bm{k}}=-m(k)\,\hat{\psi}_{\bm{k}},$ (14)

and the inversion function $m(k)$ (with dimensions of inverse length) is
defined as

$m(k)=\frac{1}{\sigma_{0}^{2}}\,\frac{\mathrm{d}\Psi_{k}(0)}{\mathrm{d}z}.$
(15)

In all our applications, we find the inversion function to be a positive
monotonically increasing function of $k$ [i.e., $m(k)>0$ and
$\mathrm{d}m/\mathrm{d}k\geq 0$]. The inversion function is related to the
transfer function (1) through

$\mathcal{F}(k)=\frac{1}{f\,\sigma_{0}^{2}\,m(k)}=\left[f\,\frac{\mathrm{d}\Psi_{k}(0)}{\mathrm{d}z}\right]^{-1},$
(16)

which shows that the transfer function, evaluated at a wavenumber $k$, is
related to the characteristic vertical scale of $\Psi_{k}(z)$.

###  2.3 The case of constant stratification

To recover the well-known case of the uniformly stratified surface
quasigeostrophic model (Held et al., 1995), set $\sigma=\sigma_{0}$. Then
solving the vertical structure equation (9) along with boundary conditions
(10) and (11) yields the exponentially decaying vertical structure,

$\Psi_{k}(z)=\mathrm{e}^{\sigma_{0}\,k\,z}.$ (17)

Substituting $\Psi_{k}(z)$ into the definition of the inversion function (15),
we obtain

$m(k)=k/\sigma_{0},$ (18)

and hence [through the inversion relation (13)] a linear-in-wavenumber
inversion relation of

$\hat{\theta}_{\bm{k}}=-(k/\sigma_{0})\,\hat{\psi}_{\bm{k}}.$ (19)

In appendix A, we show that $m(k)\rightarrow k/\sigma_{0}$ as
$k\rightarrow\infty$ for arbitrary stratification $\sigma(z)$. Therefore, at
sufficiently small horizontal scales, surface quasigeostrophic dynamics
behaves as in constant stratification regardless of the functional form of
$\sigma(z)$.

##  3 Surface quasigeostrophic turbulence

Suppose a two-dimensional barotropic fluid is forced in the wavenumber
interval $[k_{1},k_{2}]$. In such a fluid, Kraichnan (1967) argued that two
inertial ranges will form: one inertial range for $k<k_{1}$ where kinetic
energy cascades to larger scales and another inertial range for $k>k_{2}$
where enstrophy cascades to smaller scales. Kraichnan’s argument depends on
three properties of two-dimensional vorticity dynamics. First, that there are
two independent conserved quantities; namely, the kinetic energy and the
enstrophy. Second, that turbulence is sufficiently local in wavenumber space
so that the only available length scale is $k^{-1}$. Third, that the inversion
relation between the vorticity and the streamfunction is scale invariant.

There are two independent conserved quantities in surface quasigeostrophic
dynamics, as in Kraichnan’s two-dimensional fluid; namely the total energy,
$E$, and the surface potential enstrophy, $P$. However, the second and third
properties of two-dimensional vorticity dynamics do not hold for surface
quasigeostrophic dynamics. Even if the turbulence is local in wavenumber
space, there are two available length scales at each wavenumber $k$; namely,
$k^{-1}$ and $[m(k)]^{-1}$. Moreover, the inversion relation (14) is generally
not scale invariant.333A function $m(k)$ is scale invariant if $m(\lambda
k)=\lambda^{s}m(k)$ for all $\lambda$, where $s$ is a real number. In
particular, note that power laws, $m(k)=k^{\alpha}$, are scale invariant.
Therefore, the arguments in Kraichnan (1967) do not hold in general for
surface quasigeostrophic dynamics.

Even so, in the remainder of this section we show that there is a net inverse
cascade of total energy and a net forward cascade of surface potential
enstrophy even if there are no inertial ranges in the turbulence. Then we
consider the circumstances under which we expect inertial ranges to form.
Finally, assuming the existence of an inertial range, we derive the spectra
for the cascading quantities. We begin, however, with some definitions.

###  3.1 Quadratic quantities

The two quadratic quantities needed for the cascade argument are the volume-
integrated total mechanical energy per mass per unit area,

$\displaystyle E$
$\displaystyle=\frac{1}{2\,A}\int_{V}\left(\left\lvert\bm{\nabla}\psi\right\rvert^{2}+\frac{1}{\sigma^{2}}\left\lvert\frac{\partial\psi}{\partial
z}\right\rvert^{2}\right)\mathrm{d}V$ (20)
$\displaystyle=-\frac{1}{2}\,\overline{\psi|_{z=0}\,\theta}=\frac{1}{2}\sum_{\bm{k}}m(k)\,\left\lvert\hat{\psi}_{\bm{k}}\right\rvert^{2},$

and the upper surface potential enstrophy,

$P=\frac{1}{2}\,\overline{\theta^{2}}=\frac{1}{2}\sum_{\bm{k}}\left[m(k)\right]^{2}\left\lvert\hat{\psi}_{\bm{k}}\right\rvert^{2},$
(21)

where the overline denotes an area average over the periodic domain. Both
quantities are time-independent in the absence of forcing and dissipation, as
can be seen by multiplying the time-evolution equation (6) by either
$-\psi|_{z=0}$ or $\theta$ and taking an area average.

Two other quadratics we use are the surface kinetic energy

$K=\frac{1}{2}\overline{\left\lvert\bm{\nabla}\psi\right\rvert_{z=0}^{2}}=\frac{1}{2}\sum_{\bm{k}}k^{2}\,\left\lvert\hat{\psi}_{\bm{k}}\right\rvert^{2}$
(22)

and the surface streamfunction variance

$S=\frac{1}{2}\overline{\left(\psi|_{z=0}\right)^{2}}=\frac{1}{2}\sum_{\bm{k}}\left\lvert\hat{\psi}_{\bm{k}}\right\rvert^{2}.$
(23)

Moreover, the isotropic spectrum $\mathscr{A}(k)$ of a quadratic quantity $A$
is defined by

$A=\int_{0}^{\infty}\mathscr{A}(k)\,\mathrm{d}k,$ (24)

so that the isotropic spectra of $E,P,K,$ and $S$ are given by
$\mathscr{E}(k),\mathscr{P}(k),\mathscr{K}(k),$ and $\mathscr{S}(k)$. The
isotropic spectra are then related by

$\mathscr{P}(k)=m(k)\,\mathscr{E}(k)=\left[m(k)\right]^{2}\,\mathscr{S}(k)$
(25)

and

$\mathscr{K}(k)=k^{2}\,\mathscr{S}(k).$ (26)

For $\mathscr{A}(k)=\mathscr{E}(k)$ or $\mathscr{A}(k)=\mathscr{P}(k)$, there
is a time-evolution equation of the form (Gkioulekas and Tung, 2007)

$\frac{\partial\mathscr{A}(k)}{\partial t}+\frac{\partial\Pi_{A}(k)}{\partial
k}=F_{A}(k)-D_{A}(k),$ (27)

where $\Pi_{A}(k)$ is the transfer of the spectral density $\mathscr{A}(k)$
from $(0,k)$ to $(k,\infty)$, and $D_{A}(k)$ and $F_{A}(k)$ are the
dissipation and forcing spectra of $A$, respectively. In an inertial range
where $A$ is the cascading quantity, then $\Pi_{A}(k)=\varepsilon_{A}$ where
$\varepsilon_{A}$ is a constant and thus $\partial\Pi_{A}(k)/\partial k=0$.

###  3.2 The inverse and forward cascade

For a fluid with the variable stratification inversion relation (14) that is
forced in the wavenumber interval $[k_{1},k_{2}]$, Gkioulekas and Tung (2007)
prove the following two inequalities for stationary turbulence,

$\displaystyle\int_{0}^{k}\frac{\mathrm{d}m(k^{\prime})}{\mathrm{d}k^{\prime}}\,\Pi_{E}(k^{\prime})\,\mathrm{d}k^{\prime}<0,\,\,\text{for
all }k>k_{2},$ (28)
$\displaystyle\int_{k}^{\infty}\frac{\mathrm{d}m(k^{\prime})}{\mathrm{d}k^{\prime}}\,\frac{\Pi_{P}(k^{\prime})}{[m(k^{\prime})]^{2}}\,\mathrm{d}k^{\prime}>0,\,\,\text{for
all }k<k_{1}.$ (29)

These two inequalities do not require the existence of inertial ranges, only
that the inversion function $m(k)$ is an increasing function of $k$.
Therefore, if $\mathrm{d}m(k)/\mathrm{d}k>0$, then there is a net inverse
cascade of total energy and a net forward cascade of surface potential
enstrophy.

###  3.3 When do inertial ranges form?

The lack of scale invariance along with the presence of two length scales,
$k^{-1}$ and $[m(k)]^{-1}$, prevents the use of the Kraichnan (1967) argument
to establish the existence of an inertial range. However, suppose that in a
wavenumber interval, $[k_{a},k_{b}]$, the inversion function takes the power
law form

$m(k)\approx m_{\alpha}\,k^{\alpha},$ (30)

where $m_{\alpha}>0$ and $\alpha>0$. Then, in this wavenumber interval, the
inversion relation takes the form of the $\alpha$-turbulence inversion
relation (Pierrehumbert et al., 1994),

$\hat{\xi}_{\bm{k}}=-k^{\alpha}\,\hat{\psi}_{\bm{k}},$ (31)

with $\xi=\theta/m_{\alpha}$. The inversion relation (31) is then scale
invariant in the wavenumber interval $[k_{a},k_{b}]$. Moreover, $k^{-1}$ is
the only available length scale if the turbulence is sufficiently local in
wavenumber space. It follows that if the wavenumber interval $[k_{a},k_{b}]$
is sufficiently wide (i.e., $k_{a}\ll k_{b}$), then Kraichnan’s argument
applies to the turbulence over this wavenumber interval and inertial ranges
are expected to form.

###  3.4 The Tulloch and Smith (2006) argument

If we assume the existence of inertial ranges, then we can adapt the cascade
argument of Tulloch and Smith (2006) to general surface quasigeostrophic
fluids to obtain predictions for the cascade spectra.

In the inverse cascade inertial range, we must have
$\Pi_{E}(k)=\varepsilon_{E}$ where $\varepsilon_{E}$ is a constant. Assuming
locality in wavenumber space, we have

$\varepsilon_{E}\sim\frac{k\,\mathscr{E}(k)}{\tau(k)},$ (32)

where $\tau(k)$ is a spectrally local timescale444A spectrally local timescale
is appropriate so long as $m(k)$ grows less quickly than $k^{2}$. Otherwise, a
non-local timescale must be used (Kraichnan, 1971; Watanabe and Iwayama,
2004).. If we further assume that the timescale $\tau(k)$ is determined by the
kinetic energy spectrum, $\mathscr{K}(k)$, then dimensional consistency
requires

$\tau(k)\sim\left[k^{3}\,\mathscr{K}(k)\right]^{-1/2}.$ (33)

Substituting this timescale into equation (32) and using the relationship
between the energy spectrum, $\mathscr{E}(k)$, and the streamfunction variance
spectrum, $\mathscr{S}(k)$, in equations (25) and (26), we obtain the total
energy spectrum in the inverse cascade inertial range,

$\mathscr{E}(k)\sim\varepsilon_{E}^{2/3}\,k^{-7/3}\,\left[m(k)\right]^{1/3}.$
(34)

Analogously, in the forward cascade inertial range, we must have
$\Pi_{P}(k)=\varepsilon_{P}$ where $\varepsilon_{P}$ is a constant. A similar
argument yields the surface potential enstrophy spectrum in the forward
cascade inertial range,

$\mathscr{P}(k)\sim\varepsilon_{P}^{2/3}\,k^{-7/3}\,\left[m(k)\right]^{2/3}.$
(35)

The predicted spectra (34) and (35) are not uniquely determined by dimensional
analysis. Rather than assuming that the spectrally local timescale $\tau(k)$
is determined by the kinetic energy spectrum, $\mathscr{K}(k)$, we can assume
that $\tau(k)$ is determined by the total energy spectrum, $\mathscr{E}(k)$,
or the surface potential enstrophy spectrum, $\mathscr{P}(k)$.555These
assumptions lead to timescales of
$\tau(k)\sim\left[k^{4}\,\mathscr{E}(k)\right]^{-1/2}$ and
$\tau(k)\sim\left[k^{3}\,\mathscr{P}(k)\right]^{-1/2}$, respectively. Either
choice will result in cascade spectra distinct from (34) and (35). However, by
assuming that the timescale $\tau(k)$ is determined by the kinetic energy
spectrum, the resulting cascade spectra agree with the $\alpha$-turbulence
predictions of Pierrehumbert et al. (1994) when the inversion function takes
the power law form (30).

For later reference, we provide the expressions for the inverse and forward
cascade surface kinetic energy spectra. Using either the inverse cascade
spectrum (34) or forward cascade spectrum (35) along with the relations
between the various spectra [equations (25) and (26)], we obtain

$\mathscr{K}(k)\sim\varepsilon_{E}^{2/3}\,k^{-1/3}\,\left[m(k)\right]^{-2/3}$
(36)

in the inverse cascade and

$\mathscr{K}(k)\sim\varepsilon_{P}^{2/3}\,k^{-1/3}\,\left[m(k)\right]^{-4/3}$
(37)

in the forward cascade.

Finally, we note that the vorticity spectrum,

$\mathscr{Z}(k)=k^{2}\,\mathscr{K}(k),$ (38)

is an increasing function of $k$ if $m(k)$ is flatter than $k^{5/4}$. In
particular, at small scales, we expect $m(k)\sim k$ (section 22.3), implying a
vorticity spectrum of $\mathscr{Z}(k)\sim k^{1/3}$. Such an increasing
vorticity spectrum implies high Rossby numbers and the breakdown of
geostrophic balance at small scales.

##  4 Idealized stratification profiles

In this section we provide analytical solutions for $m(k)$ in the cases of an
increasing and decreasing piecewise constant stratification profiles as well
as in the case of exponential stratification. These idealized stratification
profiles then provide intuition for the inversion function’s functional form
in the case of an arbitrary stratification profile, $\sigma(z)$.

###  4.1 Piecewise constant stratification

Figure 1: Log-log plots of the inversion function, $m(k)$ [panels (a), (d),
and (g)], for three stratification profiles [panels (b), (e), and (h)] and the
resulting streamfunctions at the two horizontal length scales of 50 km
(dashed) and 100 km (solid) [for panels (c) and (i)] or 2 km and 10 km [panel
(f)]. In the first two inversion function plots [panels (a) and (d)], the thin
solid diagonal lines represent the two linear asymptotic states of
$k/\sigma_{0}$ and $k/\sigma_{\mathrm{pyc}}$. The vertical solid line is the
mixed-layer length scale $L_{\mathrm{mix}}$, given by equation (41), whereas
the vertical dotted line is the pycnocline length scale $L_{\mathrm{pyc}}$,
given by equation (42). The power $\alpha$, where
$m(k)/k^{\alpha}\approx\mathrm{constant}$, is computed by fitting a straight
line to the log-log plot of $m(k)$ between $2\pi/L_{\mathrm{mix}}$ and
$2\pi/L_{\mathrm{pyc}}$. This straight line is shown as a grey dashed line in
panels (a) and (d). In panel (g), the thin diagonal line is the linear small-
scale limit, $m(k)\approx k/\sigma_{0}$, whereas the thin horizontal line is
the constant large-scale limit, $m(k)=2/(\sigma_{0}^{2}\,h_{\mathrm{exp}})$.
Finally, the solid vertical lines in panel (g) indicate the horizontal length
scale $L_{\mathrm{exp}}=2\pi/k_{\mathrm{exp}}$ [equation (48)] induced by the
exponential stratification. Further details on the stratification profiles are
in the text.

Consider the piecewise constant stratification profile, given by

$\sigma(z)=\begin{cases}\sigma_{0}\,&\text{for }-h<z\leq 0\\\
\sigma_{\mathrm{pyc}}\,&\text{for }\infty<z\leq-h.\end{cases}$ (39)

This stratification profile consists of an upper layer of thickness $h$ with
constant stratification $\sigma_{0}$ overlying an infinitely deep layer with
constant stratification $\sigma_{\mathrm{pyc}}$. If
$\sigma_{0}<\sigma_{\mathrm{pyc}}$, then this stratification profile is an
idealization of a weakly stratified mixed-layer overlying an ocean of stronger
stratification. See panels (b) and (e) in figure 1 for an illustration.

For this stratification profile, an analytical solution is possible, with the
solution provided in appendix B. The resulting inversion function is

$m(k)=\frac{k}{\sigma_{0}}\left[\frac{\cosh\left(\sigma_{0}hk\right)+\left(\frac{\sigma_{\mathrm{pyc}}}{\sigma_{0}}\right)\sinh\left(\sigma_{0}hk\right)}{\sinh\left(\sigma_{0}hk\right)+\left(\frac{\sigma_{\mathrm{pyc}}}{\sigma_{0}}\right)\cosh\left(\sigma_{0}hk\right)}\right].$
(40)

At small horizontal scales, $2\pi/k\ll L_{\mathrm{mix}}$, where

$L_{\mathrm{mix}}=2\,\pi\,\sigma_{0}\,h,$ (41)

the inversion function takes the form $m(k)\approx k/\sigma_{0}$, as expected
from the uniformly stratified theory (Held et al., 1995). At large horizontal
scales, $2\pi/k\gg L_{\mathrm{pyc}}$, where

$L_{\mathrm{pyc}}=2\,\pi\begin{cases}\,\sigma_{\mathrm{pyc}}\,h\,&\text{if
}\sigma_{0}\leq\sigma_{\mathrm{pyc}}\\\
\sigma_{\mathrm{0}}^{2}\,h/\sigma_{\mathrm{pyc}}\,&\text{if
}\sigma_{0}>\sigma_{\mathrm{pyc}},\end{cases}$ (42)

then the inversion function takes the form $m(k)\approx
k/\sigma_{\mathrm{pyc}}$, because at large horizontal scales, the ocean will
seem to have constant stratification $\sigma_{\mathrm{pyc}}$.

The functional form of the inversion function at horizontal scales between
$L_{\mathrm{mix}}$ and $L_{\mathrm{pyc}}$ depends on whether $\sigma(z)$ is a
decreasing or increasing function. If $\sigma(z)$ is a decreasing function,
with $\sigma_{0}<\sigma_{\mathrm{pyc}}$, then we obtain a mixed-layer like
stratification profile and the inversion function steepens to a super linear
wavenumber dependence at these scales. An example is shown in figure 1(a)-(b).
Here, the stratification abruptly jumps from a value of $\sigma_{0}\approx 14$
to $\sigma_{\mathrm{pyc}}=100$ at $z\approx-79$ m. Consequently, the inversion
function takes the form $m(k)\sim k^{1.57}$ between $2\pi/L_{\mathrm{pyc}}$
and $2\pi/L_{\mathrm{mix}}$. In contrast, if
$\sigma_{0}>\sigma_{\mathrm{pyc}}$ then the inversion function flattens to a
sublinear wavenumber dependence for horizontal scales between
$L_{\mathrm{mix}}$ and $L_{\mathrm{pyc}}$. An example is shown in figure
1(d)-(e), where the stratification abruptly jumps from $\sigma_{0}\approx 14$
to $\sigma_{\mathrm{pyc}}\approx 2$ at $z\approx-79$ m. In this case, the
inversion function has a sublinear wavenumber dependence, $m(k)\sim k^{0.43}$,
between $2\pi/L_{\mathrm{pyc}}$ and $2\pi/L_{\mathrm{mix}}$.

Figure 2: Results of three pseudo-spectral simulations, forced at
approximately 100 km, with $1024^{2}$ horizontal grid points. See appendix C
for a description of the numerical model. The first simulation [panels (a),
(d), and (g)] corresponds to the stratification profile and inversion function
shown in figure 1(a)-(b), the second simulations [panels (b), (e), and (h)]
corresponds to the stratification profile and inversion function shown in
figure 1(d)-(e), and the third simulation corresponds to the stratification
profile and inversion function shown in figure 1(g)-(h). Plots (a), (b), and
(c) are snapshots of the surface potential vorticity, $\theta$, normalized by
its maximum value in the snapshot. Plots (d), (e), and (f) are snapshots of
the horizontal speed $\left\lvert\bm{u}\right\rvert$ normalized by its maximum
value in the snapshot. Plots (g), (h), and (i) show the model kinetic energy
spectrum (solid black line) along with the prediction given by equation (37)
(dashed black line). We also provide linear fits to the model kinetic energy
spectrum (dash-dotted red line) and to the predicted spectrum (dotted blue
line).

By fitting a power law, $k^{\alpha}$, to the inversion function, we do not
mean to imply that $m(k)$ indeed takes the form of a power law. Instead, the
purpose of obtaining the estimated power $\alpha$ is to apply the intuition
gained from $\alpha$-turbulence (Pierrehumbert et al., 1994; Smith et al.,
2002; Sukhatme and Smith, 2009; Burgess et al., 2015; Foussard et al., 2017)
to surface quasigeostrophic turbulence. In $\alpha$-turbulence, an active
scalar $\xi$, defined by the power law inversion relation (31), is materially
conserved in the absence of forcing and dissipation [that is, $\xi$ satisfies
the time-evolution equation (6) with $\theta$ replaced by $\xi$]. The scalar
$\xi$ can be thought of as a generalized vorticity; if $\alpha=2$ we recover
the vorticity of two-dimensional barotropic model. If $\alpha=1$, $\xi$ is
proportional to the surface buoyancy in the uniformly stratified surface
quasigeostrophic model. To discern how $\alpha$ modifies the dynamics, we
consider a point vortex $\xi\sim\delta(r)$, where $r$ is the horizontal
distance from the vortex and $\delta(r)$ is the Dirac delta. If $\alpha=2$, we
obtain $\psi(r)\sim\log(r)/2\pi$; otherwise, if $0<\alpha<2$, we obtain
$\psi(r)\sim-C_{\alpha}/r^{2-\alpha}$ where $C_{\alpha}>0$ (Iwayama and
Watanabe, 2010). Therefore, larger $\alpha$ leads to vortices with a longer
interaction range whereas smaller $\alpha$ leads to a shorter interaction
range.

More generally, $\alpha$ controls the spatial locality of the resulting
turbulence. In two-dimensional turbulence ($\alpha=2$), vortices induce flows
reaching far from the vortex core and the combined contributions of distant
vortices dominates the local fluid velocity. These flows are characterized by
thin filamentary $\xi$-structures due to the dominance of large-scale strain
(Watanabe and Iwayama, 2004). As we decrease $\alpha$, the turbulence becomes
more spatially local, the dominance of large-scale strain weakens, and a
secondary instability becomes possible in which filaments roll-up into small
vortices; the resulting turbulence is distinguished by vortices spanning a
wide range of horizontal scales, as in uniform stratification surface
quasigeostrophic turbulence (Pierrehumbert et al., 1994; Held et al., 1995).
As $\alpha$ is decreased further the $\xi$ field becomes spatially diffuse
because the induced velocity, which now has small spatial scales, is more
effective at mixing small-scale inhomogeneities in $\xi$ (Sukhatme and Smith,
2009).

These expectations are confirmed in the simulations shown in figure 2. The
simulations are set in a doubly periodic square with side length 400 km and
are forced at a horizontal scale of 100 km. Large-scale dissipation is
achieved through a linear surface buoyancy damping whereas an exponential
filter is applied at small scales. In the case of a mixed-layer like
stratification, with $\sigma_{0}<\sigma_{\mathrm{pyc}}$, the $\theta$ field
exhibits thin filamentary structures (characteristic of the $\alpha=2$ case)
as well as vortices spanning a wide range of horizontal scales (characteristic
of the $\alpha=1$ case). In contrast, although the
$\sigma_{0}>\sigma_{\mathrm{pyc}}$ simulation exhibits vortices spanning a
wide range of scales, no large-scale filaments are evident. Instead, we see
that the surface potential vorticity is spatially diffuse. These contrasting
features are consequences of the induced horizontal velocity field. The mixed-
layer like case has a velocity field dominated by large-scale strain, which is
effective at producing thin filamentary structures. In contrast the velocity
field in the $\sigma_{0}>\sigma_{\mathrm{pyc}}$ case consists of narrow
meandering currents, which are effective at mixing away small-scale
inhomogeneities.

Both the predicted [equation (37)] and diagnosed surface kinetic energy
spectra are plotted in figure 2. In the mixed-layer like case, with
$\sigma_{0}<\sigma_{\mathrm{pyc}}$, the predicted and diagnosed spectrum are
close, although the diagnosed spectrum is steeper at large scales (a too steep
spectrum is also observed in the $\alpha=1$ and $\alpha=2$ cases, see
Schorghofer, 2000). In the $\sigma_{0}>\sigma_{\mathrm{pyc}}$ case, the large-
scale spectrum agrees with the predicted spectrum. However, at smaller scales,
the model spectrum is significantly steeper.

Figure 3: The spectral density transfer functions for surface potential
enstrophy (a) and surface kinetic energy (b) normalized by their absolute
maximum for the three simulations in figure 2.

The derivation of the predicted spectra in section 3 assumed the existence of
an inertial range, which in this case means $\Pi_{P}(k)=$ constant. To verify
whether this assumption holds, we show in figure 3(a) the spectral transfer of
surface potential enstrophy for both the $\sigma_{0}<\sigma_{\mathrm{pyc}}$
and the $\sigma_{0}>\sigma_{\mathrm{pyc}}$ cases. In the mixed-layer like
case, with $\sigma_{0}<\sigma_{\mathrm{pyc}}$, an approximate inertial range
forms with some deviations at larger scales. However, in the
$\sigma_{0}>\sigma_{\mathrm{pyc}}$ case, $\Pi_{P}$ is an increasing function
at small scales, which indicates that the spectral density of surface
potential enstrophy, $\mathscr{P}(k)$, is diverging at these scales. That is,
at small scales, there is a depletion of $\mathscr{P}(k)$ and this depletion
is causing the steepening of the kinetic energy spectrum at small scales in
figure 2.

The steepening of the model spectrum at small scales may be a consequence of
the inversion function being steeper at small scales than at large scales. We
tested this hypothesis with two additional simulations having a prescribed
inversion function of the form

$m(k)=\begin{cases}k^{\alpha_{1}}\quad\text{for }k<k_{0}\\\
C\,k^{\alpha_{2}}\quad\text{for }k>k_{0},\end{cases}$ (43)

where $\alpha_{1}$,$\alpha_{2}$, and $k_{0}$ are positive numbers, and $C$ is
chosen to ensure that $m(k)$ is continuous. In the first simulation, we chose
$\alpha_{1}=3/2$ and $\alpha_{2}=1/2$ so that $m(k)$ is flatter at small
scales than at large scales (as in the mixed-layer like case). We obtained an
approximate inertial range and the model spectrum is close to the predicted
surface kinetic energy spectrum (37). In the second simulation, we chose
$\alpha_{1}=1/2$ and $\alpha_{2}=3/2$ so that $m(k)$ is steeper at small
scales than at large scales (as in the $\sigma_{0}>\sigma_{\mathrm{pyc}}$
case). We found that $\Pi_{P}(k)$ is an increasing function (so that no
inertial range exists) and we obtained a model surface kinetic energy spectrum
that is much steeper than predicted. It is not clear why the inertial range
theory fails in this case, and the failure may be a consequence of the finite
model resolution.

###  4.2 An exponentially stratified ocean

Now consider the exponential stratification profile

$\sigma=\sigma_{0}\,\mathrm{e}^{z/h_{\mathrm{exp}}}.$ (44)

Substituting the stratification profile (44) into the vertical structure
equation (9) with boundary conditions (10) and (11) yields the vertical
structure

$\Psi_{k}(z)=\mathrm{e}^{z/h_{\mathrm{exp}}}\,\frac{I_{1}\left(\mathrm{e}^{z/h_{\mathrm{exp}}}\sigma_{0}\,h_{\mathrm{exp}}\,k\right)}{I_{1}\left(\sigma_{0}\,h_{\mathrm{exp}}\,k\right)},$
(45)

where $I_{n}(z)$ is the modified Bessel function of the first kind of order
$n$.

To obtain the inversion function, we substitute the vertical structure (45)
into the definition of the inversion function (15) to obtain

$\displaystyle m(k)=$
$\displaystyle\frac{1}{\sigma_{0}^{2}h_{\mathrm{exp}}}\,+$ (46)
$\displaystyle\frac{k}{2\sigma_{0}}\left[\frac{I_{0}\left(\sigma_{0}h_{\mathrm{exp}}k\right)}{I_{1}\left(\sigma_{0}h_{\mathrm{exp}}k\right)}+\frac{I_{2}\left(\sigma_{0}h_{\mathrm{exp}}k\right)}{I_{1}\left(\sigma_{0}h_{\mathrm{exp}}k\right)}\right].$

In the small-scale limit $k\gg 1/\left(\sigma_{0}\,h_{\mathrm{exp}}\right)$,
the inversion function becomes $m(k)\approx k/\sigma_{0}$ as in constant
stratification surface quasigeostrophic theory. In contrast, the large-scale
limit $k\ll 1/\left(\sigma_{0}\,h_{\mathrm{exp}}\right)$ gives

$m(k)\approx\frac{h_{\mathrm{exp}}}{4}\left(k_{\mathrm{exp}}^{2}+k^{2}\right),$
(47)

where $k_{\mathrm{exp}}$ is given by

$k_{\mathrm{exp}}=\frac{2\,\sqrt{2}}{\sigma_{0}\,h_{\mathrm{exp}}}.$ (48)

As $k/k_{\mathrm{exp}}\rightarrow 0$, the inversion function asymptotes to a
constant value and the vertical structure becomes independent of the
horizontal scale $2\pi/k$, with $\Psi_{k}\rightarrow\Psi_{0}$ where

$\Psi_{0}(z)=\mathrm{e}^{2z/h_{\mathrm{exp}}}.$ (49)

Further increasing the horizontal scale no longer modifies $\Psi_{k}(z)$ and
so vertical structure is arrested at $\Psi_{0}$.

An example with $h_{\mathrm{exp}}=300$ m and $\sigma_{0}=100$ is shown in
figure 1(g)-(i). At horizontal scales smaller than
$L_{\mathrm{exp}}=2\pi/k_{\mathrm{exp}}$, the inversion function rapidly
transitions to the linear small-scale limit of $m(k)\approx k/\sigma_{0}$. In
contrast, at horizontal scales larger than $L_{\mathrm{exp}}$, the large-scale
approximation (47) holds, and at sufficiently large horizontal scales, the
inversion function asymptotes to constant value of
$m(k)=h_{\mathrm{exp}}\,k_{\mathrm{exp}}^{2}/4$.

The inversion relation implied by the inversion function (47) is

$\hat{\theta}_{\bm{k}}\approx-\frac{h_{\mathrm{exp}}}{4}\left(k_{\mathrm{exp}}^{2}+k^{2}\right)\hat{\psi}_{\bm{k}},$
(50)

which is isomorphic to the inversion relation in the equivalent barotropic
model (Larichev and McWilliams, 1991), with $k_{\mathrm{exp}}$ assuming the
role of the deformation wavenumber. In this limit, the total energy and the
surface potential enstrophy are not independent quantities to leading order in
$k_{\mathrm{exp}}^{2}$. Using the relations between the various spectra
[equations (25) and (26)] with an inversion function of the form $m(k)\approx
m_{0}+m_{1}k^{2}$, we obtain $\mathscr{E}(k)\approx
m_{0}\,\mathscr{S}(k)+m_{1}\,\mathscr{K}(k)$ and $\mathscr{P}(k)\approx
m_{0}^{2}\,\mathscr{S}(k)+2\,m_{0}\,m_{1}\mathscr{K}(k)$; solving for
$\mathscr{S}(k)$ and $\mathscr{K}(k)$ then yields

$\mathscr{S}(k)\approx\frac{2\,m_{0}\mathscr{E}(k)-\mathscr{P}(k)}{m_{0}^{2}},$
(51)

and

$\mathscr{K}(k)\approx\frac{\mathscr{P}(k)-m_{0}\,\mathscr{E}(k)}{m_{0}\,m_{1}},$
(52)

which are now the two independent quantities. Then, using an argument
analogous to that in Larichev and McWilliams (1991), we find that

$\mathscr{S}(k)\sim k^{-11/3}$ (53)

in the inverse cascade inertial range whereas

$\mathscr{K}(k)\sim k^{-3}$ (54)

in the forward cascade inertial range.

The implied dynamics are extremely local; a point vortex,
$\theta(r)\sim\delta(r)$, leads to an exponentially decaying streamfunction,
$\psi(r)\sim\exp(-k_{\mathrm{exp}}r)/\sqrt{k_{\mathrm{exp}}r}$ (Polvani et
al., 1989). Therefore, as for the $\sigma_{0}>\sigma_{\mathrm{pyc}}$ case
above, we expect a spatially diffuse surface potential vorticity field and no
large-scale strain. However, unlike the $\sigma_{0}>\sigma_{\mathrm{pyc}}$
case, the presence of a distinguished length scale, $L_{\mathrm{exp}}$, leads
to the emergence of plateaus of homogenized surface potential vorticity
surrounded by kinetic energy ribbons (Arbic and Flierl, 2003). Both of these
features can be seen in figure 2.

The $k^{-3}$ surface kinetic energy spectrum (54) is only expected to hold at
horizontal scales larger than $2\,\pi\,\sigma_{0}\,h_{\mathrm{exp}}$; at
smaller scales we should recover the $k^{-5/3}$ spectrum expected from
uniformly stratified surface quasiogeostrophic theory. Figure 2(i) shows that
there is indeed a steepening of the kinetic energy spectrum at horizontal
scales larger than 20 km, although the model spectrum is somewhat steeper than
the predicted $k^{-3}$. Similarly, although the spectrum flattens at smaller
scales, the small-scale spectrum is also slightly steeper than the predicted
$k^{-5/3}$.

We can also examine the spectral transfer functions of $\mathscr{P}(k)$ and
$\mathscr{K}(k)$. At large scales, we expect an inertial range in surface
kinetic energy, so $\Pi_{K}(k)=$ constant, whereas at small scales, we expect
an inertial range in surface potential enstrophy, so $\Pi_{P}(k)=$ constant.
However, figure 3 shows that although both $\Pi_{K}(k)$ and $\Pi_{P}(k)$
become approximately flat at small scales, we observe significant deviations
at larger scales.

###  4.3 More general stratification profiles

These three idealized cases provide intuition for how the inversion function
behaves for an arbitrary stratification profile, $\sigma(z)$. Generally, if
$\sigma(z)$ is decreasing over some depth, then the inversion function will
steepen to a super linear wavenumber dependence over a range of horizontal
wavenumber whose vertical structure function significantly impinges on these
depths. A larger difference in stratification between these depths leads to a
steeper inversion function. Analogously, if $\sigma(z)$ is increasing over
some depth, then the inversion function will flatten to a sublinear wavenumber
dependence, with a larger difference in stratification leading to a flatter
inversion function. Finally, if $\sigma(z)$ is much smaller at depth than near
the surface, the inversion function will flatten to become approximately
constant, and we recover an equivalent barotopic like regime, similar to the
exponentially stratified example.

##  5 Application to the ECCOv4 ocean state estimate

We now show that, over the mid-latitude North Atlantic, the inversion function
is seasonal at horizontal scales between 1-100 km, transitioning from
$m(k)\sim k^{3/2}$ in winter to $m(k)\sim k^{1/2}$ in summer. To compute the
inversion function $m(k)$, we obtain the stratification profile
$\sigma(z)=N(z)/f$ at each location from the Estimating the Circulation and
Climate of the Ocean version 4 release 4 (ECCOv4, Forget et al., 2015) state
estimate. We then compute $\Psi_{k}(z)$ using the vertical structure equation
(9) and then use the definition of the inversion function (15) to obtain
$m(k)$ at each wavenumber $k$.

###  5.1 The three horizontal length-scales

In addition to $L_{\mathrm{mix}}$ and $L_{\mathrm{pyc}}$ [defined in equations
(41) and (42)], we introduce the horizontal length scale, $L_{H}$, the full-
depth horizontal scale, defined by

$L_{H}=2\,\pi\,\sigma_{\mathrm{ave}}\,H,$ (55)

where $\sigma_{\mathrm{ave}}$ is the vertical average of $\sigma$ and $H$ is
the local ocean depth. The bottom boundary condition becomes important to the
dynamics at horizontal scales larger than $\approx L_{H}$.

Figure 4: Panels (a), (b), and (c) show the horizontal length scales $L_{H}$,
$L_{\textrm{mix}}$, and $L_{\textrm{pyc}}$ as computed from 2017 January mean
ECCOv4 stratification profiles, $\sigma(z)=N(z)/f$, over the North Atlantic.
The green ‘x’ in panel (a) shows the location chosen for the inversion
functions of figure 6 and the model simulations of figure 7. Panel (d) shows
$\alpha$, defined by $m(k)/k^{\alpha}\approx\mathrm{constant}$, over the North
Atlantic. We compute $\alpha$ by fitting a straight line to a log-log plot of
$m(k)$ between $2\pi/L_{\mathrm{mix}}$ and $2\pi/L_{\mathrm{pyc}}$. Panel (e)
is a histogram of the computed values of $\alpha$ over the North Atlantic. We
exclude from this histogram grid cells with $L_{H}<150$ km; these are
primarily continental shelves and high-latitude regions. In these excluded
regions, our chosen bottom boundary condition (56) may be influencing the
computed value of $\alpha$.

Figure 5: Panels (a), (b), (c), and (e) are as in figure 4(a)-(d), but
computed from 2017 July mean stratification profiles. The calculation of
$L_{0}$ in panel (d) is explained in the text. In panel (f), we show $\alpha$
but measured between $2\pi/(50\,\mathrm{km})$ and $2\pi/L_{0}$.

We compute all three length scales using ECCOv4 stratification profiles over
the North Atlantic, with results displayed in figures 4(a)-(c) and 5(a)-(c)
for January and July, respectively. To compute the mixed-layer horizontal
scale, $L_{\mathrm{mix}}=2\,\pi\sigma_{0}\,h_{\mathrm{mix}}$, we set
$\sigma_{0}$ equal to the stratification at the uppermost grid cell. The
mixed-layer depth, $h_{\mathrm{mix}}$, is then defined as follows. We first
define the pycnocline stratification, $\sigma_{\mathrm{pyc}}$, to be the
maximum of $\sigma(z)$. The mixed-layer depth, $h_{\mathrm{mix}}$, is then the
depth at which
$\sigma(-h_{\mathrm{mix}})=\sigma_{0}+\left(\sigma_{\mathrm{pyc}}-\sigma_{0}\right)/4$.
Finally, the pycnocline horizontal scale, $L_{\mathrm{pyc}}$, is computed as
$L_{\mathrm{pyc}}=2\,\pi\,\sigma_{\mathrm{pyc}}\,h_{\mathrm{pyc}}$, where
$h_{\mathrm{pyc}}$ is the depth of the stratification maximum
$\sigma_{\mathrm{pyc}}$.

Figures 4(a) and 5(a) show that $L_{H}$ is not seasonal, with typical mid-
latitude open ocean values between $400-700$ km. On continental shelves, as
well as high-latitudes, $L_{H}$ decreases to values smaller than $200$ km. As
we approach the equator, the full-depth horizontal scale $L_{H}$ becomes large
due to the smallness of the Coriolis parameter.

Constant stratification surface quasigeostrophic theory is only valid at
horizontal scales smaller than $L_{\mathrm{mix}}$. Figure 4(b) shows that the
wintertime $L_{\mathrm{mix}}$ is spatially variable with values ranging
between $1-15$ km. In contrast, figure 5(b) shows that the summertime
$L_{\mathrm{mix}}$ is less than 2 km over most of the midlatitude North
Atlantic.

Finally, we expect to observe a superlinear inversion function for horizontal
scales between $L_{\mathrm{mix}}$ and $L_{\mathrm{pyc}}$. The latter,
$L_{\mathrm{pyc}}$, is shown in figures 4(c) and 5(c). Typical mid-latitude
values range between $70-110$ km in winter but decrease to $15-30$ km in
summer. In figure 4(c), there is a region close to the Gulf stream with
anomalously high values of $L_{\mathrm{pyc}}$. In this region, the
stratification maximum, $\sigma_{\mathrm{pyc}}$ is approximately half as large
as the surrounding region, but its depth, $h_{\mathrm{pyc}}$, is about five
times deeper, resulting in larger values of
$L_{\mathrm{pyc}}=2\,\pi\,\sigma_{\mathrm{pyc}}\,h_{\mathrm{pyc}}$.

###  5.2 The inversion function at a single location

Before computing the inversion function over the North Atlantic, we focus on a
single location. However, we must first address what boundary conditions to
use in solving the vertical structure equation $\eqref{eq:Psi_equation}$ for
$\Psi_{k}(z)$. We cannot use the infinite bottom boundary condition (11)
because the ocean has a finite depth. However, given that figures 4(a) and
5(a) show that the bottom boundary condition should not effect the inversion
function at horizontal scales smaller than 400 km in the mid-latitude North
Atlantic open ocean, we choose to use the no-slip bottom boundary condition

$\displaystyle\Psi_{k}(-H)=0.$ (56)

The alternate free-slip boundary condition

$\displaystyle\frac{\mathrm{d}\Psi_{k}(-H)}{\mathrm{d}z}=0$ (57)

gives qualitatively identical results for horizontal scales smaller than 400
km, which are the scales of interest in this study [see appendix A for the
large-scale limit of $m(k)$ under these boundary conditions].666 The no-slip
boundary condition (56) is appropriate over strong bottom friction (Arbic and
Flierl, 2004) or steep topography (LaCasce, 2017) whereas the free-slip
boundary condition (57) is appropriate over a flat bottom.

Figure 6: As in figure 1 but for the mid-latitude North Atlantic location
($38^{\circ}$ N, $45^{\circ}$ W) in January [(a)-(c)] and July [(d)-(f)]. This
location is marked by a green ‘x’ in figure 4(a). Only the upper 750 m of the
stratification profiles and vertical structures are shown in panels (b), (c),
(e) and (f).

Figure 7: Two pseudo-spectral simulations differing only in the chosen
stratification profile $\sigma(z)=N(z)/f$. Both simulations use a monthly
averaged 2017 stratification at the mid-latitude North Atlantic location
($38^{\circ}$ N,$45^{\circ}$ W) [see the green ’x’ in figure 4a] in January
[(a), (c), (e)] and July [(b), (d), (f)]. The stratification profiles are
obtained from the Estimating the Circulation and Climate of the Ocean version
4 release 4 (ECCOv4, Forget et al., 2015) state estimate. Otherwise as in
figure 2.

Figure 6 shows the computed inversion function in the mid-latitude North
Atlantic at ($38^{\circ}$ N, $45^{\circ}$ W) [see the green ‘x’ in figure
4(a)]. In winter, at horizontal scales smaller than $L_{\mathrm{mix}}\approx
5$ km, we recover the linear $m(k)\approx k/\sigma_{0}$ expected from constant
stratification surface quasigeostrophic theory. However, for horizontal scales
between $L_{\mathrm{mix}}\approx 5$ km and $L_{\mathrm{pyc}}\approx 70$ km,
the inversion function, $m(k)$, becomes as steep as a $k^{3/2}$ power law.
Figure 7 shows a snapshot of the surface potential vorticity and the
geostrophic velocity from a surface quasigeostrophic model using the
wintertime inversion function. The surface potential vorticity snapshot is
similar to the idealized mixed-layer snapshot of figure 2(a), which is also
characterized by $\alpha\approx 3/2$ (but at horizontal scales between 7-50
km). Both simulations exhibit a preponderance of small-scale vortices as well
as thin filaments of surface potential vorticity. As expected, the kinetic
energy spectrum [figure 7(e)] transitions from an $\alpha\approx 3/2$ regime
to an $\alpha=1$ regime near $L_{\mathrm{mix}}=5$ km. Moreover, as shown in
figure 8, an approximate inertial range is evident between the forcing and
dissipation scales.

Figure 8: The spectral density transfer functions for surface potential
enstrophy (a) and surface kinetic energy (b) normalized by their absolute
maximum for the two simulations in figure 7.

In summer, the mixed-layer horizontal scale, $L_{\mathrm{mix}}$, becomes
smaller than 1 km and the pycnocline horizontal scale, $L_{\mathrm{pyc}}$,
decreases to 20 km. We therefore obtain a super linear regime, with $m(k)$ as
steep as $k^{1.2}$, but only for horizontal scales between 1-20 km. Thus,
although there is a range of wavenumbers for which $m(k)$ steepens to a super
linear wavenumber dependence in summer, this range of wavenumbers is narrow,
only found at small horizontal scales, and the steepening is much less
pronounced than in winter. At horizontal scales larger than
$L_{\mathrm{pyc}}$, the summertime inversion function flattens, with the
$m(k)$ increasing like a $k^{1/2}$ power law between 50-400 km. This
flattening is due to the largely decaying nature of ocean stratification below
the stratification maximum.

As expected from a simulation with a sublinear inversion function at large
scales, the surface potential vorticity appears spatially diffuse [figure
7(d)] and comparable to the $\sigma_{0}>\sigma_{\mathrm{pyc}}$ and the
exponential simulations [figure 2(b)-(c)]. However, despite having a sublinear
inversion function, the July simulation is dynamically more similar to the
exponential simulation rather than the $\sigma_{0}>\sigma_{\mathrm{pyc}}$
simulation. The July simulation displays approximately homogenized regions of
surface potential vorticity surrounded by surface kinetic energy ribbons. As a
result, the surface kinetic energy spectrum does not follow the predicted
spectrum (37).

###  5.3 The inversion function over the North Atlantic

We now present power law approximations to the inversion function $m(k)$ over
the North Atlantic in winter and summer. In winter, we obtain the power
$\alpha$, where $m(k)/k^{\alpha}\approx\mathrm{constant}$, by fitting a
straight line to $m(k)$ on a log-log plot between $2\pi/L_{\mathrm{mix}}$ and
$2\pi/L_{\mathrm{pyc}}$. A value of $\alpha=1$ is expected for constant
stratification surface quasigeostrophic theory. A value of $\alpha=2$ leads to
an inversion relation similar to two-dimensional barotropic dynamics. However,
in general, we emphasize that $\alpha$ is simply a crude measure of how
quickly $m(k)$ is increasing; we do not mean to imply that $m(k)$ in fact
takes the form of a power law. Nevertheless, the power $\alpha$ is useful
because, as $\alpha$-turbulence suggests (and the simulations in section 4
confirm), the rate of increase of the inversion function measures the spatial
locality of the resulting flow.

Figure 9: Panel (a) is as in figure 4(e), but with the additional restriction
that $L_{H}<750$ km to filter out the non-seasonal equatorial region. In panel
(b), we instead plot $\alpha$ as obtained by fitting a straight line to a log-
log plot of $m(k)$ between $2\pi/(50\,\mathrm{km})$ and $2\pi/L_{0}$ with the
same restrictions as in panel (a).

Figure 4(d) shows that we generally have $\alpha\approx 3/2$ in the wintertime
open ocean. Deviations appear at high-latitudes (e.g., the Labrador sea and
southeast of Greenland) and on continental shelves where we find regions of
low $\alpha$. However, both of these regions have small values of $L_{H}$ so
that our chosen no-slip bottom boundary condition (56) may be influencing the
computed $\alpha$ there.

A histogram of the computed values of $\alpha$ [figure 4(e)] confirms that
$\alpha\approx 1.53\pm 0.08$ in the wintertime mid-latitude open ocean. This
histogram only includes grid cells with $L_{H}>150$ km, which ensures that the
no-slip bottom boundary condition (56) is not influencing the computed
distribution.

An inversion function of $m(k)\sim k^{3/2}$ implies a surface kinetic energy
spectrum of $k^{-4/3}$ upscale of small-scale forcing [equation (36)] and a
spectrum of $k^{-7/3}$ downscale of large-scale forcing [equation (37)]. As we
expect wintertime surface buoyancy anomalies to be forced both by large-scale
baroclinic instability and by small-scale mixed-layer baroclinic instability,
the realized surface kinetic energy spectrum should be between $k^{-4/3}$ and
$k^{-7/3}$. Such a prediction is consistent with the finding that wintertime
North Atlantic geostrophic surface velocities are mainly due to surface
buoyancy anomalies (Lapeyre, 2009; González-Haro and Isern-Fontanet, 2014) and
observational evidence of a $k^{-2}$ wintertime spectrum (Callies et al.,
2015).

The universality of the $m(k)\sim k^{3/2}$ regime over the mid-latitudes is
expected because it arises from a mechanism universally present over the mid-
latitude ocean in winter; namely, the deepening of the mixed-layer. However, a
comment is required on why this regime also appears at low latitudes where we
do not observe deep wintertime mixed-layers. At low latitudes, the $m(k)\sim
k^{3/2}$ regime emerges because there is a large scale-separation between
$L_{\mathrm{mix}}$ and $L_{\mathrm{pyc}}$. The smallness of the low latitude
Coriolis parameter $f$ cancels out the shallowness of the low latitude mixed-
layer depth resulting in values of $L_{\mathrm{mix}}$ comparable to the
remainder of the mid-latitude North Atlantic, as seen in figure 4(b). However,
no similar cancellation occurs for $L_{\mathrm{pyc}}$ which reaches values of
$\approx 500$ km due to the smallness of the Coriolis parameter $f$ at low
latitudes. As a consequence, there is a non-seasonal $m(k)\sim k^{3/2}$ regime
at low latitudes for horizontal scales between $10-500$ km.

The analogous summertime results are presented in figure 5(e) and figure 9(a).
Near the equator, we obtain values close to $\alpha\approx 3/2$, as expected
from the weak seasonality there. In contrast, the midlatitudes generally
display $\alpha\approx 1.2-1.3$ but this superlinear regime is only present at
horizontal scales smaller than $L_{\mathrm{pyc}}\approx 20-30$ km. Figure 9(a)
shows a histogram of the measured $\alpha$ values but with the additional
restriction that $L_{H}<750$ km to filter out the near equatorial region
(where $\alpha\approx 3/2$).

The summertime inversion function shown in figure 6(d) suggests that the
inversion function flattens at horizontal scales larger than 50 km, with
$m(k)$ increasing like a $k^{1/2}$ power law. We now generalize this
calculation to the summertime midlatitude North Atlantic by fitting a straight
line to $m(k)$ on a log-log plot between $2\pi/(50\,\mathrm{km})$ and
$2\pi/L_{0}$ where $L_{0}$ is defined by

$m\left(\frac{2\,\pi}{L_{0}}\right)=m_{0}=\left[\int_{-H}^{0}\sigma^{2}(s)\mathrm{d}s\right]^{-1}$
(58)

and $m_{0}$ is defined by the second equality. In this case, we solve for
$m(k)$ using the free-slip boundary condition (57). We made this choice
because $m(k)$ must cross $m_{0}$ in the large-scale limit if we apply the
free-slip boundary condition (57). In contrast, $m(k)$ asymptotes to $m_{0}$
from above if we apply the no-slip boundary condition (56). See appendix A for
more details. In any case, if we use the free-slip boundary condition (57),
then $L_{0}$ is a horizontal length scale at which the flattening of $m(k)$
ceases and $m(k)$ instead begins to steepen in order to attain the required
$H\,k^{2}$ dependence at large horizontal scales [see equation (65)]. Over the
mid-latitudes North Atlantic, $L_{0}$ has typical values of 200-500 km [figure
5(d)].

When $\alpha$ is measured between 50 km and $L_{0}$, we find typical
midlatitude values close to $\alpha\approx 1/2$ [figure 5(f)]. A histogram of
these $\alpha$ values is provided in figure 9(b), where we only consider grid
cells satisfying $L_{H}>150$ km and $L_{H}<750$ km (the latter condition
filters out near equatorial grid cells). The distribution is broad with a mean
of $\alpha=0.56\pm 0.15$ and a long tail of $\alpha>0.8$ values. Therefore,
$m(k)$ flattens considerably in response to the decaying nature of summertime
upper ocean stratification. It is not clear, however, whether the resulting
dynamics will be similar to the $\sigma_{0}>\sigma_{\mathrm{pyc}}$ case or the
exponentially stratified case in section 4. As we have seen, the summertime
simulation (in figure 7) displayed characteristics closer to the idealized
exponential case than the $\sigma_{0}>\sigma_{\mathrm{pyc}}$ case.
Nevertheless, the low summertime values of $\alpha$ indicate that buoyancy
anomalies generate shorter range velocity fields in summer than in winter.

Isern-Fontanet et al. (2014) and González-Haro et al. (2020) measured the
inversion function empirically, through equation (1), and found that the
inversion function asymptotes to a constant at large horizontal scales (270 km
near the western coast of Australia and 100 km in the Mediterranean Sea). They
suggested this flattening is due to the dominance of the interior
quasigeostrophic solution at large scales (a consequence of equation 29 in
Lapeyre and Klein, 2006). We instead suggest this flattening is intrinsic to
surface quasigeostrophy. In our calculation, the inversion function does not
become constant at horizontal scales smaller than 400 km. However, if the
appropriate bottom boundary condition is the no-slip boundary condition (56),
then the inversion asymptotes to a constant value at horizontal scales larger
than $L_{H}$ (appendix A).

##  6 Discussion and conclusion

As reviewed in the introduction, surface geostrophic velocities over the Gulf
Stream, the Kuroshio, and the Southern Ocean are primarily induced by surface
buoyancy anomalies in winter (Lapeyre, 2009; Isern-Fontanet and Hascoët, 2014;
González-Haro and Isern-Fontanet, 2014; Qiu et al., 2016; Miracca-Lage et al.,
2022). However, the kinetic energy spectra found in observations and numerical
models are too steep to be consistent with uniformly stratified surface
quasigeostrophic theory (Blumen, 1978; Callies and Ferrari, 2013). By
generalizing surface quasigeostrophic theory to account for variable
stratification, we have shown that surface buoyancy anomalies can generate a
variety of dynamical regimes depending on the stratification’s vertical
structure. Buoyancy anomalies generate longer range velocity fields over
decreasing stratification [$\sigma^{\prime}(z)\leq 0$] and shorter range
velocity fields over increasing stratification [$\sigma^{\prime}(z)\geq 0$].
As a result, the surface kinetic energy spectrum is steeper over decreasing
stratification than over increasing stratification. An exception occurs if
there is a large difference between the surface stratification and the deep
ocean stratification (as in the exponential stratified example of section 4).
In this case, we find regions of approximately homogenized surface buoyancy
surrounded by kinetic energy ribbons (similar to Arbic and Flierl, 2003) and
this spatial reorganization of the flow results in a steep kinetic energy
spectrum. By applying the variable stratification theory to the wintertime
North Atlantic and assuming that mixed-layer instability acts as a narrowband
small-scale surface buoyancy forcing, we find that the theory predicts a
surface kinetic energy spectrum between $k^{-4/3}$ and $k^{-7/3}$, which is
consistent with the observed wintertime $k^{-2}$ spectrum (Sasaki et al.,
2014; Callies et al., 2015; Vergara et al., 2019). There remains the problem
that mixed-layer instability may not be localized at a certain horizontal
scale but is forcing the surface flow at a wide range of scales (Khatri et
al., 2021). In this case we suggest that the main consequence of this
broadband forcing is again to flatten the $k^{-7/3}$ spectrum.

Over the summertime North Atlantic, buoyancy anomalies generate a more local
velocity field and the surface kinetic energy spectrum is flatter than in
winter. This contradicts the $k^{-3}$ spectrum found in observations and
numerical models (Sasaki et al., 2014; Callies et al., 2015). However,
observations also suggest that the surface geostrophic velocity is no longer
dominated by the surface buoyancy induced contribution, suggesting the
importance of interior potential vorticity for the summertime surface velocity
(González-Haro and Isern-Fontanet, 2014; Miracca-Lage et al., 2022). As such,
the surface kinetic energy predictions of the present model, which neglects
interior potential vorticity, are not valid over the summertime North
Atlantic.

The situation in the North Pacific is broadly similar to that in the North
Atlantic. In the Southern Ocean, however, the weak depth-averaged
stratification leads to values of $L_{H}$ close to 150-200 km. As such, the
bottom boundary becomes important at smaller horizontal scales than in the
North Atlantic. Regardless of whether the appropriate bottom boundary
condition is no-slip (56) or free-slip (57), in both cases, the resulting
inversion function implies a steepening to a $k^{-3}$ surface kinetic energy
spectrum (appendix A). The importance of the bottom boundary in the Southern
Ocean may explain the observed steepness of the surface kinetic energy spectra
[between $k^{-2.5}$ to $k^{-3}$ (Vergara et al., 2019)] even though the
surface geostrophic velocity seems to be largely due to surface buoyancy
anomalies throughout the year (González-Haro and Isern-Fontanet, 2014).

The claims made in this article can be explicitly tested in a realistic high-
resolution ocean model; this can be done by finding regions where the surface
streamfunction as reconstructed from sea surface height is highly correlated
to the surface streamfunction as reconstructed from sea surface buoyancy (or
temperature, as in González-Haro and Isern-Fontanet, 2014). Then, in regions
where both streamfunctions are highly correlated, the theory predicts that the
inversion function, as computed from the stratification [equation (15)],
should be identical to the inversion function computed through the surface
streamfunction and buoyancy fields [equations (1) and (16)]. Moreover, in
these regions, the model surface kinetic energy spectrum must be between the
inverse cascade and forward cascade kinetic energy spectra [equations (36) and
(37)].

Finally the vertical structure equation (9) along with the inversion relation
(13) between $\hat{\theta}_{\bm{k}}$ and $\hat{\psi}_{\bm{k}}$ suggest the
possibility of measuring the buoyancy frequency’s vertical structure, $N(z)$,
using satellites observations. This approach, however, is limited to regions
where the surface geostrophic velocity is largely due to surface buoyancy
anomalies. By combining satellite measurements of sea surface temperature and
sea surface height, we can use the inversion relation (13) to solve for the
inversion function. Then we obtain $N(z)$ by solving the inverse problem for
the vertical structure equation (9). How practical this approach is to
measuring the buoyancy frequency’s vertical structure remains to be seen.

###### Acknowledgements.

We sincerely thank Bob Hallberg, Isaac Held, and Sonya Legg for useful
comments on earlier drafts of this manuscript. We also thank Guillaume Lapeyre
and an anonymous reviewer whose comments and recommendations greatly helped
with our presentation. This report was prepared by Houssam Yassin under award
NA18OAR4320123 from the National Oceanic and Atmospheric Administration, U.S.
Department of Commerce. The statements, findings, conclusions, and
recommendations are those of the authors and do not necessarily reflect the
views of the National Oceanic and Atmospheric Administration, or the U.S.
Department of Commerce. The ECCO data (ECCO Consortium et al., 2021a, b) is
available from NASA’s Physical Oceanography Distributed Active Archive Center
(https://podaac.jpl.nasa.gov/ECCO). The pyqg model is available on Github
(https://github.com/pyqg/pyqg). [A] The small- and large-scale limits

###  .1 The small-scale limit

Let $h$ be a characteristic vertical length scale associated with $\sigma(z)$
near $z=0$. Then, in the small-scale limit, $k\,\sigma_{0}\,h\gg 1$, the
infinite bottom boundary condition (11) is appropriate. With the substitution

$\displaystyle\Psi(z)=\sigma(z)\,P(z),$ (59)

we transform the vertical structure equation (9) into a Schrödinger equation,

$\frac{\mathrm{d}^{2}P}{\mathrm{d}z^{2}}=\left[-\frac{1}{\sigma}\frac{\mathrm{d}^{2}\sigma}{\mathrm{d}z^{2}}+2\left(\frac{1}{\sigma}\frac{\mathrm{d}\sigma}{\mathrm{d}z}\right)^{2}+k^{2}\,\sigma^{2}\right]P,$
(60)

with a lower boundary condition

$\sigma\,P\rightarrow 0\quad\text{ as }\quad z\rightarrow-\infty.$ (61)

In the limit $k\,\sigma_{0}\,h\gg 1$, the solution to the Schrödinger equation
equation (60) is given by

$\Psi_{k}(z)\approx\sqrt{\frac{\sigma(z)}{\sigma_{0}}}\,\exp\left({k\,\int_{0}^{z}\sigma(s)\mathrm{d}s}\right).$
(62)

On substituting $\Psi_{k}(z)$ into the definition of the inversion function
(15), we obtain $m(k)\approx k/\sigma_{0}$ to leading order in
$(k\sigma_{0}h)^{-1}$. Therefore, the inversion relation in the small-scale
limit coincides with the familiar inversion relation of constant
stratification surface quasigeostrophic theory (Blumen, 1978; Held et al.,
1995).

###  .2 The large-scale free-slip limit

Let $k_{H}=2\pi/L_{H}$, where the horizontal length scale $L_{H}$ is defined
in equation (55). Then, in the large-scale limit, $k/k_{H}\ll 1$, we assume a
solution of the form

$\Psi_{k}(z)=\Psi_{k}^{(0)}(z)+\left(\frac{k}{k_{H}}\right)^{2}\,\Psi_{k}^{(1)}(z)+\cdots.$
(63)

Substituting the series expansion (63) into the vertical structure equation
(9) and applying the free-slip bottom boundary condition (57) yields

$\Psi_{k}(z)\approx
A\left[1+k^{2}\int_{-H}^{z}\sigma^{2}(s)\,\left(s+H\right)\mathrm{d}\,s+\cdots\right],$
(64)

where $A$ is a constant determined by the upper boundary condition (10). To
leading order in $k/k_{H}$, the large-scale vertical structure is independent
of depth.

Substituting the solution (64) into the definition of the inversion function
(15) gives

$m(k)\approx H\,k^{2}.$ (65)

Therefore, over a free-slip bottom boundary, the large-scale dynamics resemble
two-dimensional vorticity dynamics, generalizing the result of Tulloch and
Smith (2006) to arbitrary stratification $\sigma(z)$.

###  .3 The large-scale no-slip limit

Substituting the expansion (63) into the vertical structure equation (9) and
applying the no-slip lower boundary condition (56) yields

$\displaystyle\begin{split}\Psi_{k}&(z)\approx
B\Bigg{[}\int_{-H}^{z}\,\sigma^{2}(s)\,\mathrm{d}\,s\,+\\\
&k^{2}\,\int_{-H}^{z}\sigma^{2}(s_{3})\int_{-H}^{s_{3}}\int_{-H}^{s_{2}}\sigma^{2}(s_{1})\,\mathrm{d}s_{1}\,\mathrm{d}s_{2}\,\mathrm{d}s_{3}\Bigg{]},\end{split}$
(66)

where $B$ is a constant determined by the upper boundary condition (10).
Substituting the solution (66) into the definition of the inversion function
(15) gives

$m(k)\approx m_{1}\left(k_{\sigma}^{2}+k^{2}\right),$ (67)

where $k_{\sigma}=\sqrt{m_{0}/m_{1}}$ is analogous to the deformation
wavenumber, the constant $m_{0}$ is given by

$m_{0}=\left[\int_{-H}^{0}\sigma^{2}(s)\mathrm{d}s\right]^{-1},$ (68)

and $m_{1}$ is some constant determined by integrals of $\sigma(z)$. If
$\sigma(z)$ is positive then both $m_{0}$ and $m_{1}$ are also positive.
Therefore, over a no-slip bottom boundary, the large-scale dynamics resemble
those of the equivalent barotropic model.

[B]

Inversion function for piecewise constant stratification

We seek a solution to the vertical structure equation (9) for the piecewise
constant stratification (39) with upper boundary condition (10) and the
infinite lower boundary condition (11). The solution has the form

$\Psi_{k}(z)=\cosh\left(\sigma_{0}\,k\,z\right)+a_{2}\sinh\left(\sigma_{0}\,k\,z\right)$
(69)

for $-h\leq z\leq 0$, and

$\Psi_{k}(z)=a_{3}\,e^{\sigma_{\mathrm{pyc}}k(z+h)}$ (70)

for $-\infty<z<-h$. To determine $a_{2}$ and $a_{3}$, we require $\Psi_{k}(z)$
to be continuous across $z=-h$ and that its derivative satisfy

$\frac{1}{\sigma_{0}^{2}}\,\frac{\mathrm{d}\Psi_{k}(-h^{+})}{\mathrm{d}z}=\frac{1}{\sigma_{\mathrm{pyc}}^{2}}\,\frac{\mathrm{d}\Psi_{k}(-h^{-})}{\mathrm{d}z},$
(71)

where the $-$ and $+$ superscripts indicate limits from the below and above
respectively. Solving for $a_{2}$ and substituting the vertical structure (69)
into the definition of the inversion function (15) then yields $m(k)$.

[C]

The numerical model

We solve the time-evolution equation (6) using the pseudo-spectral pyqg model
(Abernathey et al., 2019). To take stratification into account, we use the
inversion relation (14). Given a stratification profile $\sigma(z)$ from
ECCOv4, we first interpolate the ECCOv4 stratification profile with a cubic
spline onto a vertical grid with 350 vertical grid points. We then numerically
solve the vertical structure equation (9), along with boundary conditions (10)
and either (56) or (57), and obtain the vertical structure at each wavevector
$\bm{k}$. Using the definition of the inversion function (15) then gives
$m(k)$.

We apply a large-scale forcing, $F$, between the (non-dimensional) wavenumbers
$3.5<k<4.5$ in all our simulations, corresponding to horizontal length scales
88 - 114 km. Otherwise, the forcing $F$ is as in Smith et al. (2002). The
dissipation term can be written as

$D=r_{d}\,\theta+\mathrm{ssd}$ (72)

where $r_{d}$ is a damping rate and $\mathrm{ssd}$ is small-scale dissipation.
Small-scale dissipation is through an exponential surface potential enstrophy
filter as in Arbic and Flierl (2003). Simulations have a horizontal resolution
of $1024^{2}$ horizontal grid points.

## References

  * Abernathey et al. (2019) Abernathey, R., and Coauthors, 2019: pyqg/pyqg: v0.3.0. Zenodo, https://doi.org/10.5281/zenodo.3551326.
  * Arbic and Flierl (2003) Arbic, B. K., and G. R. Flierl, 2003: Coherent vortices and kinetic energy ribbons in asymptotic, quasi two-dimensional f-plane turbulence. Physics of Fluids, 15, 2177–2189, https://doi.org/10.1063/1.1582183.
  * Arbic and Flierl (2004) Arbic, B. K., and G. R. Flierl, 2004: Baroclinically Unstable Geostrophic Turbulence in the Limits of Strong and Weak Bottom Ekman Friction: Application to Midocean Eddies. J. Phys. Oceanogr., 34, 2257–2273, https://doi.org/10.1175/1520-0485(2004)034¡2257:BUGTIT¿2.0.CO;2.
  * Blumen (1978) Blumen, W., 1978: Uniform Potential Vorticity Flow: Part I. Theory of Wave Interactions and Two-Dimensional Turbulence. J. Atmos. Sci., 35, 774–783, https://doi.org/10.1175/1520-0469(1978)035¡0774:UPVFPI¿2.0.CO;2.
  * Boccaletti et al. (2007) Boccaletti, G., R. Ferrari, and B. Fox-Kemper, 2007: Mixed Layer Instabilities and Restratification. J. Phys. Oceanogr., 37, 2228–2250, https://doi.org/10.1175/JPO3101.1.
  * Bretherton (1966) Bretherton, F. P., 1966: Critical layer instability in baroclinic flows. Quart. J. Roy. Meteor. Soc., 92, 325–334, https://doi.org/10.1002/qj.49709239302.
  * Burgess et al. (2015) Burgess, B. H., R. K. Scott, and T. G. Shepherd, 2015: Kraichnan–Leith–Batchelor similarity theory and two-dimensional inverse cascades. J. Fluid Mech., 767, 467–496, https://doi.org/10.1017/jfm.2015.26.
  * Callies and Ferrari (2013) Callies, J., and R. Ferrari, 2013: Interpreting Energy and Tracer Spectra of Upper-Ocean Turbulence in the Submesoscale Range (1–200 km). J. Phys. Oceanogr., 43, 2456–2474, https://doi.org/10.1175/JPO-D-13-063.1.
  * Callies et al. (2015) Callies, J., R. Ferrari, J. M. Klymak, and J. Gula, 2015: Seasonality in submesoscale turbulence. Nature, 6, https://doi.org/10.1038/ncomms7862.
  * Callies et al. (2016) Callies, J., G. Flierl, R. Ferrari, and B. Fox-Kemper, 2016: The role of mixed-layer instabilities in submesoscale turbulence. J. Fluid Mech., 788, 5–41, https://doi.org/10.1017/jfm.2015.700.
  * Charney (1971) Charney, J. G., 1971: Geostrophic Turbulence. J. Atmos. Sci., 28, 1087–1095, https://doi.org/10.1175/1520-0469(1971)028¡1087:GT¿2.0.CO;2.
  * ECCO Consortium et al. (2021a) ECCO Consortium, I. Fukumori, O. Wang, I. Fenty, G. Forget, P. Heimbach, and R. M. Ponte, 2021a: Ecco central estimate (version 4 release 4). Date accessed 28 January 2022, https://ecco.jpl.nasa.gov/drive/files.
  * ECCO Consortium et al. (2021b) ECCO Consortium, I. Fukumori, O. Wang, I. Fenty, G. Forget, P. Heimbach, and R. M. Ponte, 2021b: Synopsis of the ecco central production global ocean and sea-ice state estimate (version 4 release 4). https://doi.org/10.5281/zenodo.4533349.
  * Forget et al. (2015) Forget, G., J.-M. Campin, P. Heimbach, C. N. Hill, R. M. Ponte, and C. Wunsch, 2015: ECCO version 4: an integrated framework for non-linear inverse modeling and global ocean state estimation. Geoscientific Model Development, 8, 3071–3104, https://doi.org/10.5194/gmd-8-3071-2015.
  * Foussard et al. (2017) Foussard, A., S. Berti, X. Perrot, and G. Lapeyre, 2017: Relative dispersion in generalized two-dimensional turbulence. J. Fluid Mech., 821, 358–383, https://doi.org/10.1017/jfm.2017.253.
  * Gkioulekas and Tung (2007) Gkioulekas, E., and K. K. Tung, 2007: A new proof on net upscale energy cascade in two-dimensional and quasi-geostrophic turbulence. J. Fluid Mech., 576, 173–189, https://doi.org/10.1017/S0022112006003934.
  * González-Haro and Isern-Fontanet (2014) González-Haro, C., and J. Isern-Fontanet, 2014: Global ocean current reconstruction from altimetric and microwave SST measurements. J. Geophys. Res.: Oceans, 119, 3378–3391, https://doi.org/10.1002/2013JC009728.
  * González-Haro et al. (2020) González-Haro, C., J. Isern-Fontanet, P. Tandeo, and R. Garello, 2020: Ocean Surface Currents Reconstruction: Spectral Characterization of the Transfer Function Between SST and SSH. J. Geophys. Res.: Oceans, 125, https://doi.org/10.1029/2019JC015958.
  * Held et al. (1995) Held, I. M., R. T. Pierrehumbert, S. T. Garner, and K. L. Swanson, 1995: Surface quasi-geostrophic dynamics. J. Fluid Mech., 282, 1–20, https://doi.org/10.1017/S0022112095000012.
  * Isern-Fontanet and Hascoët (2014) Isern-Fontanet, J., and E. Hascoët, 2014: Diagnosis of high-resolution upper ocean dynamics from noisy sea surface temperatures. J. Geophys. Res.: Oceans, 119, 121–132, https://doi.org/10.1002/2013JC009176.
  * Isern-Fontanet et al. (2014) Isern-Fontanet, J., M. Shinde, and C. González-Haro, 2014: On the Transfer Function between Surface Fields and the Geostrophic Stream Function in the Mediterranean Sea. J. Phys. Oceanogr., 44, 1406–1423, https://doi.org/10.1175/JPO-D-13-0186.1.
  * Isern‐Fontanet et al. (2006) Isern‐Fontanet, J., B. Chapron, G. Lapeyre, and P. Klein, 2006: Potential use of microwave sea surface temperatures for the estimation of ocean currents. Geophys. Res. Lett., 33, https://doi.org/10.1029/2006GL027801.
  * Isern‐Fontanet et al. (2008) Isern‐Fontanet, J., G. Lapeyre, P. Klein, B. Chapron, and M. W. Hecht, 2008: Three-dimensional reconstruction of oceanic mesoscale currents from surface information. J. Geophys. Res.: Oceans, 113, https://doi.org/10.1029/2007JC004692.
  * Iwayama and Watanabe (2010) Iwayama, T., and T. Watanabe, 2010: Green’s function for a generalized two-dimensional fluid. Physical Review E, 82, https://doi.org/10.1103/PhysRevE.82.036307.
  * Khatri et al. (2021) Khatri, H., S. M. Griffies, T. Uchida, H. Wang, and D. Menemenlis, 2021: Role of mixed-layer instabilities in the seasonal evolution of eddy kinetic energy spectra in a global submesoscale permitting simulation. Geophys. Res. Lett., 48, e2021GL094777, https://doi.org/10.1029/2021GL094777.
  * Kraichnan (1967) Kraichnan, R. H., 1967: Inertial Ranges in Two‐Dimensional Turbulence. Physics of Fluids, 10, 1417–1423, https://doi.org/10.1063/1.1762301.
  * Kraichnan (1971) Kraichnan, R. H., 1971: Inertial-range transfer in two and three‐dimensional turbulence. J. Fluid Mech., 47, 525–535, https://doi.org/10.1017/S0022112071001216.
  * LaCasce (2017) LaCasce, J. H., 2017: The Prevalence of Oceanic Surface Modes. Geophys. Res. Lett., 44, 11 097–11 105, https://doi.org/10.1002/2017GL075430.
  * LaCasce and Mahadevan (2006) LaCasce, J. H., and A. Mahadevan, 2006: Estimating subsurface horizontal and vertical velocities from sea-surface temperature. J. Mar. Res., 64, 695–721, https://doi.org/10.1357/002224006779367267.
  * Lapeyre (2009) Lapeyre, G., 2009: What Vertical Mode Does the Altimeter Reflect? On the Decomposition in Baroclinic Modes and on a Surface-Trapped Mode. J. Phys. Oceanogr., 39, 2857–2874, https://doi.org/10.1175/2009JPO3968.1.
  * Lapeyre (2017) Lapeyre, G., 2017: Surface Quasi-Geostrophy. Fluids, 2, 7, https://doi.org/10.3390/fluids2010007.
  * Lapeyre and Klein (2006) Lapeyre, G., and P. Klein, 2006: Dynamics of the Upper Oceanic Layers in Terms of Surface Quasigeostrophy Theory. J. Phys. Oceanogr., 36, 165–176, https://doi.org/10.1175/JPO2840.1.
  * Larichev and McWilliams (1991) Larichev, V. D., and J. C. McWilliams, 1991: Weakly decaying turbulence in an equivalent‐barotropic fluid. Physics of Fluids A: Fluid Dynamics, 3, 938–950, https://doi.org/10.1063/1.857970.
  * Lilly (1989) Lilly, D. K., 1989: Two-Dimensional Turbulence Generated by Energy Sources at Two Scales. J. Atmos. Sci., 46, 2026–2030, https://doi.org/10.1175/1520-0469(1989)046¡2026:TDTGBE¿2.0.CO;2.
  * Maltrud and Vallis (1991) Maltrud, M. E., and G. K. Vallis, 1991: Energy spectra and coherent structures in forced two-dimensional and beta-plane turbulence. J. Fluid Mech., 228, 321–342, https://doi.org/10.1017/S0022112091002720.
  * Mensa et al. (2013) Mensa, J. A., Z. Garraffo, A. Griffa, T. M. Özgökmen, A. Haza, and M. Veneziani, 2013: Seasonality of the submesoscale dynamics in the Gulf Stream region. Ocean Dynamics, 63, 923–941, https://doi.org/10.1007/s10236-013-0633-1.
  * Miracca-Lage et al. (2022) Miracca-Lage, M., C. González-Haro, D. C. Napolitano, J. Isern-Fontanet, and P. S. Polito, 2022: Can the Surface Quasi-Geostrophic (SQG) Theory Explain Upper Ocean Dynamics in the South Atlantic? J. Geophys. Res.: Oceans, 127, e2021JC018 001, https://doi.org/10.1029/2021JC018001.
  * Pierrehumbert et al. (1994) Pierrehumbert, R. T., I. M. Held, and K. L. Swanson, 1994: Spectra of local and nonlocal two-dimensional turbulence. Chaos, Solitons & Fractals, 4, 1111–1116, https://doi.org/10.1016/0960-0779(94)90140-6.
  * Polvani et al. (1989) Polvani, L. M., N. J. Zabusky, and G. R. Flierl, 1989: Two-layer geostrophic vortex dynamics. Part 1. Upper-layer V-states and merger. J. Fluid Mech., 205, 215–242, https://doi.org/10.1017/S0022112089002016.
  * Qiu et al. (2020) Qiu, B., S. Chen, P. Klein, H. Torres, J. Wang, L.-L. Fu, and D. Menemenlis, 2020: Reconstructing Upper-Ocean Vertical Velocity Field from Sea Surface Height in the Presence of Unbalanced Motion. J. Phys. Oceanogr., 50, 55–79, https://doi.org/10.1175/JPO-D-19-0172.1.
  * Qiu et al. (2016) Qiu, B., S. Chen, P. Klein, C. Ubelmann, L.-L. Fu, and H. Sasaki, 2016: Reconstructability of Three-Dimensional Upper-Ocean Circulation from SWOT Sea Surface Height Measurements. J. Phys. Oceanogr., 46, 947–963, https://doi.org/10.1175/JPO-D-15-0188.1.
  * Sasaki et al. (2014) Sasaki, H., P. Klein, B. Qiu, and Y. Sasai, 2014: Impact of oceanic-scale interactions on the seasonal modulation of ocean dynamics by the atmosphere. Nature Communications, 5, https://doi.org/10.1038/ncomms6636.
  * Schorghofer (2000) Schorghofer, N., 2000: Energy spectra of steady two-dimensional turbulent flows. Physical Review E, 61, 6572–6577, https://doi.org/10.1103/PhysRevE.61.6572.
  * Smith et al. (2002) Smith, K. S., G. Boccaletti, C. C. Henning, I. Marinov, C. Y. Tam, I. M. Held, and G. K. Vallis, 2002: Turbulent diffusion in the geostrophic inverse cascade. J. Fluid Mech., 469, 13–48, https://doi.org/10.1017/S0022112002001763.
  * Stammer (1997) Stammer, D., 1997: Global Characteristics of Ocean Variability Estimated from Regional TOPEX/POSEIDON Altimeter Measurements. J. Phys. Oceanogr., 27, 1743–1769, https://doi.org/10.1175/1520-0485(1997)027¡1743:GCOOVE¿2.0.CO;2.
  * Sukhatme and Smith (2009) Sukhatme, J., and L. M. Smith, 2009: Local and nonlocal dispersive turbulence. Phys. Fluids, 21, 056603, https://doi.org/10.1063/1.3141499.
  * Tulloch and Smith (2006) Tulloch, R., and K. S. Smith, 2006: A theory for the atmospheric energy spectrum: Depth-limited temperature anomalies at the tropopause. Proc. Natl. Acad. Sci. (USA), 103, 14 690–14 694, https://doi.org/10.1073/pnas.0605494103.
  * Vergara et al. (2019) Vergara, O., R. Morrow, I. Pujol, G. Dibarboure, and C. Ubelmann, 2019: Revised Global Wave Number Spectra From Recent Altimeter Observations. J. Geophys. Res.: Oceans, 124, 3523–3537, https://doi.org/10.1029/2018JC014844.
  * Watanabe and Iwayama (2004) Watanabe, T., and T. Iwayama, 2004: Unified Scaling Theory for Local and Non-local Transfers in Generalized Two-dimensional Turbulence. Journal of the Physical Society of Japan, 73, 3319–3330, https://doi.org/10.1143/JPSJ.73.3319.
  * Wunsch (1997) Wunsch, C., 1997: The Vertical Partition of Oceanic Horizontal Kinetic Energy. J. Phys. Oceanogr., 27, 1770–1794, https://doi.org/10.1175/1520-0485(1997)027¡1770:TVPOOH¿2.0.CO;2.