# The buoyancy staircase limit in surface quasigeostrophic turbulence

Houssam Yassin<EMAIL_ADDRESS>Program in Atmospheric and Oceanic
Sciences, Princeton University, Princeton, NJ 08544, USA

###### Abstract

Surface buoyancy gradients over a quasigeostrophic fluid permit the existence
of surface-trapped Rossby waves. The interplay of these Rossby waves with
surface quasigeostrophic turbulence results in latitudinally inhomogeneous
mixing that, under certain conditions, culminates in a surface buoyancy
staircase: a meridional buoyancy profile consisting of mixed-zones punctuated
by sharp buoyancy gradients, with eastward jets centred at the sharp gradients
and weaker westward flows in between. In this article, we investigate the
emergence of this buoyancy staircase limit in surface quasigeostrophic
turbulence and we examine the dependence of the resulting dynamics on the
vertical stratification. Over decreasing stratification
[$\mathrm{d}N(z)/\mathrm{d}z\leq 0$, where $N(z)$ is the buoyancy frequency],
we obtain flows with a longer interaction range (than in uniform
stratification) and highly dispersive Rossby waves. In the staircase limit, we
find straight jets that are perturbed by eastward propagating along jet waves,
similar to two-dimensional barotropic $\beta$-plane turbulence. In contrast,
over increasing stratification [$\mathrm{d}N(z)/\mathrm{d}z\geq 0$], we obtain
flows with shorter interaction range and weakly dispersive Rossby waves. In
the staircase limit, we find sinuous jets with large latitudinal meanders
whose shape evolves in time due to the westward propagation of weakly
dispersive along jet waves. These along jet waves have larger amplitudes over
increasing stratification than over decreasing stratification, and, as a
result, the ratio of domain-averaged zonal to meridional speeds is two to
three times smaller over increasing stratification than over decreasing
stratification. Finally, we find that, for a given Rhines wavenumber, jets
over increasing stratification are closer together than jets over decreasing
stratification.

###### keywords:

##  1 Introduction

Perturbations to a barotropic (i.e., depth-independent) fluid with a
background potential vorticity gradient, $\beta>0$, propagate westward as
Rossby waves. In a turbulent flow, the non-linear interplay between Rossby
waves and turbulence results in the latitudinally inhomogeneous mixing of
potential vorticity, which, through a positive dynamical feedback,
spontaneously reorganizes the flow into one characterized by eastward jets
(Dritschel & McIntyre, 2008). The ultimate limit of such inhomogeneous mixing,
which can be achieved for sufficiently large values of $\beta$, is a potential
vorticity staircase: a piecewise constant potential vorticity profile
consisting well-mixed regions separated by isolated discontinuities, with
eastward jets centred at the discontinuities and westward flows in between
(Danilov & Gurarie, 2004; Dunkerton & Scott, 2008; Scott & Dritschel, 2012,
2019).

Analogously, a buoyancy gradient at the surface of a quasigeostrophic fluid
supports the existence of surface-trapped Rossby waves that are less
dispersive than their barotropic counterparts (Held et al., 1995; Lapeyre,
2017). The purpose of this article is to investigate the formation of zonal
jets in the presence of a background surface buoyancy gradient and to examine
the realizability of surface buoyancy staircases in the surface
quasigeostrophic model. Although the present study is the first to
systematically investigate the emergence of surface quasigeostrophic jets,
there are previous studies which make use of the uniformly stratified surface
quasigeostrophic model with a background buoyancy gradient. These include
Smith et al. (2002), who derive the dependence of the diffusion coefficient of
a passive tracer in the presence a background buoyancy gradient. Another is
Sukhatme & Smith (2009), who, in their investigation of $\alpha$-turbulence
models with a background gradient, note that, because of the decreased
interaction range, surface quasigeostrophic jets in uniform stratification
should be narrower than their counterparts in the barotropic model. Finally,
Lapeyre (2017) demonstrates that jets can indeed form in the uniformly
stratified surface quasigeostrophic model.

We also investigate how surface quasigeostrophic jets depend on the underlying
vertical stratification. Yassin & Griffies (2022) show that the vertical
stratification modifies the interaction range of vortices in the surface
quasigeostrophic model. Suppose we have an infinitely deep fluid governed by
the time-evolution of geostrophic buoyancy anomalies at its upper boundary.
Then if the stratification is decreasing [$N^{\prime}(z)\leq 0$, where $N(z)$
is buoyancy frequency] towards the fluid’s surface (that is, the upper
boundary), then the interaction range is longer than in the uniformly
stratified model and the resulting turbulence is characterized by thin
buoyancy filaments — analogous to the thin vorticity filaments in two-
dimensional barotropic turbulence. Conversely, if the stratification is
increasing [$N^{\prime}(z)\geq 0$] towards the surface, then the interaction
range is shorter than in uniform stratification, and the buoyancy field
appears spatially diffuse and lacks thin filamentary structures. In this
article, we find that the interaction range is related to Rossby wave
dispersion: flows with a longer interaction range have more dispersive Rossby
waves whereas flows with a shorter interaction range have less dispersive
Rossby waves. One of our aims is to characterize the dependence of surface
quasigeostrophic jets on the functional form of the vertical stratification.

There are two motivations behind the present work. The first is its potential
relevance to the upper ocean. Buoyancy anomalies at the ocean’s surface are
governed by the surface quasigeostrophic model (Lapeyre & Klein, 2006; LaCasce
& Mahadevan, 2006; Isern‐Fontanet et al., 2006). Both numerical
(Isern‐Fontanet et al., 2008; Lapeyre, 2009; Qiu et al., 2016, 2020; Miracca-
Lage et al., 2022) as well as observational (González-Haro & Isern-Fontanet,
2014) studies indicate that a significant fraction of the surface geostrophic
velocity is induced by sea surface buoyancy anomalies, especially over
wintertime extratropical currents. Moreover, upper ocean turbulence has been
found to be anisotropic (Maximenko et al., 2005; Scott et al., 2008), with
significant differences in anisotropy between major extratropical currents and
other regions in the ocean (Wang et al., 2019). However, our neglect of the
planetary $\beta$ effect, as well our assumption of vanishing interior
potential vorticity, may limit the direct relevance of this study to the upper
ocean.

The second motivation is that the variable stratification surface
quasigeostrophic model is a simple two-dimensional model in which we can
investigate how jet dynamics depend on the stratification’s vertical
structure. Another such model is the equivalent barotropic model for which the
deformation radius represents the rigidity of the free surface. Small values
of the deformation radius lead to a pliable free surface allowing a
significant degree of horizontal divergence. The resulting flow then has an
exponentially short interaction range, with a horizontal attenuation on the
order of the deformation radius (Polvani et al., 1989), and with approximately
non-dispersive Rossby waves. Consequently, for a finite deformation radius, we
obtain jets whose width is on the order of the deformation radius with a fixed
meandering shape (Scott et al., 2022). In contrast, for the variable
stratification surface quasigeostrophic model, rather than just specifying a
constant (i.e., the deformation wavenumber), one instead has to specify the
stratification’s functional form, $N(z)$. Over decreasing stratification
$[N^{\prime}(z)<0]$, because of the longer interaction range and the more
dispersive waves, we obtain jets similar to the two-dimensional barotropic
model. Conversely, over increasing stratification $[N^{\prime}(z)>0]$, the
shorter interaction range along with the weakly dispersive waves lead to
sinuous jets whose shape evolves in time through the propagation of weakly
dispersive along jet waves. Moreover, because of these along jet waves, a
smaller fraction of the total energy is contained in the zonal mode over
increasing stratification (with a shorter interaction range) than over
decreasing stratification (with a longer interaction range).

The remainder of this article is organized as follows. Section 2 introduces
the variable stratification surface quasigeostrophic model and shows how the
stratification’s vertical structure controls both the interaction range of
point vortices as well as the dispersion of surface-trapped Rossby waves.
Then, in section 3, we introduce two wavenumbers, $k_{\varepsilon}$ and
$k_{r}$, whose ratio, $k_{\varepsilon}/k_{r}$, forms the key non-dimensional
parameter of this study; here, $k_{\varepsilon}$ is a wavenumber depending on
the energy injection rate whereas $k_{r}$ is a wavenumber depending on surface
damping rate. This non-dimensional number is a generalization of the non-
dimensional number used in previous studies (Danilov & Gurarie, 2002;
Sukoriansky et al., 2007; Scott & Dritschel, 2012). By considering an
idealized buoyancy staircase, we also investigate how the Rhines wavenumber
relates to the jet spacing under decreasing, increasing, and uniform
stratification. Section 4 then presents numerical experiments detailing the
emergence of the staircase limit as $k_{\varepsilon}/k_{r}$ is increased for
various stratification profiles. In addition, we also present experiments
where we fix the external parameters and vary the vertical stratification
alone. Finally, we conclude in section 5.

##  2 The interaction range and wave dispersion

###  2.1 Equations of motion

Consider an infinitely deep fluid with zero interior potential vorticity. The
geostrophic streamfunction, $\psi$, then satisfies

$\frac{\partial}{\partial
z}\left(\frac{1}{\sigma^{2}}\frac{\partial\psi}{\partial
z}\right)+\nabla^{2}\psi=0$ (1)

in the fluid interior, $z\in(-\infty,0)$. The horizontal Laplacian is denoted
by $\nabla^{2}=\p_{x}^{2}+\p_{y}^{2}$ and the non-dimensional stratification
is given by

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

where $N(z)$ is the buoyancy frequency and $f$ is the constant local value of
the Coriolis parameter. Time-evolution is determined by the material
conservation of surface potential vorticity (Bretherton, 1966),

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

at the upper boundary, $z=0$, where $\sigma_{0}=\sigma(0)$. Explicitly, the
time-evolution equation is

$\frac{\partial\theta}{\partial
t}+J(\psi,\theta)+\Lambda\,\partial_{x}\theta=F-D,$ (4)

at $z=0$, where
$J(\psi,\theta)=\partial_{x}\psi\,\partial_{y}\theta-\partial_{x}\theta\,\partial_{y}\,\psi$
represents the advection of $\theta$ by the geostrophic velocity,
$\boldsymbol{u}=\hat{\boldsymbol{z}}\times\boldsymbol{\nabla}\psi$. The
frequency, $\Lambda$, is given by

$\Lambda=\frac{1}{\sigma_{0}^{2}}\frac{\mathrm{d}U}{\mathrm{d}z}\Big{|}_{z=0},$
(5)

where $U(z)$ is a background zonal geostrophic flow. Without loss of
generality, we have assumed that $U(0)=0$ in the time-evolution equation (4)
to eliminate a constant advective term. The dissipation, $D$, consists of
linear damping and small-scale dissipation,

$D=r\,\theta+\mathrm{ssd},$ (6)

where $r$ is the damping rate. The forcing, $F$, and the small-scale
dissipation, $\mathrm{ssd}$, are described in section 4.

The surface buoyancy anomaly, $b|_{z=0}$, is related to the surface potential
vorticity, $\theta$, through

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

Therefore, the time-evolution equation (4) equivalently states that surface
buoyancy anomalies are materially conserved in the absence of forcing and
dissipation. In addition, the frequency, $\Lambda$, corresponds to a
meridional buoyancy gradient,

$\frac{\mathrm{d}B}{\mathrm{d}y}\Big{|}_{z=0}=-f\,\sigma_{0}^{2}\,\Lambda,$
(8)

where $B(y,z)$ is the buoyancy field that is in geostrophic balance with
background zonal velocity, $U(z)$.

If we further assume a doubly periodic domain in the horizontal, then we can
expand the streamfunction as

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

where $\boldsymbol{x}=(x,y)$ is the horizontal position vector, $z$ is the
vertical coordinate, $\boldsymbol{k}=(k_{x},k_{y})$ is the horizontal
wavevector, $k=\left\lvert\boldsymbol{k}\right\rvert$ is the horizontal
wavenumber, and $t$ is the time coordinate. The non-dimensional wavenumber-
dependent vertical structure, $\Psi_{k}(z)$, is determined by the boundary
value problem (Yassin & Griffies, 2022)

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

with the upper boundary condition

$\Psi_{k}(0)=1,$ (11)

and lower boundary condition

$\Psi_{k}\rightarrow 0\quad\textrm{as }\quad z\rightarrow-\infty.$ (12)

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

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

The corresponding Fourier expansion of the surface potential vorticity is
given by

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

where

$\hat{\theta}_{\boldsymbol{k}}=-m(k)\,\hat{\psi}_{\boldsymbol{k}},$ (15)

and the function $m(k)$ is given by

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

The function $m(k)$ relates $\hat{\theta}_{\boldsymbol{k}}$ to
$\hat{\psi}_{\boldsymbol{k}}$ in the Fourier space inversion relation (15) and
so we call $m(k)$ the _inversion function_.

To recover the well-known case of the uniformly stratified quasigeostrophic
model (Held et al., 1995), set $\sigma(z)=\sigma_{0}$. Then the vertical
structure equation (10) along with boundary conditions (11) and (12) yield the
exponentially decaying vertical structure
$\Psi_{k}(z)=\exp\left({\sigma_{0}\,k\,z}\right)$. On substituting
$\Psi_{k}(z)$ into equation (16), we obtain a linear inversion function

$m(k)=\frac{k}{\sigma_{0}}$ (17)

and hence [from the inversion relation (15)] a linear-in-wavenumber inversion
relation
$\hat{\theta}_{\boldsymbol{k}}=-(k/\sigma_{0})\,\hat{\psi}_{\boldsymbol{k}}$.

###  2.2 The inversion function and spatial locality

Figure 1: The inversion functions, $m(k)$ [in panel (a)] for two
stratification profiles [panel (b)] given by the piecewise stratification
profile (18). One stratification profile is increasing
[$\sigma^{\prime}(z)\geq 0$, blue], with $\sigma_{0}=1$,
$\sigma_{\mathrm{pyc}}=0.15$, $h_{\mathrm{mix}}=0.01$, and
$h_{\mathrm{lin}}=0.05$. The other stratification profile is decreasing
[$\sigma^{\prime}(z)\leq 0$, red] with $\sigma_{0}=1$,
$\sigma_{\mathrm{pyc}}=10$, $h_{\mathrm{mix}}=0.01$, and
$h_{\mathrm{lin}}=0.05$. The thin black line is given by $k/\sigma_{0}$ where
$\sigma_{0}=1$, whereas the blue and red lines are given by
$k/\sigma_{\mathrm{pyc}}$ with $\sigma_{\mathrm{pyc}}=0.15$ for the thin blue
line and $\sigma_{\mathrm{pyc}}=10$ for the thin red line.

The inversion function $m(k)$, which is determined by the stratification’s
vertical structure, controls the spatial locality of the resulting turbulence.
We illustrate this point with the following piecewise stratification profile,

$\sigma(z)=\begin{cases}\sigma_{0}\quad&\text{for}\quad-
h_{\mathrm{mix}}<z<0\\\
\sigma_{0}+\Delta\sigma\left(\frac{z+h_{\mathrm{mix}}}{h_{\mathrm{lin}}}\right)\quad&\text{for}\quad-(h_{\mathrm{mix}}+h_{\mathrm{lin}})<z<-h_{\mathrm{mix}}\\\
\sigma_{\mathrm{pyc}}\quad&\text{for}\quad-\infty<z<-(h_{\mathrm{mix}}+h_{\mathrm{lin}}),\end{cases}$
(18)

where $\Delta\sigma=\sigma_{0}-\sigma_{\mathrm{pyc}}$. At small horizontal
scales, where $k\gg k_{s}$, and

$k_{s}=1/\left(\sigma_{0}\,h_{\mathrm{mix}}\right),$ (19)

then $m(k)\approx k/\sigma_{0}$, as in the uniformly stratified model of Held
et al. (1995). Likewise, in the large-scale limit, where $k\ll
k_{\mathrm{pyc}}$, and

$k_{\mathrm{pyc}}=\begin{cases}1/\left(\sigma_{\mathrm{pyc}}\,h_{\mathrm{mix}}\right)\quad&\text{for}\quad\Delta\sigma\leq
0\\\
\sigma_{\mathrm{pyc}}/\left(\sigma_{0}^{2}\,h_{\mathrm{mix}}\right)\quad&\text{for}\quad\Delta\sigma>0,\end{cases}$
(20)

then $m(k)\approx k/\sigma_{\mathrm{pyc}}$. However, for wavenumbers between
$k_{\mathrm{pyc}}\lesssim k\lesssim k_{s}$, the inversion function takes an
approximate power law form

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

where $m_{0}>0$ and $\alpha\geq 0$. The power $\alpha$ depends on the ratio
$\sigma_{\mathrm{pyc}}/\sigma_{0}$ between the deep and surface
stratification. If the stratification decreases towards the surface
[$\sigma^{\prime}(z)\leq 0$, or $\sigma_{\mathrm{pyc}}/\sigma_{0}>1$] then
$\alpha>1$, with $\sigma_{\mathrm{pyc}}/\sigma_{0}\rightarrow\infty$ sending
$\alpha\rightarrow 2$. In contrast, if the stratification increases towards
the surface [$\sigma^{\prime}(z)\geq 0$, or
$\sigma_{\mathrm{pyc}}/\sigma_{0}<1$] then $\alpha<1$, with
$\sigma_{\mathrm{pyc}}/\sigma_{0}\rightarrow 0$ sending $\alpha\rightarrow 0$.
Thus, for wavenumbers $k_{\mathrm{pyc}}\lesssim k\lesssim k_{s}$, the
inversion relation (15) has the approximate form

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

where $\hat{\xi}_{\boldsymbol{k}}=\hat{\theta}_{\boldsymbol{k}}/m_{0}$, which
is the inversion relation for $\alpha$-turbulence (Pierrehumbert et al., 1994;
Smith et al., 2002; Sukhatme & Smith, 2009). Figure 1 provides two examples,
one with decreasing stratification (with $\alpha\approx 1.50$) and another
with increasing stratification (with $\alpha\approx 0.49$).

To see how the parameter $\alpha$ modifies the resulting dynamics, consider a
point vortex at the origin, given by
$\xi=\delta(\left\lvert\boldsymbol{x}\right\rvert)$, where
$\left\lvert\boldsymbol{x}\right\rvert$ is the horizontal distance from the
vortex centre, and $\delta(\left\lvert\boldsymbol{x}\right\rvert)$ is the
Dirac delta. If $\alpha=2$, then the streamfunction induced by the point
vortex is logarithmic,
$\psi(\left\lvert\boldsymbol{x}\right\rvert)=\log(\left\lvert\boldsymbol{x}\right\rvert)/(2\pi)$.
If $0<\alpha<2$, then
$\psi(\left\lvert\boldsymbol{x}\right\rvert)=-C_{\alpha}/\left\lvert\boldsymbol{x}\right\rvert^{2-\alpha}$
where $C_{\alpha}>0$ is a constant (Iwayama & Watanabe, 2010). Smaller
$\alpha$ leads to vortices with velocities decaying more quickly with the
horizontal distance $\left\lvert\boldsymbol{x}\right\rvert$, and hence a
shorter interaction range. Thus, the vertical stratification modifies the
relationship between a surface buoyancy anomaly and its induced velocity
field: a surface buoyancy anomaly over decreasing stratification
[$\sigma^{\prime}(z)\leq 0$] generates a longer range velocity field than an
identical buoyancy anomaly over increasing stratification
[$\sigma^{\prime}(z)\geq 0$].

###  2.3 Wave dispersion in variable stratification

The background gradient term, $\Lambda$, in the time-evolution equation (4)
allows for the propagation of surface-trapped Rossby waves. Substituting a
wave solution of the form
$\psi(x,z,t)=\Psi_{k}(z)\,\exp\left[{\mathrm{i}\left(\boldsymbol{k}\cdot\boldsymbol{r}-\omega
t\right)}\right]$, where the vertical structure $\Psi_{k}(z)$ satisfies the
boundary value problem (10)–(12), into the time-evolution equation (4) yields
the angular frequency

$\omega(\boldsymbol{k})=-\frac{\Lambda\,k_{x}}{m(k)}.$ (23)

Given the relationship (8) between the meridional surface buoyancy gradient
$\mathrm{d}B/\mathrm{d}y|_{z=0}$ and the frequency $\Lambda$, a poleward
decreasing buoyancy gradient ($f\mathrm{d}{B}/\mathrm{d}{y}<0$) implies
westward propagating $(\omega<0)$ Rossby waves.

The dispersion relation (23) shows that Rossby wave dispersion is coupled to
the flow’s interaction range and hence the stratification’s vertical
structure. If we approximate the inversion function as a power law (21)
between $k_{\mathrm{pyc}}\lesssim k\lesssim k_{s}$, then the zonal phase
speed, $c=\omega/k_{x}$, becomes $c\sim 1/k^{\alpha}$. Therefore, at these
horizontal scales, Rossby waves are more dispersive over decreasing
stratification (with $\alpha>1$) than over increasing stratification (with
$\alpha<1$). In the limit that $\sigma_{0}\gg\sigma_{\mathrm{pyc}}$ in which
$\alpha\rightarrow 0$, then $c\approx$ constant, and so Rossby waves become
non-dispersive.

##  3 From edge waves to surface-trapped jets

The emergence of jets in barotropic $\beta$-plane turbulence is due to two
properties of the potential vorticity (Dritschel & McIntyre, 2008; Scott &
Dritschel, 2019). The first is the resilience of strong latitudinal potential
vorticity gradients to mixing (i.e., ”Rossby wave elasticity”, Dritschel &
McIntyre, 2008). Regions with weak latitudinal potential vorticity gradients
are preferentially mixed, weakening the gradient in these regions and
enhancing the gradient in regions where the latitudinal potential vorticity
gradient is already strong (Dritschel & Scott, 2011). The ultimate limit of
such latitudinally inhomogeneous mixing is a potential vorticity staircase
(Danilov & Gryanik, 2004; Dritschel & McIntyre, 2008; Scott & Dritschel,
2012), which consists of uniform regions of potential vorticity punctuated by
sharp potential vorticity gradients. The second property is that, through
potential vorticity inversion, strong (positive) latitudinal gradients in
potential vorticity correspond to eastward jets. Therefore, inverting a
potential vorticity staircase produces a flow with eastward zonal jets centred
at the sharp frontal zones, with weaker westward flows in between (Scott &
Dritschel, 2019).

However, the limit of a potential vorticity staircase is only achieved for
sufficiently large values of the non-dimensional number
$k_{\varepsilon}/k_{\mathrm{Rh}}$ (Scott & Dritschel, 2012), which is a ratio
of the forcing intensity wavenumber, $k_{\varepsilon}$, to the Rhines
wavenumber, $k_{\mathrm{Rh}}$. The forcing intensity wavenumber is given by
(Maltrud & Vallis, 1991)

$k_{\varepsilon}=(\beta^{3}/\varepsilon_{\mathcal{K}})^{1/5},$ (24)

where $\varepsilon_{\mathcal{K}}$ is the kinetic energy injection rate in the
barotropic model, and is obtained by setting the turbulent strain rate equal
to the Rossby wave frequency (Vallis & Maltrud, 1993). The Rhines wavenumber
is given by (Rhines, 1975)

$k_{\mathrm{Rh}}=\sqrt{\beta/U_{\mathrm{rms}}},$ (25)

where $U_{\mathrm{rms}}$ is the rms velocity. Scott & Dritschel (2012) found
that the ratio $k_{\varepsilon}/k_{\mathrm{Rh}}$ controls the structure of
zonal jets in barotropic $\beta$-plane turbulence; as
$k_{\varepsilon}/k_{\mathrm{Rh}}$ is increased, the zonal jet strength
increases and the potential vorticity gradient at the jet core becomes larger,
with the staircase limit approached as $k_{\varepsilon}/k_{\mathrm{Rh}}\sim
O(10)$.

Jet formation in surface quasigeostrophic turbulence proceeds similarly, with
the surface buoyancy (which is proportional to $\theta$) taking the role of
the potential vorticity and the frequency, $\Lambda$, taking the role of the
potential vorticity gradient, $\beta$. In this section, we first derive a non-
dimensional number analogous to $k_{\varepsilon}/k_{\mathrm{Rh}}$ for surface
quasigeostrophy. Then we consider how vertical stratification (and the non-
locality parameter $\alpha$) modifies jet structure in the buoyancy staircase
limit, as well as how it modifies the relationship between the Rhines
wavenumber and the jet spacing.

Before proceeding, we comment on two differences between two-dimensional
barotropic turbulence and its surface quasigeostrophic counterpart. First, in
the absence of forcing and dissipation, the kinetic energy,

$\mathcal{K}=-\frac{1}{2}\,\overline{\psi\nabla^{2}\psi}=\frac{1}{2}\,\overline{\left\lvert
u\right\rvert^{2}},$ (26)

is a conserved constant in two-dimensional barotropic turbulence (the overline
denotes an area average). With a constant kinetic energy injection rate,
$\varepsilon_{\mathcal{K}}$, and a linear damping rate, $r$, the equilibrium
kinetic energy is $\mathcal{K}=\varepsilon_{\mathcal{K}}/2r$. By definition,
the rms velocity is given by $U_{\mathrm{rms}}=\sqrt{2\mathcal{K}}$. Combining
this expression with the definition of the kinetic energy (26) and
substituting into the definition of the Rhines wavenumber (25) yields a Rhines
wavenumber expressed in terms of external parameters alone,

$k_{\mathrm{Rh}}=\beta^{1/2}(r/\varepsilon_{\mathcal{K}})^{1/4}.$ (27)

In contrast, in surface quasigeostrophy, the total energy,

$\mathcal{E}=-\frac{1}{2}\overline{\psi|_{z=0}\,\theta},$ (28)

is a conserved constant in the absence of forcing and dissipation and there is
no general relationship between the rms velocity, $U_{\mathrm{rms}}$, and the
equilibrium total energy, $\mathcal{E}=\varepsilon/2r$, where $\varepsilon$ is
the total energy injection rate in the surface quasigeostrophic model.
Therefore, we are not generally able to express the Rhines wavenumber in terms
of the external parameters $\varepsilon$, $\Lambda$, and $r$. Second, because
$\mathcal{E}$ and $\mathcal{K}$ have different dimensions, the kinetic energy
injection in the barotropic model, $\varepsilon_{\mathcal{K}}$, has different
dimensions than the total energy injection rate in the surface
quasigeostrophic model, $\varepsilon$. In particular, $\varepsilon$ has
dimensions of $L^{2}/T^{3}$.

###  3.1 The forcing intensity wavenumber

To obtain the forcing intensity wavenumber, $k_{\varepsilon}$, we compare the
Rossby wave frequency (23) to the turbulent strain rate, $\omega_{s}(k)$. If
the inversion function is not approximately constant (i.e., $\alpha\neq 0$)
then the strain rate is (Yassin & Griffies, 2022)

$\omega_{s}(k)\sim\varepsilon^{1/3}\,k^{4/3}\,\left[m(k)\right]^{-1/3}.$ (29)

In particular, if $m(k)=m_{0}\,k^{\alpha}$, then $\omega_{s}(k)\sim
m_{0}^{1/3}\varepsilon^{1/3}\,k^{\left(4-\alpha\right)/3}$. Setting the
absolute value of the Rossby wave frequency for waves with $k=k_{x}$ equal to
the turbulent strain rate (29) yields the condition

$k_{\varepsilon}\left[m(k_{\varepsilon})\right]^{2}\sim\frac{\left\lvert\Lambda\right\rvert^{3}}{\varepsilon}.$
(30)

A solution to this equation always exists because $\mathrm{d}m/\mathrm{d}k\geq
0$. If the inversion function takes the power law form (21), then we obtain

$k_{\varepsilon}=\left(\frac{\left\lvert\Lambda\right\rvert^{3}}{m_{0}^{2}\,\varepsilon}\right)^{1/\left(2\alpha+1\right)},$
(31)

which is equivalent to a wavenumber derived in Smith et al. (2002).

###  3.2 The damping rate wavenumber and the Rhines wavenumber

Suppose the inversion function takes an approximate power law form,
$m(k)\approx m_{0}\,k^{\alpha}$, near the energy containing wavenumbers. Then
the generalization of the Rhines wavenumber at these wavenumbers is

$k_{\mathrm{Rh}}=\left(\frac{\Lambda}{m_{0}\,U_{\mathrm{rms}}}\right)^{1/\alpha}.$
(32)

However, unlike in two-dimensional barotropic turbulence where
$U_{\mathrm{rms}}=\sqrt{2\mathcal{K}}=\sqrt{\varepsilon_{\mathcal{K}}/r}$, we
do not have a general relationship between $U_{\mathrm{rms}}$ and the external
parameters $r$ and $\varepsilon$ in surface quasigeostrophic turbulence. To
obtain a second wavenumber that depends on the damping rate, $r$, we follow
Smith et al. (2002). From dimensional considerations, the energy spectrum at
small wavenumbers is

$E_{\Lambda}(k)\sim\Lambda^{2}\,k^{-(\alpha+3)}/m_{0}.$ (33)

Then, defining $k_{r}$ as the wavenumber at which the inverse cascade halts,
we obtain

$\frac{\varepsilon}{2r}\approx\int_{k_{r}}^{\infty}E(k)\,\mathrm{d}k\approx\left(\frac{\Lambda^{2}/m_{0}}{\alpha+2}\right)k_{r}^{-(\alpha+2)},$
(34)

where the second equality follows because the integral is dominated by its
peak at low wavenumbers. Solving for $k_{r}$ and neglecting any non-
dimensional coefficients, we obtain

$k_{r}=\left(\frac{\Lambda^{2}\,r}{m_{0}\,\varepsilon}\right)^{1/(\alpha+2)}.$
(35)

Note that the damping rate wavenumber, $k_{r}$, has the same dependence on
$\Lambda$, $\varepsilon$, and $r$ as the Rhines wavenumber, $k_{\mathrm{Rh}}$,
only if $\alpha=2$.

###  3.3 Surface potential vorticity inversion

A perfect surface potential vorticity staircase consists of mixed zones of
halfwidth $b$, where $\mathrm{d}\theta/\mathrm{d}y=-\Lambda$, separated by
jump discontinuities at which $\mathrm{d}\theta/\mathrm{d}y=\infty$. We find
it more conveniant to work with the relative surface potential vorticity,
$\theta$, rather than the total surface potential vorticity,
$\theta+\Lambda\,y$. In this case, if the total surface potential vorticity,
$\theta+\Lambda\,y$, is a perfect staircase with step width $2b$, then the
relative surface potential vorticity, $\theta$, is a $2b$-periodic sawtooth
wave.

Our first question is whether such a staircase is possible for general $m(k)$.
To answer this question, we consider the velocity field induced by a jump
discontinuity in $\theta$. For a jump discontinuity in an infinite domain,

$\theta=\begin{cases}\Delta\theta\quad&\text{for}\quad 0<y<\infty\\\
0\quad&\text{for}\quad-\infty<y<0,\end{cases}$ (36)

the zonal velocity is given by

$u=\frac{\Delta\theta}{2\pi}\int_{-\infty}^{\infty}\frac{\mathrm{e}^{\mathrm{i}\,k_{y}y}}{m\left(\left\lvert
k_{y}\right\rvert\right)}\mathrm{d}k_{y}.$ (37)

If $m(k)=m_{0}\,k^{\alpha}$, then this expression is proportional to
$\left\lvert y\right\rvert^{\alpha-1}$ if $\alpha\neq 1$ and logarithmic
otherwise, and so the zonal velocity diverges at $y=0$ if $\alpha\leq 1$.
Consequently, we expect that a perfect staircase should not be possible over
constant or increasing stratification due to the divergence of the zonal
velocity at a jump discontinuity.

Figure 2: Panel (a) shows a sloping sawtooth function (thick black line) along
with its derivative (thin black line). Panel (b) shows the normalized zonal
velocity induced by the sloping sawtooth function in panel (a) for various
values of the parameter $\alpha$. Panel (c) shows the normalized zonal
velocity induced by the sawtooth function in (a) in the increasing (blue line)
and decreasing stratifications (red line) shown in figure 1.

We therefore consider the more general case of a sloping staircase, where
there is a finite frontal zone of width $2a$ between the mixed zones. In this
case, $\theta$ is a $2(a+b)$-periodic sloping sawtooth wave (see figure 2),
and is given by the periodic extension of

$\theta=\Lambda\begin{cases}-\left[y-{(a+b)}\right]\quad&\text{for}\quad{a}<y<{a+b}\\\
\frac{\,b}{a}y\quad&\text{for}\quad\left\lvert y\right\rvert\leq{a}\\\
-\left[y+{(a+b)}\right]\quad&\text{for}\quad-{(a+b)}<y<-{a}.\end{cases}$ (38)

The meridional gradient $\mathrm{d}\theta/\mathrm{d}y$ is then a piecewise
constant $2(a+b)$-periodic function

$\frac{\mathrm{d}\theta}{\mathrm{d}y}=\Lambda\begin{cases}-1\quad&\text{for}\quad{a}<y<{a+b}\\\
\frac{\,b}{a}\quad&\text{for}\quad\left\lvert y\right\rvert\leq{a}\\\
-1\quad&\text{for}\quad-{(a+b)}<y<-{a}.\end{cases}$ (39)

Therefore the gradient in the frontal zones exceeds the gradient in the mixed
zones by a factor of $b/a$, which approaches infinity as
$b/a\rightarrow\infty$ in the sawtooth wave limit.

The zonal velocity, $u=-\partial_{y}\psi$, is obtained by using the inversion
relation (15) to solve for the streamfunction. Alternatively, taking the
meridional derivative of surface potential vorticity (3) gives

$\frac{\partial\theta}{\partial y}=\frac{1}{\sigma_{0}^{2}}\frac{\partial
u}{\partial z}\Big{|}_{z=0}.$ (40)

Then in Fourier space [$\partial_{y}\rightarrow\mathrm{i}k_{y}$ and
$\sigma_{0}^{-2}\partial_{z}|_{z=0}\rightarrow m(k)$] we obtain

$\hat{u}_{\boldsymbol{k}}=\frac{1}{m(k)}\left(\mathrm{i}\,k_{y}\,\hat{\theta}_{\boldsymbol{k}}\right),$
(41)

which shows that the induced zonal velocity is obtained by smoothing
$\mathrm{d}\theta/\mathrm{d}y$ by the function $m(k)$. An immediate
consequence is that the east-west asymmetry in the zonal velocity is
fundamentally due to the east-west asymmetry in the gradient
$\mathrm{d}\theta/\mathrm{d}y$.

Figure 2 shows an example of sloping sawtooth $\theta$ profile along with the
induced zonal velocities. For a power law inversion function,
$m(k)=m_{0}k^{\alpha}$, the parameter $\alpha$ modifies the zonal velocity in
two ways. First, in more local flows (with smaller $\alpha$), the zonal
velocity decays more rapidly away from the jet centre, as expected. Second,
the degree of smoothing increases with $\alpha$, and so more local regimes
(with smaller $\alpha$) are more east-west asymmetric, with the ratio
$\left\lvert u_{\mathrm{min}}\right\rvert/u_{\mathrm{max}}$ taking smaller
values for smaller $\alpha$. Figure 3(b) shows $\left\lvert
u_{\mathrm{min}}\right\rvert/u_{\mathrm{max}}$ as a function of $a/b$ for
$\alpha\in\\{1/2,\,1,\,3/2,\,2\\}$. For $\alpha=2$, we obtain $\left\lvert
u_{\mathrm{min}}\right\rvert/u_{\mathrm{max}}\rightarrow 1/2$ in the limit
$a/b\rightarrow 0$ so that eastward jets are only twice as strong as westward
flows in the perfect staircase limit (Danilov & Gurarie, 2004; Dritschel &
McIntyre, 2008). At $\alpha=3/2$, we find $\left\lvert
u_{\mathrm{min}}\right\rvert/u_{\mathrm{max}}\approx 0.29$ in the
$a/b\rightarrow 0$ limit so that eastward jets are now more than three time as
strong as westward flows. Once $\alpha\leq 1$, then the maximum jet velocity
diverges as $\alpha\rightarrow 0$ [figure 3(a)] and so $\left\lvert
u_{\mathrm{min}}\right\rvert/u_{\mathrm{max}}\rightarrow 0$ as $a/b\rightarrow
0$.

If $m(k)$ is not a power law, then the results are similar so long as $m(k)$
can be approximated by a power law at small wavenumbers. Figure 2 shows the
induced velocity for the inversion functions computed from idealized
stratifications profiles (shown in figure 1). Because these inversion
functions can be approximated by power laws $m(k)\approx k^{0.49}$ and
$m(k)\approx k^{1.50}$ at small wavenumbers, the induced velocity fields
nearly coincide with the velocity fields computed from power law inversion
functions with $\alpha=0.5$ and $\alpha=1.5$.

Figure 3: Properties of zonal velocity profiles induced by sloping sawtooth
profiles (38) of $\theta$ as a function of the non-dimensional frontal zone
width $a/b$ separating the mixed zones for four values of $\alpha$. Panel (a)
shows the maximum zonal velocity, panel (b) shows the ratio of westward speed
to eastward speed, panel (c) shows the rms zonal velocity, and panel (d) shows
the product $L_{j}k_{\mathrm{Rh}}$ where $L_{j}=a+b$ is the halfwidth
separation (the distance between $U_{\mathrm{mix}}$ and $U_{\mathrm{min}}$)
and $k_{\mathrm{Rh}}$ is the Rhines wavenumber (32).

Finally, we examine how the Rhines wavenumber, $k_{\mathrm{Rh}}$, relates to
jet spacing. Let

$L_{j}=a+b$ (42)

be the half-separation between the jets, i.e., the half distance between
consecutive zonal velocity maxima. For two-dimensional barotropic turbulence
(i.e., the $\alpha=2$ case), we have $L_{j}=45^{1/4}/k_{\mathrm{Rh}}\approx
2.59/k_{\mathrm{Rh}}$ in the staircase limit (i.e, for $a/b\rightarrow 0$,
Dritschel & McIntyre, 2008; Scott & Dritschel, 2012). This result is found by
solving for the zonal velocity induced by a staircase with halfwidth
$L_{j}=b$, taking the rms of the zonal velocity, and then substituting into
the definition of the generalized Rhines wavenumber (32). As figure 3(d)
shows, because the velocity field induced by a perfect staircase depends on
the inversion function, $m(k)$, the relationship between $L_{j}$ and
$k_{\mathrm{Rh}}$ also depends on the inversion function. For $m(k)=k^{3/2}$,
an analogous calculation gives $L_{j}\approx 2.35/k_{\mathrm{Rh}}$ in the
staircase limit. For $\alpha=1$, even though the maximum velocity diverges at
$a/b\rightarrow 0$, the rms velocity asymptotes to a constant value, and so we
obtain a half jet-separation of $L_{j}\approx 1.73/k_{\mathrm{Rh}}$ (figure
3). Finally in the $\alpha=1/2$ case, although the rms speed has not converged
by $a/b=10^{-6}$, the product $L_{j}\,k_{\mathrm{Rh}}$ is approaching values
close to zero.

##  4 Numerical Simulations

###  4.1 The numerical model

We use the pyqg pseudo-spectral model (Abernathey et al., 2019) which solves
the time-evolution equation (4) in a square domain with side length $L=2\pi$.
Time-stepping is through a third-order Adam-Bashforth scheme with small-scale
dissipation achieved through a scale-selective exponential filter (Smith et
al., 2002; Arbic & Flierl, 2003),

$\mathrm{ssd}=\begin{cases}1\quad\text{for}\quad k\leq k_{0}\\\
e^{-a(k-k_{0})^{4}}\quad\text{for}\quad k>k_{0},\end{cases}$ (43)

with $a=23.6$ and $k_{0}=0.65k_{\mathrm{Nyq}}$ where $k_{\mathrm{Nyq}}=\pi$ is
the Nyquist wavenumber. The forcing is isotropic, centred at wavenumber
$k_{f}=80$, and normalized so that the energy injection rate is
$\varepsilon=1$ (see appendix B in Smith et al., 2002). However, the effective
energy injection rate, $\varepsilon_{\mathrm{eff}}$, is smaller than
$\varepsilon$ due to dissipation. To determine $\varepsilon_{\mathrm{eff}}$
from numerical simulations, we use
$\varepsilon_{\mathrm{eff}}=2\,r\,\mathcal{E}$ where $\mathcal{E}$ is the
equilibrated total energy diagnosed from the model. In what follows, we report
values of $k_{\varepsilon}/k_{r}$ using $\varepsilon_{\mathrm{eff}}$ instead
of $\varepsilon$. The model is integrated forward in time until at least
$t=5/r$ to allow the fluid to reach equilibrium. All model runs use $1024^{2}$
horizontal grid points.

###  4.2 For what values of $k_{\varepsilon}/k_{r}$ do jets form?

Figure 4: Diagnostics from five series of simulations as a function of the
non-dimensional number $k_{\varepsilon}/k_{r}$. The first three series of
simulations have inversion function $m(k)=k^{\alpha}$ with
$\alpha\in\\{\nicefrac{{1}}{{2}},\,1,\,\nicefrac{{3}}{{2}}\\}$. For the other
two series, the inversion functions are shown in figure 8. Panel (a) shows the
ratio of energy in the zonal mode to total energy. Panel (b) shows the ratio
of domain averaged zonal speed to domain averaged meridional speed. Panel (c)
shows the ratio of westward zonal speed to eastward zonal speed. Panel (d)
shows the relationship between the halfwidth jet spacing, $L_{j}$, and the
Rhines wavenumber, $k_{\mathrm{Rh}}$.

For our first set of simulations, we vary $k_{\varepsilon}/k_{r}$ over the
values shown in figure 4. We do so by fixing $k_{r}=8$ and varying
$k_{\varepsilon}$. For a given value of $k_{\varepsilon}$, we choose $\Lambda$
and $r$ so as to maintain $k_{r}=8$ (the energy injection rate, $\varepsilon$,
is fixed at unity for all model runs). Given $k_{r}$ and $k_{\varepsilon}$, we
rearrange the definition of $k_{r}$ (35) to solve for $\gamma=r\,\Lambda^{2}$,

$\gamma=m_{0}\,\varepsilon\,k_{r}^{\alpha+2},$ (44)

then solve for $r$ in the implicit equation (30) for $k_{\varepsilon}$ ,

$r=\frac{\gamma}{\left(\varepsilon\,k_{\varepsilon}\,[m(k_{\varepsilon})]^{2}\right)^{2/3}},$
(45)

and finally use the definition $\gamma=r\,\Lambda^{2}$ to solve for $\Lambda$.

####  4.2.1 Power law inversion functions

We first describe the results from three series of simulations with power law
inversion functions, $m(k)=k^{\alpha}$, with
$\alpha\in\\{\nicefrac{{1}}{{2}},\,1,\,\nicefrac{{3}}{{2}}\\}$. Summary
diagnostics from these simulations are shown in figure 4. In panel (a), we
observe that the ratio of energy in the zonal mode to total energy,
$\mathcal{E}_{\mathrm{zonal}}/\mathcal{E}$, increases with
$k_{\varepsilon}/k_{r}$, and that the majority of the total energy is in the
zonal mode for sufficiently large $k_{\varepsilon}/k_{r}$. For a fixed
$k_{\varepsilon}/k_{r}$, more of the total energy is zonal in more non-local
flows (with larger $\alpha$) than in more local flows (with smaller $\alpha$);
for $\alpha=\nicefrac{{3}}{{2}}$, we have
$\mathcal{E}_{\mathrm{zonal}}/\mathcal{E}\approx 1$ by
$k_{\varepsilon}/k_{r}\approx 6$ as compared to $k_{\varepsilon}/k_{r}\approx
12$ for $\alpha=1$. Moreover, for $\alpha=\nicefrac{{1}}{{2}}$, we find that
$\mathcal{E}_{\mathrm{zonal}}/\mathcal{E}$ asymptotes to approximately 0.9
once $k_{\varepsilon}/k_{r}\approx 18$ with little subsequent change for
larger values of $k_{\varepsilon}/k_{r}$. In panel (b), we observe a striking
contrast in the ratio $\overline{\left\lvert
u\right\rvert}/\overline{\left\lvert v\right\rvert}$ between different values
of $\alpha$ (the overline denotes a domain average). For
$\alpha=\nicefrac{{3}}{{2}}$, the domain averaged zonal speed,
$\overline{\left\lvert u\right\rvert}$, is approximately eight times larger
than the domain averaged meridional speed, $\overline{\left\lvert
v\right\rvert}$, for large $k_{\varepsilon}/k_{r}$. In contrast, for
$\alpha=\nicefrac{{1}}{{2}}$, $\overline{\left\lvert u\right\rvert}$ only
exceeds $\overline{\left\lvert v\right\rvert}$ by a multiple of two for large
$k_{\varepsilon}/k_{r}$.

Figure 5: Snapshots of the relative surface potential vorticity, $\theta$, for
simulations with power law inversion functions, $m(k)=k^{\alpha}$. In each
snapshot, the $\theta$ field is normalized by its maximum value in the
snapshot. Only one quarter of the domain is shown (i.e., $512^{2}$ grid
points).

Figure 6: The zonal mean total surface potential vorticity,
$\theta+\Lambda\,y$, in black and the zonal mean zonal velocity, $U$, in grey.

Next, we examine the jet structure for different $\alpha$ as a function of
$k_{\varepsilon}/k_{r}$. Figure 5 shows $\theta$-snapshots from model runs
with $m(k)=k^{\alpha}$. For each value of $\alpha$, two model runs are shown:
one where jets have just become visible in the $\theta$-snapshot and another
with the largest value of $k_{\varepsilon}/k_{r}$, which we expect to be
closest to the staircase limit. The jets are visible in these snapshots as the
regions with strong gradients. Because these are $\theta$-snapshots rather
than $(\theta+\Lambda\,y)$-snapshots, the $(\theta+\Lambda\,y)$-staircase is
instead a $\theta$-sawtooth, and the mixed zones between the jets are
approximately linear in $\theta$. We confirm this to be the case in figure 6,
where the zonal averages of the total surface potential vorticity,
$\theta+\Lambda\,y$, and the zonal velocity are shown. For the
$\alpha=\nicefrac{{3}}{{2}}$ and $\alpha=1$ cases, we observe an approximate
staircase structure with nearly uniform mixed zones separated by frontal
zones, and with jets centred at sharp $\theta$ gradients. As expected from the
idealized staircases of section 3, close to the staircase limit, the
$\alpha=1$ jets are narrower than the $\alpha=\nicefrac{{3}}{{2}}$ jets, and
the ratio of maximum westward speed to maximum eastward speed,
$|U_{\mathrm{min}}|/|U_{\mathrm{max}}|$, is smaller at $\alpha=1$ than at
$\alpha=\nicefrac{{3}}{{2}}$.

In contrast to the $\alpha=\nicefrac{{3}}{{2}}$ and the $\alpha=1$ series, the
$\alpha=\nicefrac{{1}}{{2}}$ series approaches the staircase limit slowly with
$k_{\varepsilon}/k_{r}$. The $\alpha=\nicefrac{{1}}{{2}}$ staircase remains
smooth even at $k_{\varepsilon}/k_{r}=42$ [figure 6(c)]. The ratio of frontal
zone width to mixed zone width, $a/b$, is between $0.5$ and $0.65$ for
$\alpha=\nicefrac{{1}}{{2}}$ jets. In contrast, this ratio is between $0.15$
and $0.2$ for the $\alpha=\nicefrac{{3}}{{2}}$ and $\alpha=1$ jets. In part,
the broadness of the $\alpha=\nicefrac{{1}}{{2}}$ frontal zones is a
consequence of zonal averaging in the presence of large amplitude undulations.
However, it is evident from the $\theta$-snapshots of figure 5 that the
$\alpha=\nicefrac{{1}}{{2}}$ frontal zones are indeed broader than the
$\alpha=\nicefrac{{3}}{{2}}$ and $\alpha=1$ frontal zones [e.g., compare
panels (a) and (d) with (f) in figure 5], even without zonal averaging.

We now examine how the generalized Rhines wavenumber, $k_{\mathrm{Rh}}$,
relates to the jet spacing. From figure 3(d), a ratio of $a/b\approx 0.2$
leads to a $L_{j}\,k_{\mathrm{Rh}}\approx 2.2$ for
$\alpha=\nicefrac{{3}}{{2}}$ and $L_{j}\,k_{\mathrm{Rh}}\approx 2.0$ for
$\alpha=1$. But as figure 4(d) shows, we find values closer to
$L_{j}\,k_{\mathrm{Rh}}\approx 3$ for both of these cases. In contrast, for
the $\alpha=\nicefrac{{1}}{{2}}$ jets, figure 3(d) predicts $1.98\lesssim
L_{j}\,k_{\mathrm{Rh}}\lesssim 2.5$ for the observed range of $0.5\lesssim
a/b\lesssim 0.65$, but we find $L_{j}\,k_{\mathrm{Rh}}\approx 1.5$ for
$k_{\varepsilon}/k_{r}\geq 18$, which is smaller than predicted.

Returning to figure 5, we observe that there are undulations along the jets,
with smaller values of $\alpha$ corresponding to larger amplitude undulations.
These undulations propagate as waves and are less dispersive for smaller
$\alpha$, propagating eastward for $\alpha=\frac{3}{2}$, westward for
$\alpha=\nicefrac{{1}}{{2}}$, and are nearly stationary for $\alpha=1$.
Moreover, the waves in the $\alpha=\nicefrac{{1}}{{2}}$ case maintain their
shape as they propagate for a significant fraction of the domain, although
they eventually disperse or merge with other along jet waves. That we obtain
larger amplitude along jet undulations for smaller $\alpha$ is a consequence
of the more local inversion operator (15) at smaller $\alpha$. A jet in a
highly local flow (with small $\alpha$) is “a coherent structure that hangs
together strongly while being easy to push sideways” (McIntyre, 2008, in the
context of equivalent barotropic jets). However, although both an equivalent
barotropic jet and an $\alpha=\nicefrac{{1}}{{2}}$ jet exhibit large
meridional undulations, the undulations in the equivalent barotropic case are
frozen in place (because of a vanishing group velocity at large scales,
McIntyre, 2008) and so the equivalent barotropic jet behaves like a meandering
river with a fixed shape. In contrast, the $\alpha=\nicefrac{{1}}{{2}}$ jet
behaves like a flexible string whose shape evolves in time with the
propagation of weakly dispersive waves. Another difference between the two
cases is that an equivalent barotropic jet has a width given by the
deformation radius. In contrast, there is no analogous characteristic scale
for $\alpha=\nicefrac{{1}}{{2}}$ jets and, in principle, the jets should
become infinitely thin as $k_{\varepsilon}/k_{r}\rightarrow\infty$.

Figure 7: The total energy spectrum, $E(k)$, as a function of the wavenumber,
$k=k_{x}^{2}+k_{x}^{2}$, for three simulations with power law inversion
functions, $m(k)=k^{\alpha}$. The values of $k_{\varepsilon}/k_{r}$ are 18.0
for panel (a), 21.0 for panel (b), and 42.3 for panel (c).

Energy spectra for the three power law simulations are shown in figure 7. The
energy spectrum obtained from dimensional analysis (33) gives a
$k^{-\alpha-3}$ wavenumber dependence, which leads to the familiar $k^{-5}$
spectrum for beta-plane barotropic turbulence ($\alpha=2$). Although early
investigations (Chekhlov et al., 1996; Huang et al., 2000; Danilov & Gryanik,
2004) found a $k^{-5}$ spectrum in barotropic $\beta$-plane turbulence, Scott
& Dritschel (2012) instead found a shallower $k^{-4}$ spectrum in the
staircase limit (suggested earlier by Danilov & Gryanik, 2004; Danilov &
Gurarie, 2004), which they explained as a consequence of the sharp
discontinuities of the staircase. Generalizing their argument to the present
case, a one dimensional $\theta(y)$ series with discontinuities implies a
Fourier series with coefficients decaying as $k^{-1}$, leading to a
$\theta^{2}$ spectrum of $k^{-2}$, and hence an energy spectrum

$E(k)\sim k^{-2}\left[m(k)\right]^{-1}.$ (46)

If $m(k)\sim k^{\alpha}$, then we obtain a spectrum $E(k)\sim k^{-\alpha-2}$,
which yields the $k^{-4}$ spectrum observed in Scott & Dritschel (2012), where
$\alpha=2$. For $\alpha=\nicefrac{{3}}{{2}}$, $\alpha=1$, and
$\alpha=\nicefrac{{1}}{{2}}$, the predicted spectrum is proportional to
$k^{-3.5}$, $k^{-3}$, and $k^{-2.5}$, respectively. The diagnosed spectra
shown in figure 7 are consistent with these shallow spectra, instead of energy
spectrum (33) obtained from dimensional considerations.

####  4.2.2 Inversion functions from $\sigma(z)$

Figure 8: Inversion functions [panel (a) and (c)] along with their
corresponding stratification profiles [panels (b) and (d), respectively]. The
stratification profiles are given by the piecewise function (18). For panel
(a), we have $\sigma_{0}=1$, $\sigma_{\mathrm{pyc}}=0.1$,
$h_{\mathrm{mix}}=0.01$, and $h_{\mathrm{lin}}=0.05$. For panel (c), we have
$\sigma_{0}=0.133$, $\sigma_{\mathrm{pyc}}=1$, $h_{\mathrm{mix}}=0.125$, and
$h_{\mathrm{lin}}=0.2$. The thin grey lines in panels (a) and (c) are given by
$k/\sigma_{0}$ and $k/\sigma_{\mathrm{pyc}}$.

We also ran two series of simulations where we specified a piecewise
stratification profile (18), and then obtained $m(k)$ by solving the boundary
value problem (10)–(12) at each wavenumber. The stratification profiles and
the resulting inversion functions are shown in figure 8. One case consists of
an increasing stratification profile [$\sigma^{\prime}(z)\geq 0$] with
$\sigma_{0}=1$, $\sigma_{\mathrm{pyc}}=0.1$, $h_{\mathrm{mix}}=0.01$ and
$h_{\mathrm{lin}}=0.05$. The resulting $m(k)$ is approximately linear for
$k\gtrsim 70$ and transitions to an approximate sub-linear wavenumber
dependence $m(k)\sim k^{0.40}$ for wavenumbers $5\lesssim k\lesssim 50$. The
second case consists of a decreasing stratification profile
[$\sigma^{\prime}(z)\leq 0$] with $\sigma_{0}=0.13$,
$\sigma_{\mathrm{pyc}}=1$, $h_{\mathrm{mix}}=0.125$ and
$h_{\mathrm{lin}}=0.2$. The resulting $m(k)$ is approximately linear at
wavenumbers $k\gtrsim 60$ and transitions to an approximate super linear
wavenumber dependence $m(k)\sim k^{1.50}$ between $3\lesssim k\lesssim 60$.

As seen in figure 4, the $\sigma^{\prime}(z)\leq 0$ case is similar to the
$\alpha=\nicefrac{{3}}{{2}}$ case, with the various diagnostics close to the
$\alpha=\nicefrac{{3}}{{2}}$ counterpart. In contrast, there are significant
differences between the $\sigma^{\prime}(z)\geq 0$ simulations and the
$\alpha=\nicefrac{{1}}{{2}}$ simulations. In the $\sigma^{\prime}(z)\geq 0$
series, the ratio of energy in the zonal mode to total energy continues to
increase as $k_{\varepsilon}/k_{r}$ is increased, whereas it asymptotes to a
constant in the $\alpha=\nicefrac{{1}}{{2}}$ series. Moreover, the ratio of
domain average zonal speed to domain averaged meridional speed,
$\overline{|u|}/\overline{|v|}$, is generally larger in the $\sigma\geq 0$
series than in the $\alpha=\nicefrac{{1}}{{2}}$ series. Finally, for the
largest values of $k_{\varepsilon}/k_{r}$, the product
$L_{j}\,k_{\mathrm{Rh}}$ reaches smaller values in the $\sigma\geq 0$
simulations than in the $\alpha=\nicefrac{{1}}{{2}}$ simulations.

Figure 9: Snapshots of the relative surface potential vorticity, $\theta$,
normalized by its maximum value in the snapshot, for simulations with
inversion functions shown in figure 8. Only one quarter of the domain is shown
(i.e., $512^{2}$ grid points).

These differences can be explained by the snapshots of figure 9 as well as the
zonal averages of figure 6. As expected from the model diagnostics, both the
snapshots and the zonal average from the $\sigma^{\prime}\leq 0$ simulation
are qualitatively similar to the $\alpha=\nicefrac{{3}}{{2}}$ simulation. In
contrast, the $\sigma^{\prime}\geq 0$ snapshot is evidently closer to the
staircase limit than the $\alpha=\nicefrac{{1}}{{2}}$ snapshot: the mixed
zones are more homogeneous and the frontal zones are sharper. The zonal
average of the $\sigma^{\prime}\geq 0$ simulation in figure 6 also shows how
the $\sigma^{\prime}\geq 0$ simulation is closer to the staircase limit than
the $\alpha=\nicefrac{{1}}{{2}}$ simulation, although, again, zonal averaging
in the presence of large amplitude undulations is artificially smoothing the
jets. Therefore, the differences in the diagnostics between the
$\sigma^{\prime}\geq 0$ series and the $\alpha=\nicefrac{{1}}{{2}}$ series
stem from the more rapid approach (i.e., at smaller $k_{\varepsilon}/k_{r}$)
of the $\sigma^{\prime}\geq 0$ series to the staircase limit.

###  4.3 Simulations with fixed parameters

The dependence of the non-dimensional number $k_{\varepsilon}/k_{r}$ on the
external parameters $\varepsilon,\Lambda$, and $r$ depends on the functional
form of $m(k)$. For example, if $m(k)\sim k^{\alpha}$, then

$k_{\varepsilon}/{k_{r}}=|\Lambda|^{\frac{4-\alpha}{(2\alpha+1)(\alpha+2)}}\,\varepsilon^{\frac{\alpha-1}{(1+2\alpha)(\alpha+2)}}\,r^{\frac{-1}{\alpha+2}}.$
(47)

Because the forcing intensity wavenumber, $k_{\varepsilon}$, is obtained by
solving the implicit equation for $k_{\varepsilon}$ (30), an analogous
expression for $k_{\varepsilon}/k_{r}$ is not possible for general $m(k)$.
However, at sufficiently large $k_{\varepsilon}$, the inversion function
asymptotes to $m(k_{\varepsilon})\approx k_{\varepsilon}/\sigma_{0}$ and so,
using $\alpha$-turbulence expression for $k_{\varepsilon}$ (31) with
$\alpha=1$, we obtain

$k_{\varepsilon}/k_{r}\approx|\Lambda|^{\frac{\alpha}{\alpha+2}}\,\varepsilon^{\frac{1-\alpha}{3\alpha+6}}\,r^{\frac{-1}{\alpha+2}}\,m_{0}^{\frac{1}{\alpha+2}}\,\sigma_{0}^{2/3}$
(48)

for large $k_{\varepsilon}$, where $\alpha$ is the approximate power law
dependence of $m(k)$ near $k_{r}$.

Therefore, simulations with identical $k_{\varepsilon}/k_{r}$ but distinct
inversion functions cannot be directly compared because they have different
values of $\Lambda$ and $r$. Here, we investigate how the stratification
modifies jet structure as all other parameters are held fixed. We therefore
run two additional series of simulations with the stratification profiles and
inversion functions shown in figure 1. The stratification profiles were chosen
so that they both have identical stratification at the upper boundary. One
case corresponds to an increasing stratification profile, $\sigma^{\prime}\geq
0$, with an approximate power law dependence of $m(k)\sim k^{0.49}$ at small
wavenumbers. The second case consists of a decreasing stratification profile,
$\sigma^{\prime}\leq 0$, with a $m(k)\sim k^{1.50}$ at small wavenumbers.
Aside from the different stratification profiles, these two series of
simulations are run under the same conditions as the constant stratification
($\alpha=1$) simulations of section 4.2, with identical values of $\Lambda$,
$\varepsilon$, and $r$.

Figure 10: As in figure 4, but the $\sigma^{\prime}\geq 0$ and
$\sigma^{\prime}\leq 0$ series now only differ from the $\sigma^{\prime}=0$
(i.e., $\alpha=1$) series only in the vertical stratification (and hence the
inversion function).

Summary diagnostics are shown in figure 10. We see that, at a fixed value of
$\Lambda$ and $r$, more of the total energy is in the zonal mode in the
$\sigma^{\prime}(z)\leq 0$ simulation than in the constant stratification
simulation, which in turn is larger than the $\sigma^{\prime}(z)\geq 0$
simulation (and similarly for the ratio of area averaged zonal to meridional
speeds, $\overline{\left\lvert u\right\rvert}/\overline{\left\lvert
v\right\rvert}$). Therefore, increased non-locality (larger $\alpha$) promotes
anisotropy in the velocity field and leads to larger zonal velocities relative
to meridional velocities. Indeed, figure 11 shows $\theta$ snapshots from
these simulations; the more local, $\sigma^{\prime}\geq 0$, simulations have
larger meridional undulations along the jets. Moreover, compared to the
$k_{\varepsilon}/k_{r}=15$ constant stratification simulation in figure 5(c),
the $\sigma^{\prime}(z)\leq 0$ simulation in figure 11(a) is closer to the
staircase limit whereas the frontal zones in the $\sigma^{\prime}(z)\geq 0$
simulation [figure 11(c)] remain broad. Finally, we show values of the product
$L_{j}\,k_{\mathrm{Rh}}$, relating the Rhines wavenumber to the half spacing
between the jets, in figure 10(c). These values are similar to those in shown
in figure 4(c).

Figure 11: Snapshots of relative surface potential vorticity, $\theta$, where
$\theta$ is normalized by its maximum value in the snapshot. Panels (a) and
(c) are from simulations with identical $\Lambda$, $r$, and $\varepsilon$ as
the $\alpha=1$ simulation shown in figure 5(c), whereas panels (b) and (d) are
from simulations with identical $\Lambda$, $r$, and $\varepsilon$ as the
$\alpha=1$ simulation shown in figure 5(d). Only one quarter of the domain is
shown (i.e., $512^{2}$ grid points).

##  5 Conclusion

We have examined the emergence of staircase-like buoyancy structures in
surface quasigeostrophic turbulence with a mean background buoyancy gradient.
We found that the stratification’s vertical structure controls the locality of
the inversion operator and the dispersion of surface-trapped Rossby waves. As
we go from decreasing stratification profiles [$\sigma^{\prime}(z)\leq 0$] to
increasing stratification profiles [$\sigma^{\prime}(z)\geq 0$], the inversion
operator becomes more local and Rossby wave less dispersive. In all cases, we
find that the non-dimensional ratio, $k_{\varepsilon}/k_{r}$, controls the
extent of inhomogeneous buoyancy mixing. Larger $k_{\varepsilon}/k_{r}$
correspond to sharper buoyancy gradients at jet centres with larger peak jet
velocities that are separated by more homogeneous mixed-zones. Moreover, we
found that the staircase limit is reached at smaller $k_{\varepsilon}/k_{r}$
in more non-local flows; the staircase limit is reached by
$k_{\varepsilon}/k_{r}\approx 15$ for our $\sigma\leq 0$ simulations compared
to $k_{\varepsilon}/k_{r}\approx 25$ for our $\sigma\geq 0$ simulations.

In addition, once the staircase limit is reached, the dynamics of the jets
depends on the locality of the inversion operator and, hence, on the
stratification’s vertical structure. In flows with a more non-local inversion
operator [or decreasing stratification, $\sigma^{\prime}(z)\leq 0$], we obtain
straight jets that are perturbed by dispersive, eastward propagating, along
jet waves. In contrast, for more local flows [or over increasing
stratification, $\sigma^{\prime}(z)\geq 0$], we obtain jets with latitudinal
meanders on the order of the jet spacing. The shape of these jets evolves in
time as these meanders propagate westwards as weakly dispersive waves.

The inversion operator’s locality is also reflected in two more aspects of the
dynamics. First, the domain-averaged zonal speed exceeds the domain-averaged
meridional speed by approximately a factor of eight in our most non-local
simulations, whereas this ratio is merely two in our most local simulations.
This observation is consistent with the fact that jets are narrower and
exhibit larger latitudinal meanders in more local flows. Second, for a given
Rhines wavenumber, jets in more local flows are closer together. Indeed, we
found $L_{j}\,k_{\mathrm{Rh}}\approx 3-4$ in our most non-local simulations,
where $L_{j}$ is the jet half spacing, as compared to
$L_{j}\,k_{\mathrm{Rh}}\approx 0.5-1.5$ in our most local simulations.

Several open questions remain. First, we have not examined the dynamics of the
along jet waves. As we observed, these waves propagate eastwards in our most
non-local simulations [with $\sigma^{\prime}(z)\leq 0$] but westwards for our
most local simulations [with $\sigma^{\prime}(z)\geq 0$]. These waves are not
described by the dispersion relation (23); rather, the relevant model is that
of freely propagating edge waves along a buoyancy discontinuity (McIntyre,
2008). However, the difficulty here is that a jump discontinuity in the
buoyancy field results in infinite velocities over constant or increasing
stratification. In addition, the relationship of the along jet waves in the
staircase limit to the non-linear zonons found by Sukoriansky et al. (2008)
remains unclear.

The divergence of the velocity at a buoyancy discontinuities raises a second
question. Is there a limit to how close the staircase limit can be approached?
In barotropic dynamics, the velocity remains finite at a jump continuity in
the vorticity, and, in this case, Scott & Dritschel (2012) report that a
vorticity staircase case can be approached arbitrarily. Whether this result
continues to hold for arbitrarily sharp buoyancy gradients and arbitrarily
large zonal velocities is not clear. Because the rms velocity seems to
converge for arbitrarily sharp staircases, even for the most local inversion
relations we considered, there may not be any energetic reason precluding
arbitrarily sharp buoyancy gradients.

Finally, there remains the question of how relevant these results are for the
upper ocean, which, in addition to surface buoyancy gradients, has interior
potential vorticity gradients as well. In particular, our neglect of the
$\beta$-effect limits the direct relevance of this model to the upper ocean.
Whether surface buoyancy staircases can emerge under more realistic oceanic
conditions requires further investigation.

[Acknowledgements]The author thanks Stephen Griffies for useful discussions
throughout the project as well as helpful comments that improved the clarity
of manuscript.

[Funding]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 author and do not necessarily reflect the views of the National Oceanic
and Atmospheric Administration, or the U.S. Department of Commerce.

[Declaration of interests] The author reports no conflict of interest.

[Data availability statement]The data that support the findings of this study
are available within the article.

[Author ORCID]H. Yassin, https://orcid.org/0000-0003-1742-745X

## References

  * Abernathey et al. (2019) Abernathey et al. 2019 pyqg/pyqg: v0.3.0. https://doi.org/10.5281/zenodo.3551326.
  * Arbic & Flierl (2003) Arbic, Brian K. & Flierl, Glenn R. 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.
  * 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.
  * Chekhlov et al. (1996) Chekhlov, Alexei, Orszag, Steven A., Sukoriansky, Semion, Galperin, Boris & Staroselsky, Ilya 1996 The Effect of Small-Scale Forcing on Large-Scale Structures in Two-Dimensional Flows. Physica D: Nonlinear Phenomena 98, 321–334, https://doi.org/10.1016/0167-2789(96)00102-9.
  * Danilov & Gryanik (2004) Danilov, S. & Gryanik, V. M. 2004 Barotropic Beta-Plane Turbulence in a Regime with Strong Zonal Jets Revisited. J. Atmos. Sci. 61, 2283–2295, https://doi.org/10.1175/1520-0469(2004)061<2283:BBTIAR>2.0.CO;2.
  * Danilov & Gurarie (2002) Danilov, Sergey & Gurarie, David 2002 Rhines scale and spectra of the $\beta$-plane turbulence with bottom drag. Physical Review E 65, 067301, https://doi.org/10.1103/PhysRevE.65.067301.
  * Danilov & Gurarie (2004) Danilov, Sergey & Gurarie, David 2004 Scaling, spectra and zonal jets in beta-plane turbulence. Phys. Fluids 16, 2592–2603.
  * Dritschel & McIntyre (2008) Dritschel, D. G. & McIntyre, M. E. 2008 Multiple Jets as PV Staircases: The Phillips Effect and the Resilience of Eddy-Transport Barriers. J. Atmos. Sci. 65, 855–874, https://doi.org/10.1175/2007JAS2227.1.
  * Dritschel & Scott (2011) Dritschel, D. G. & Scott, R. K. 2011 Jet sharpening by turbulent mixing. Philos. Trans. Roy. Soc. A 369, 754–770, https://doi.org/10.1098/rsta.2010.0306.
  * Dunkerton & Scott (2008) Dunkerton, Timothy J. & Scott, Richard K. 2008 A Barotropic Model of the Angular Momentum-Conserving Potential Vorticity Staircase in Spherical Geometry. J. Atmos. Sci. 65, 1105–1136, https://doi.org/10.1175/2007JAS2223.1.
  * González-Haro & Isern-Fontanet (2014) González-Haro, C. & Isern-Fontanet, J. 2014 Global ocean current reconstruction from altimetric and microwave SST measurements. J. Geophys. Res.: Oceans 119, 3378–3391, https://doi.org/10.1002/2013JC009728.
  * Held et al. (1995) Held, Isaac M., Pierrehumbert, Raymond T., Garner, Stephen T. & Swanson, Kyle L. 1995 Surface quasi-geostrophic dynamics. J. Fluid Mech. 282, 1–20, https://doi.org/10.1017/S0022112095000012.
  * Huang et al. (2000) Huang, Huei-Ping, Galperin, Boris & Sukoriansky, Semion 2000 Anisotropic spectra in two-dimensional turbulence on the surface of a rotating sphere. Phys. Fluids 13, 225–240, https://doi.org/10.1063/1.1327594.
  * Isern‐Fontanet et al. (2006) Isern‐Fontanet, J., Chapron, B., Lapeyre, G. & Klein, P. 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, Jordi, Lapeyre, Guillaume, Klein, Patrice, Chapron, Bertrand & Hecht, Matthew W. 2008 Three-dimensional reconstruction of oceanic mesoscale currents from surface information. J. Geophys. Res.: Oceans 113, https://doi.org/10.1029/2007JC004692.
  * Iwayama & Watanabe (2010) Iwayama, Takahiro & Watanabe, Takeshi 2010 Green’s function for a generalized two-dimensional fluid. Physical Review E 82, https://doi.org/10.1103/PhysRevE.82.036307.
  * LaCasce & Mahadevan (2006) LaCasce, J. H. & Mahadevan, A. 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, Guillaume 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, Guillaume 2017 Surface Quasi-Geostrophy. Fluids 2, 7, https://doi.org/10.3390/fluids2010007.
  * Lapeyre & Klein (2006) Lapeyre, G. & Klein, P. 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.
  * Maltrud & Vallis (1991) Maltrud, M. E. & Vallis, G. K. 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.
  * Maximenko et al. (2005) Maximenko, Nikolai A., Bang, Bohyun & Sasaki, Hideharu 2005 Observational evidence of alternating zonal jets in the world ocean. Geophys. Res. Lett. 32, https://doi.org/10.1029/2005GL022728.
  * McIntyre (2008) McIntyre, M. E. 2008 Potential-vorticity inversion and the wave-turbulence jigsaw: some recent clarifications. In Advances in Geosciences, , vol. 15, pp. 47–56. Copernicus GmbH, https://doi.org/10.5194/adgeo-15-47-2008.
  * Miracca-Lage et al. (2022) Miracca-Lage, Mariana, González-Haro, Cristina, Napolitano, Dante Campagnoli, Isern-Fontanet, Jordi & Polito, Paulo Simionatto 2022 Can the Surface Quasi-Geostrophic (SQG) Theory Explain Upper Ocean Dynamics in the South Atlantic? J. Geophys. Res.: Oceans 127, e2021JC018001, https://doi.org/10.1029/2021JC018001.
  * Pierrehumbert et al. (1994) Pierrehumbert, Raymond T., Held, Isaac M. & Swanson, Kyle L. 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., Zabusky, N. J. & Flierl, G. R. 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, Bo, Chen, Shuiming, Klein, Patrice, Torres, Hector, Wang, Jinbo, Fu, Lee-Lueng & Menemenlis, Dimitris 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, Bo, Chen, Shuiming, Klein, Patrice, Ubelmann, Clement, Fu, Lee-Lueng & Sasaki, Hideharu 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.
  * Rhines (1975) Rhines, Peter B. 1975 Waves and turbulence on a beta-plane. J. Fluid Mech. 69, 417–443, https://doi.org/10.1017/S0022112075001504.
  * Scott et al. (2008) Scott, Robert B., Arbic, Brian K., Holland, Christina L., Sen, Ayon & Qiu, Bo 2008 Zonal versus meridional velocity variance in satellite observations and realistic and idealized ocean circulation models. Ocean Modelling 23, 102–112, https://doi.org/10.1016/j.ocemod.2008.04.009.
  * Scott et al. (2022) Scott, R. K., Burgess, B. H. & Dritschel, D. G. 2022 On the spacing of meandering jets in the strong-stair limit. J. Fluid Mech. 930, https://doi.org/10.1017/jfm.2021.898.
  * Scott & Dritschel (2012) Scott, Richard K. & Dritschel, David G. 2012 The structure of zonal jets in geostrophic turbulence. J. Fluid Mech. 711, 576–598, https://doi.org/10.1017/jfm.2012.410.
  * Scott & Dritschel (2019) Scott, Richard K. & Dritschel, David G. 2019 Zonal Jet Formation by Potential Vorticity Mixing at Large and Small Scales. In Zonal Jets: Phenomenology, Genesis, and Physics, p. 550. Cambridge, U.K.: Cambridge University Press.
  * Smith et al. (2002) Smith, K. S., Boccaletti, G., Henning, C. C., Marinov, I., Tam, C. Y., Held, I. M. & Vallis, G. K. 2002 Turbulent diffusion in the geostrophic inverse cascade. J. Fluid Mech. 469, 13–48, https://doi.org/10.1017/S0022112002001763.
  * Sukhatme & Smith (2009) Sukhatme, Jua & Smith, Leslie M. 2009 Local and nonlocal dispersive turbulence 21, https://doi.org/10.1063/1.3141499.
  * Sukoriansky et al. (2007) Sukoriansky, Semion, Dikovskaya, Nadejda & Galperin, Boris 2007 On the Arrest of Inverse Energy Cascade and the Rhines Scale. J. Atmos. Sci. 64, 3312–3327, https://doi.org/10.1175/JAS4013.1.
  * Sukoriansky et al. (2008) Sukoriansky, Semion, Dikovskaya, Nadejda & Galperin, Boris 2008 Nonlinear Waves in Zonostrophic Turbulence. Physical Review Letters 101 (17), 178501, https://doi.org/10.1103/PhysRevLett.101.178501.
  * Vallis & Maltrud (1993) Vallis, Geoffrey K. & Maltrud, Matthew E. 1993 Generation of Mean Flows and Jets on a Beta Plane and over Topography. J. Phys. Oceanogr. 23, 1346–1362, https://doi.org/10.1175/1520-0485(1993)023<1346:GOMFAJ>2.0.CO;2.
  * Wang et al. (2019) Wang, Shihong, Qiao, Fangli, Dai, Dejun & Zhou, Xiaohui 2019 Anisotropy of the sea surface height wavenumber spectrum from altimeter observations. Scientific Reports 9, https://doi.org/10.1038/s41598-019-52328-w.
  * Yassin & Griffies (2022) Yassin, Houssam & Griffies, Stephen M. 2022 Surface quasigeostrophic turbulence in variable stratification. Journal of Physical Oceanography (in review) https://arxiv.org/abs/2110.04242.