# Evidence for multiple Ferrel-like cells on Jupiter

Keren Duer1∗, Nimrod Gavriel1∗, Eli Galanti1, Yohai Kaspi1,
Leigh N. Fletcher2, Tristan Guillot3, Scott J. Bolton4, Steven M. Levin5,
Sushil K. Atreya6, Davide Grassi7, Andrew P. Ingersoll8, Cheng Li6, Liming
Li9,
Jonathan I. Lunine10, Glenn S. Orton5, Fabiano A. Oyafuso5, J. Hunter Waite,
Jr.4

1Department of Earth and Planetary Sciences, Weizmann Institute of Science,
Rehovot, Israel
2School of Physics and Astronomy, University of Leicester, University Road,
Leicester, LE1 7RH, UK
3Universitie Cote d’Azur, OCA, Lagrange CNRS, 06304 Nice, France
4Southwest Research Institute, San Antonio, Texas, TX, USA
5Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena, CA 91109, USA
6Department of Climate and Space Sciences and Engineering, University of
Michigan, Ann Arbor, MI, USA
7Istituto di Astrofisica e Planetologia Spaziali, INAF, Rome, Italy
8California Institute of Technology, Pasadena, California, USA
9University of Houston, Houston, TX, USA
10Department of Astronomy, Cornell University, Ithaca, New York 14853, USA

* These authors contributed equally to this work

##### Key points

  * •

Measurements from multiple instruments of the Juno mission are interpreted to
reveal the meridional circulation beneath Jupiter’s clouds

  * •

16 Jet-paired deep cells, extending to at least 240 bar, are revealed between
latitudes $60^{\circ}{\rm S}$ and $60^{\circ}{\rm N}$, driven by turbulence
similar to Earth’s Ferrel cells

  * •

The findings are supported by modeling the advection of tracers due to the
cells, showing agreement with ${\rm NH}{}_{3}$ data

### Abstract

Jupiter’s atmosphere is dominated by multiple jet streams which are strongly
tied to its 3D atmospheric circulation. Lacking a rigid bottom boundary,
several models exist for how the meridional circulation extends into the
planetary interior. Here we show, collecting evidence from multiple
instruments of the Juno mission, the existence of mid-latitudinal meridional
circulation cells which are driven by turbulence, similar to the Ferrel cells
on Earth. Different than Earth, which contains only one such cell in each
hemisphere, the larger, faster rotating Jupiter can incorporate multiple
cells. The cells form regions of upwelling and downwelling, which we show are
clearly evident in Juno’s microwave data between latitudes $60^{\circ}{\rm S}$
and $60^{\circ}{\rm N}$. The existence of these cells is confirmed by
reproducing the ammonia observations using a simplistic model. This study
solves a long-standing puzzle regarding the nature of Jupiter’s sub-cloud
dynamics and provides evidence for 8 cells in each Jovian hemisphere.

### Plain Language Summary

The cloud layer of Jupiter is divided into dark and bright bands that are
shaped by strong east-west winds. Such winds in planetary atmospheres are
thought to be tied with a meridional circulation. The Juno mission collected
measurements of Jupiter’s atmosphere at various wavelengths, which penetrate
the cloud cover. Here we provide evidence, using the Juno data, of 8 deep
Jovian circulation cells in each hemisphere encompassing the east-west winds,
gaining energy from atmospheric waves, and extending at least to a depth of
hundreds of kilometers. Different than Earth, which has only 1 analogous cell
in each hemisphere, known as a Ferrel cell, Jupiter can contain more cells due
to its larger size and faster spin. To support the presented evidence, we
modeled how ammonia gas would spread under the influence of such cells and
compared it to the Juno measurements. The presented results shed light on the
unseen flow structure beneath Jupiter’s clouds.

###  1 Introduction

Over the last few decades, spacecraft and ground-based observations have
gathered data about Jupiter’s atmosphere, including measurements of cloud
reflectance (Garcıa-Melendo and Sánchez-Lavega,, 2001), winds (Porco et al.,,
2003; Salyk et al.,, 2006; Tollefson et al.,, 2017), composition (Taylor et
al.,, 2004) and lightning flashes (Little et al.,, 1999). Since 2016, the Juno
spacecraft has provided unprecedented measurements that revealed new
information on the deep dynamics of Jupiter (Bolton et al.,, 2017). Gravity
science enabled an accurate mapping of Jupiter’s gravitational field (Iess et
al.,, 2018), resulting in the inference that the zonal jets penetrate $\sim
3000$ km deep (Kaspi et al.,, 2018; Guillot et al.,, 2018), where they
possibly decay due to magnetic drag (Liu et al.,, 2008; Dietrich and Jones,,
2018; Kaspi et al.,, 2020) and may also require the presence of a stable layer
(Christensen et al.,, 2020). The Jovian Infrared Auroral Mapper (JIRAM)
provided measurements of tropospheric species distribution below the cloud
level (Grassi et al.,, 2020). The Microwave Radiometer (MWR) measurements,
inferred as brightness temperature ($T_{{\rm b}}$), revealed the deep ammonia
abundance (Li et al.,, 2017; Oyafuso et al.,, 2020), as well as lightning at a
frequency of 600 MHz (Brown et al.,, 2018). The combination of these
observations allows the essential nature of Jupiter’s deep overturning
circulation to be revealed, as the flows associated with such circulation are
directly related to cloud formation, temperature variations, lightning
occurrences, tracer distributions and turbulence.

Earth’s atmosphere is commonly referred to as possessing a three-cell
meridional structure in each hemisphere (Vallis,, 2017), which can be
recognized in the zonal-averaged velocities. Circulation cells of such nature
are thought to prevail in the atmospheres of terrestrial planets (Read et
al.,, 2018) and were observed, for example, on Mars (Lewis et al.,, 2007) and
Venus (Limaye,, 2007). On the terrestrial planets, the solid surface drag
plays a part in maintaining the circulation in the cells. However, as the
giant planets hold no such surface, the mere possibility of them possessing
meridional circulation cells remained uncertain. Earth’s midlatitudes are
governed by the Ferrel cells, which are driven by atmospheric turbulence,
creating regions of eddy momentum flux convergence at midlatitudes (Vallis,,
2017). These cells accompany the midlatitude jets and are connected to the
cloud structure in Earth’s atmosphere.

The prominent banded structure at the cloud tops of Jupiter’s atmosphere (Fig.
1a) has been observed for centuries (Vasavada and Showman,, 2005). These
reflectivity contrasts are partially aligned (mainly at low latitudes,
Ingersoll et al., (2000)) with the belts and zones (Fig. 1a), defined by the
sign of the zonal-wind vorticity
($\bar{\zeta}=-\nicefrac{{\partial\bar{u}}}{{\partial y}}$ in Fig. 1c, where
$u$ is zonal velocity, $y$ is in the meridional direction and an over-line
represents a zonal mean). Voyager measurements suggested that the zones are
associated with ascending motion, but this was limited to low latitudes due to
its equatorial trajectory and to altitudes above $0.5$ bar (Gierasch et al.,,
1986). The latitudinal profile of the zonal wind, calculated using cloud-
tracking (Garcıa-Melendo and Sánchez-Lavega,, 2001; Porco et al.,, 2003;
Tollefson et al.,, 2017), reveals that the equatorial region is characterized
by a strong eastward flow, while the midlatitudes exhibit alternating jets,
spaced $2-8^{\circ}$ apart in latitude (Fig. 1b). The midlatitude jets are
correlated with the eddy momentum flux convergence (Salyk et al.,, 2006)
($\nicefrac{{-\partial\left(\overline{u^{\prime}v^{\prime}}\right)}}{{\partial
y}}$ in Fig. 1d, where $v$ is meridional velocity and an apostrophe represents
deviations from the zonal mean, i.e. “eddy” terms), implying that the
midlatitude jets are eddy-driven (Ingersoll et al.,, 2000; Young and Read,,
2017), similar to the jets within Earth’s Ferrel cells (Schneider,, 2006;
Vallis,, 2017). To illustrate the relation between the jets and the eddies,
regions of positive (negative) vorticity gradient,
$\nicefrac{{\partial\zeta}}{{\partial
y}}=-\nicefrac{{\partial^{2}\bar{u}}}{{\partial y^{2}}}$, at midlatitudes, are
marked by light red (blue) bands (Fig. 1), where counter-clockwise (clockwise)
Ferrel-like circulation cells are expected in the northern hemisphere (NH).
Similar circulations, but in opposite directions, apply for the southern
hemisphere (SH). Evidence for vertical motion comes also from observations of
lightning flashes (Little et al.,, 1999; Porco et al.,, 2003; Brown et al.,,
2018), suggesting updrafts in cyclonic belt regions (e.g., Fig. 1e).

Additional information regarding Jupiter’s deep atmosphere can be obtained by
probing Jupiter’s interior at microwave frequencies. Juno’s MWR has 6
microwave channels (Janssen et al.,, 2017), each measuring the atmospheric
$T_{{\rm b}}$ at a different depth (Janssen et al.,, 2017; Bolton et al.,,
2017; Oyafuso et al.,, 2020; Fletcher et al.,, 2021), and collectively
covering the range between $\sim 0.7$ and $\sim 240$ bar (Fig. 1g,h, see also
supporting information - SI). $T_{{\rm b}}$ measurements are affected by both
ammonia abundance and temperature [and water in the case of the longest
wavelengths, Li et al., (2017, 2020); Fletcher et al., (2021)]. If the
latitudinal gradients of $T_{{\rm b}}$ were primarily driven by temperature
changes, then thermal wind balance implies that the midlatitude jets
strengthen from the cloud-deck to about $\sim 8$ bar, and then decay slowly
towards the interior (Fletcher et al.,, 2021). However, interpreting $T_{{\rm
b}}$ as temperature would also imply that the equatorial wind double its
magnitude below the cloud level (Bolton et al.,, 2017), which is inconsistent
with gravity constraints (Duer et al.,, 2020). Thus, the latitudinal variation
of $T_{{\rm b}}$ is probably governed by ammonia opacity, resulting in a map
of ammonia abundance (Li et al.,, 2017), and implying that the zonal winds are
nearly barotropic (Fletcher et al.,, 2021). The overall ammonia structure,
supported also by earlier observations (de Pater et al.,, 2001), reveals
stratification of ammonia with depth, although the mean ammonia profile
changes the sign of its vertical gradient at the $\sim 2-8$ bar region (Giles
et al.,, 2017; Li et al.,, 2017; de Pater et al.,, 2019). The atmospheric
depletion and stratification of ammonia is likely linked to small-scale storm
activity (Guillot et al., 2020a, ; Guillot et al., 2020b, ), where water-
ammonia hail, forming around the 1-bar level, falls below the water-cloud base
and releases ammonia and water at altitudes below 10 bar (Guillot et al.,
2020b, ). Additional measurements of ammonia come from Juno’s JIRAM, which
evaluated the ammonia distribution at a depth of $\sim 5-6$ bar (Grassi et
al.,, 2010, 2020) (Fig. 1f), indicating, as the MWR measurements, that ammonia
varies with latitude. These variations are the key observation for this study,
as ammonia anomalies (deviations from the isobaric mean) can reveal details
about Jupiter’s overturning circulation (Ingersoll et al.,, 2017; Fletcher et
al.,, 2021; Lee and Kaspi,, 2021).

Figure 1:  Observations of Jupiter’s atmosphere. (a) Image of Jupiter’s clouds
(longitudes $69-87^{\circ}$) taken by JunoCam on Dec. 26th 2019 during
perijove 24 (image credit: NASA/JPL/SwRI/MSSS/Gerald Eichstaedt/John Rogers),
with the traditional “dark” belts (“bright” zones) defined as regions of
cyclonic (anticyclonic) vorticity, identified below as ’B’ (’Z’). (b)
Jupiter’s zonally averaged zonal wind [${\pm 15\;{\rm m\,\,s^{-1}}}$] measured
by the Hubble space telescope on December 11, 2016, during Juno’s third
perijove (Tollefson et al.,, 2017). (c) The zonally averaged vorticity $[{\rm
s}^{-1}]$, calculated from the zonal wind profile (panel b). Black dots
represent local extrema in the midlatitudes. (d) Eddy momentum flux
convergence [${\pm 2\times 10^{-6}\;\rm m\,\,s^{-2}}$] calculated from 58
image pairs taken by Cassini during its Jupiter flyby in December, 2000 (Salyk
et al.,, 2006). (e) Lightning detections $[{\rm s}^{-1}]$ by Juno’s MWR during
perijoves 1-8 (yellow, left axis, Brown et al.,, 2018) and number of lightning
storms detected by the Cassini during its flyby (blue, right axis, Porco et
al.,, 2003). (f) Distribution of ammonia [volume mixing ratio] and its mean
(dashed) at a depth of $\sim 6$ bar, measured by Juno’s JIRAM during perijoves
1-15 (Grassi et al.,, 2020). (g) Nadir $T_{{\rm b}}$ [${}^{\circ}{\rm K}$]
(color) interpolated between pressure levels of $0.7$ and $240$ bar (vertical
axis), measured by Juno’s MWR during perijoves 1-12 (Oyafuso et al.,, 2020).
(h) Reconstructed MWR Brightness temperature at midlatitudes. A frequency
filter is applied according to Eq. S7. The standard deviation of each channel
and latitude is available in Fig. S3 and Fig. S4. It can be seen that $T_{{\rm
b}}$ changes its trend at the borders between cells, consistent with the
Ferrel-like cells hypothesis. (b-h) Light red (blue) bands in the midlatitudes
indicate regions of positive (negative) vorticity gradients.

###  2 Ammonia anomalies due to vertical advection

In the presence of a stable vertical ammonia concentration gradient, advection
by the vertical branches of a meridional circulation can affect the
concentration distribution, potentially leading to steady anomalies.
Therefore, the wavy structure of Jupiter’s ammonia distribution (Fig. 1f-h)
can be explained by the presence of meridional circulation cells. On Jupiter,
as condensation of ammonia is expected only at the upper levels of the
atmosphere ($0.5-1$ bar), the ammonia concentration at those levels should be
lower than at depth (Fletcher et al.,, 2020). In addition, precipitation,
small-scale turbulence, thermochemical and chemical reactions, and diffusion
are also expected to determine the vertical ammonia distribution ($M_{{\rm
a}}$) (Guillot et al., 2020b, ). The $M_{{\rm a}}$ profile estimated from the
MWR (Li et al.,, 2017) reveals a local minimum at $\sim 6$ bar (Fig. 2a). This
profile is used in this study as the background state, to explain the ammonia
anomalies.

Here, we focus on two regions with distinctly different deep dynamics: the
equatorial region (planetocentric latitudes $20^{\circ}{\rm S}$ to
$20^{\circ}{\rm N}$), where superrotation is assumed to be fueled by eddy
momentum fluxes perpendicular to the spin axis (e.g., Busse,, 2002), and
midlatitudes ($60^{\circ}{\rm S}$ to $20^{\circ}{\rm S}$ and $20^{\circ}{\rm
N}$ to $60^{\circ}{\rm N}$), where alternating jets are postulated to be
driven by horizontal eddies associated with mass-transporting meridional cells
(e.g., Salyk et al.,, 2006; Schneider and Liu,, 2009; Young et al.,, 2019).

We begin with the midlatitudes, where the meridional cells are mechanically
driven (see below) by turbulence, similar to Earth’s Ferrel cells, which form
as a consequence of atmospheric waves breaking in midlatitudes (Vallis,,
2017). Unlike the largely baroclinic midlatitudes of Earth, which result in
mostly non-mass-transporting Ferrel cells (Juckes,, 2001; Vallis,, 2017), the
predominantly barotropic flows on Jupiter (at the depth range associated with
the MWR measurements) (Kaspi et al.,, 2018, 2020; Galanti and Kaspi,, 2021)
may allow mass-transporting meridional cells (see SI). Consistently, deep
convection models of Jupiter also show barotropic flows (e.g., Busse,, 1976;
Aurnou and Olson,, 2001). The upper branch of Earth’s Ferrel cells consists of
a balance between the Coriolis force and the eddy momentum flux convergence,

$-f\bar{v}=-\frac{\partial\left(\overline{u^{{}^{\prime}}v^{{}^{\prime}}}\right)}{\partial
y},$ (1)

where $f$ is the Coriolis parameter. This upper branch balance, which is the
leading order balance of the steady state zonal mean zonal momentum equation,
is expected to hold within the equivalent cells on Jupiter (see SI). This
balance can hold down to a depth of only a few bars, as inferred from energy
considerations (Liu and Schneider,, 2010), implying flows from belts to zones
within the cloud layer of the Jovian atmosphere. In the lower branch of the
terrestrial Ferrel cells, the balance is between the Coriolis force and a
surface drag (Vallis,, 2017). Since the Jovian atmosphere lacks a bottom solid
boundary, surface drag cannot act to oppose the Coriolis force, although it
has been suggested that if the cells extend as deep as the jets (Kaspi et
al.,, 2020), the Lorentz force can act as a magnetic drag (Liu et al.,, 2008;
Liu and Schneider,, 2010; Wicht et al.,, 2019). Another possible mechanism
that allows the jets to be barotropic in the upper atmosphere and decay in the
interior is the presence of a stable layer, as was shown lately in several
studies (Debras and Chabrier,, 2019; Christensen et al.,, 2020; Wicht and
Gastine,, 2020).

Figure 2:  Schematics of Jupiter’s meridional circulation as inferred from the
ammonia distribution. (a) The vertical structure of the meridionally averaged
ammonia concentration ($M_{{\rm a}}$) [ppm], as interpreted from the $T_{{\rm
b}}$ data (Li et al.,, 2017). (b) Illustration of a midlatitude Ferrel-like
circulation cell (blue arrows) in the NH, looking from east towards the west.
The cells are accompanied by an eddy-driven barotropic jet ($\bar{u}$), which
peaks at the center of the cell (beige contours). Ammonia constant-
concentration lines are illustrated with orange shades (according to the
ammonia vertical profile from panel a). Dashed orange lines are deviations
from $M_{{\rm a}}$, driven by vertical advection. The return flow of the cell,
illustrated by a dashed blue arrow, lies at an unknown depth. An oppositely
directed upper cell, as suggested by pre-Juno measurements (Ingersoll et al.,,
2000; Showman and de Pater,, 2005; Fletcher et al.,, 2020), is demonstrated by
dashed transparent arrows. $p$ is pressure, taken as a vertical coordinate.
(c) A closer look at the region where the rising air advects ammonia-poor
fluid to an ammonia-rich layer, associated with pressure levels between $1.5$
and $6$ bar. (d) Here, rising gas drags higher ammonia concentration to a
lower ammonia concentration region, associated with pressure levels deeper
than $6$ bar. (e) A cross section of Jupiter’s equatorial upwelling
($\overline{w}$), associated with a superrotating jet ($\bar{u}$, green
contours), leading to ammonia concentration maximum. The equatorial $M_{{\rm
a}}$ (orange contours) is assumed to decrease with radius (Fig. S8). $R_{{\rm
J}}$ is Jupiter’s radius and $\vartheta$ and $r$ are the latitudinal and
radial directions, respectively.

The direction of each Ferrel-like cell corresponds to the direction of the
respective midlatitude jet. Eastward (westward) jets are located in cells of
eastward momentum flux convergence (divergence) implying (Eq. 1) a
counterclockwise (clockwise) circulation in the NH, and a clockwise
(counterclockwise) circulation in the SH (Fig. 1b,d). The upper branch of the
Ferrel-like cells may coincide with the lower branch of stacked upper cells
with an opposite circulation (Ingersoll et al.,, 2000; Showman and de Pater,,
2005) (dashed transparent lines in Fig. 2b), and therefore may share the same
balance (Eq. 1). Indications for the upper cells come from temperature and
shallow tracer distributions (Gierasch et al.,, 1986; Fletcher et al.,, 2016;
de Pater et al.,, 2019). Similar to the balance describing the deeper branch
of the lower cells, the upper branch of the upper cells requires a drag force,
which may result from breaking of atmospheric waves (Gierasch et al.,, 1986;
Ingersoll et al.,, 2021).

Figure 3: Pearson correlation coefficients as a function of latitude. The
correlations exemplify the relations presented in Fig. 2. (a) and (c),
correlations calculated between the zonal jets $\left(\bar{u}\right)$ and the
$T_{{\rm b}}$ meridional gradients $\left(\partial_{y}T_{{\rm b}}\right)$,
adjusted by the sign of the vertical gradient of $M_{{\rm a}}$, in the six MWR
channels, for the SH and NH, respectively. (b) Correlations computed between
the zonal jet velocity and $T_{{\rm b}}$ ($\bar{u}\propto-T_{{\rm b}}$). (d)
and (f), Correlations between the zonal velocity and the ammonia abundance
gradients ($\partial_{y}{\rm NH}_{3}$), measured by JIRAM, in the SH and NH,
respectively. (e) Correlations between the zonal velocity and the ammonia
abundance from JIRAM ($\bar{u}\propto{\rm NH}_{3}$). (g) and (h) Correlations
between the zonal velocity and the eddy momentum flux convergence
($\bar{u}\propto-\partial_{y}\left(\overline{u^{\prime}v^{\prime}}\right)$),
in the SH and NH, respectively. Gray dots represent correlations that are not
statistically significant (confidence level $95\%$) and gray regions show
where measurements were not available. No data is available for the eddy
fluxes in the equatorial region, as evaluating them requires measurements of
the vertical winds (see SI), which are yet to be achieved.

The background ammonia profile is skewed by the vertical branches of the cells
(dashed orange lines in Fig. 2b), maximizing the ammonia meridional gradient
where the jet velocity peaks (i.e., in the middle of the cell). This means
that a correlation (along isobars) is expected between the zonal jets and the
meridional gradient of the ammonia concentration at midlatitudes (Duer et
al.,, 2020; Fletcher et al.,, 2021). However, since the vertical gradient of
$M_{{\rm a}}$ changes with depth (Li et al.,, 2017) (Fig. 2a), the nature of
the correlation should change as well, as illustrated in Fig. 2c,d. These
simple considerations motivate the examination of the correlation between
$\bar{u}$ and $\partial_{y}m_{a}$ $\left(-\partial_{y}T_{{\rm b}}\right)$ in
midlatitudes (Fig. 3a,c, also see SI). Note that $T_{{\rm b}}$ corresponds
inversely to ammonia abundance at a certain pressure level (Li et al.,, 2017).
For a deep-wind estimate, we use the measured cloud-level winds (Tollefson et
al.,, 2017) projected inward in a direction parallel to the axis of rotation,
without any change in magnitude (in the upper $240$ bar), as implied by
gravity measurement constraints (Galanti and Kaspi,, 2021; Galanti et al.,,
2021). The correlations are performed using a $4^{\circ}$ latitudinal bin (see
SI). By this, the suggested correlations in Fig. 2 can be tested locally,
rather than over an entire hemisphere (Duer et al.,, 2020; Fletcher et al.,,
2021).

At midlatitudes the overall positive correlations for MWR channels 1-5
indicate the existence of Ferrel-like cells at depths between $1.5$ and $240$
bar (Fig. 3a,c). The positive correlation with ammonia estimates by JIRAM
(Fig. 3d,f) further strengthens the prominence of the proposed cells. Channel
6 ($\sim 0.7$ bar) exhibits negative correlations in midlatitudes (Fig. 3a,c),
implying that the deep cells do not extend higher than $\sim 1$ bar, and
support the existence of counter-rotating cells above that level (Ingersoll et
al.,, 2000; Showman and de Pater,, 2005; Fletcher et al.,, 2020). To verify
that the relation shown in Eq. 1 holds in the cells as illustrated in Fig. 2b,
regional correlations between $\bar{u}$ and
$-\partial_{y}\left(\overline{u^{{}^{\prime}}v^{{}^{\prime}}}\right)$ are
shown for midlatitudes (Fig. 3g,h). This positive correlation further
strengthens the existence of the Ferrel-like cells, where converging eddy
momentum fluxes are the source of momentum. Overall, the correlation analysis
reveals multiple deep Ferrel-like cells, extending from $\sim 1$ bar to at
least $240$ bar.

The equatorial region of Jupiter, characterized by a wide eastward jet, needs
to be treated differently. Gravity analysis reveals that Jupiter’s interior
(deeper than $\sim 3000\,\,{\rm km}$) is rotating as a rigid body (Guillot et
al.,, 2018). Extending the zonal wind along the direction of the spin axis
thus separates the equatorial region ($17^{\circ}\,{\rm S}$ to
$17^{\circ}\,{\rm N}$) from the truncated cells at midlatitudes (see SI). The
superrotating wind at low latitudes requires a source of momentum (Imamura et
al.,, 2020). Theories for such sources include meridional (Potter et al.,,
2014; Laraia and Schneider,, 2015) and vertical (Aurnou and Olson,, 2001;
Busse,, 2002; Christensen,, 2002; Heimpel et al.,, 2005; Kaspi et al.,, 2009;
Dietrich and Jones,, 2018), propagation of waves. For the vertical case,
several studies have shown that an equatorial superrotation in giant planets
can be driven by eddy momentum fluxes perpendicular to the axis of rotation
(Heimpel et al.,, 2005; Kaspi et al.,, 2009; Gastine et al.,, 2014). These
fluxes transfer momentum outwards and lead to a mean upwelling at the
equatorial region (see SI). Such an upwelling should lead to a maximum
concentration anomaly of any stably stratified matter (Fig. 2e). In the
equatorial region, the minimum in $M_{{\rm a}}$ around $\sim 6$ bar nearly
vanishes (Fig. S8), suggesting that the positive ammonia anomalies at the
equatorial region (Fig. S2a) are due to the upwelling from deep. To examine
this, at the equator, the correlation is calculated between the zonal velocity
($\bar{u}$) and the ammonia concentration itself ($-T_{{\rm b}}$ for the MWR
or ${\rm NH_{3}}$ for the JIRAM measurements). Using a regional correlation
analysis (see SI), it is apparent that the correlations are largely positive
at all depths (Fig. 3b,e), implying that an equatorial upwelling is dominant
from the cloud deck and down to at least $240$ bar. Very close to the equator,
the correlation is negative due to the local minimum in the zonal velocity
(Fig. 1b).

###  3 Model reconstruction of Jupiter’s ammonia distribution

To further validate that the positive correlations shown in Fig. 3a,c are
indeed due to the existence of meridional circulation cells, we reconstruct
the measured variations using a simplified advection-relaxation model.
Beginning with a steady-state zonal-mean conservation of species equation for
ammonia, assuming that diffusion terms are small, the leading-order balance is

$\bar{w}\left(\vartheta,r\right)\frac{\partial
m_{a}\left(\vartheta,r\right)}{\partial
r}+\bar{v}\left(\vartheta,r\right)\frac{\partial
m_{{a}}\left(\vartheta,r\right)}{r\partial\vartheta}=-G\left(r\right)\left(m_{a}\left(\vartheta,r\right)-M_{{\rm
a}}\left(r\right)\right),$ (2)

where $\bar{w}$ is the zonally averaged radial velocity, and $m_{a}$, the
variable solved for by the model, is the (zonal-mean) molar fraction of
ammonia. $M_{{\rm a}}$ is the ammonia concentration averaged over isobaric
surfaces (Fig. 2a), and $G$ is the inverse of a Newtonian relaxation
timescale. The two terms on the left-hand side represent advection by the mean
circulation, and the right-hand side term is a source term parameterized as a
simple Newtonian relaxation of ammonia. This relaxation term is assumed to
include all the processes resulting in the observed $M_{{\rm a}}$ as it acts
against local anomalies toward this mean vertical structure. To qualitatively
illustrate how the Ferrel-like cells’ footprint might appear in the ammonia
distribution map ($m_{a}$), we solve the advection-relaxation balance shown in
Eq. 2, for the midlatitudes between $1.5$ and $240$ bar (see SI). As the
balance in Eq. 2 indicates, it is assumed that the relaxation time scale
$\left(G^{-1}\right)$ is such that the advection and relaxation terms balance
each other.

Figure 4: Jupiter’s Ammonia distribution driven by an array of circulation
cells. (a) The normalized vertical zonal-mean wind ($\bar{w}$), as a function
of latitude and pressure, used in the model. Red and green contours are upward
and downward winds, respectively. (b) The ammonia anomalies reconstructed from
the data [ppm]. Here the vertical mean profile $M_{{\rm a}}$ is removed from
the ammonia map $m_{a}$. (c) The ammonia anomalies [ppm] produced by the
advection-relaxation model. (d) The full reconstructed ammonia map
$m_{a}^{{\rm(data)}}$ [ppm]. Arrows represent the direction of the cells’
vertical and meridional winds. (a-d) Red and green vertical lines are the
locations of the upward and downward branches of the cells, respectively. For
reference, light red (blue) bands indicate regions of positive (negative)
vorticity gradient as in Fig. 1. The vertical axis is truncated at $100$ bar
as $M_{{\rm a}}$ becomes largely uniform beyond this depth, thereby
suppressing footprints of advection.

The zonally averaged velocity components $\left(\bar{v},\bar{w}\right)$ of the
circulation cells, necessary for setting the advection terms of Eq. 2, can be
projected from the available wind data according to the outline illustrated in
Fig. 2b. Specifically, we relate between the circulation cells and the wind
data corresponding to the following assumptions (see also SI). The borders
between the cells are set at local extrema of the observed cloud-level
vorticity, the directions of the circulation cells are set according to the
directions of the jets in the middle of each cell, and the strength of the
circulation in each cell is set by the measured eddy momentum flux convergence
along the cell (Fig. S5). As the three terms in Eq. 2 should be proportional,
but cannot be uniquely determined, the values of $\bar{w}$, $\bar{v}$ and $G$
are normalized (Fig. 4a). This normalization means that while the model cannot
produce absolute values of winds due to unmeasured quantities, it can predict
qualitatively how these velocities would be structured spatially and what
should be their relative magnitudes, which are sufficient for assessing the
existence of the cells. Using scaling arguments, the value of Jupiter’s static
stability has recently been estimated to be in the order of
$10^{-2}\;\rm{s^{-1}}$ (Lee and Kaspi,, 2021), which can provide a further
step towards estimating the magnitude of the velocities in the cells.

The described wind scheme results in upwellings (downwellings) on the poleward
(equatorward) sides of eastward jets (Fig. 4a). The cells are reversed for
westward jets. Finally, as a benchmark for the model results, derivation of
the ammonia abundance ($m_{a}^{({\rm data})}$) from the measured $T_{{\rm b}}$
between the latitudes $60^{\circ}$S and $60^{\circ}$N is implemented (see SI,
Fig. S2). As the depth of the cells, the width of their branches and the
parameter $G$ are unknown, an optimization procedure is performed for
determining these parameters to best match the data (see SI). To ensure that
this procedure does not influence the qualitative nature of the results, Eq. 2
is also solved with a predefined physically-oriented set of parameters (Figs.
S6 and S7).

Using the above assumptions, we solve Eq. 2 to predict the ammonia map
($m_{a}^{({\rm model})}$), and compare it to $m_{a}^{({\rm data})}$ (Fig. 4).
We stress that the latitudinal variations appearing in the results (Fig. 4c),
stem only from the cloud-level wind observations without any assumption on the
meridional ammonia variation. For a clear comparison between the $m_{a}^{({\rm
data})}$ and $m_{a}^{({\rm model})}$, $M_{{\rm a}}$ is subtracted from both,
such that only anomalies are visible (Fig. 4b,c). Around 10 bar (Fig. 4b),
where $M_{{\rm a}}$ greatly increases with depth, enriched (depleted) ammonia
anomalies appear where upwellings (downwellings) are expected (Fig. 4a). These
features flip sign around the $6$-bar level, where $M_{{\rm a}}$ decreases
with depth. These elements are captured well by the advection-relaxation model
(Fig. 4c). In the SH, all $18$ anomalies apparent in the observations have a
counterpart of similar sign, shape and position in the model results,
suggesting the existence of 8 meridional circulation cells. This agreement
validates that advection by the vertical branches of the cells is the main
contributor in the creation of the observed ammonia anomalies. In the NH,
similar results are achieved, although the cells are slightly less coherent,
perhaps due to unexplained differences between the perijoves in the NH
midlatitudes [Fig. S3, Oyafuso et al., (2020); Fletcher et al., (2021)], which
might mask the cells’ footprints in the MWR data. Nevertheless, the lightning
data reinforces the existence of the NH cells, as lightning peaks are aligned
with the rising branch of the cells at the poleward side of the eastward jets
(Fig. 1e, Fig. S1), which combined with the MWR data (Fig. 4b-d) provide
indication for 8 northern cells. Additional NH centered perijoves during the
Juno extended mission may provide data to better constrain the NH cells. For
more intuition, one can look at the full ammonia map (Fig. 4d), where iso-
concentration lines are pulled up and down by the vertical winds (as
schematically illustrated in Fig. 2b), emphasizing the locations of the 16
eddy-driven cells evident in the MWR data.

###  4 Discussion

The identified array of alternating cells in midlatitudes, along with the
equatorial upwelling, are key features in the meridional overturning
circulation of the Jovian atmosphere (Fig. 5). The cell’s depth that can be
inferred from the MWR measurements is limited to the sensing range ($\sim 240$
bar), and while the midlatitudinal cells are mechanically driven, as the
Ferrel cells in Earth’s troposphere, they are likely to extend deeper into the
planet, as suggested by multiple theoretical studies (e.g., Liu and
Schneider,, 2010; Christensen et al.,, 2020). Similarly, deep meridional
cells, which are mechanically driven, have been suggested to exist on the Sun
(Miesch and Hindman,, 2011).

This study provides an explanation for the observed meridional ammonia
anomalies, given the meridionally averaged vertical ammonia profile. The
consistency of these results suggest that the $T_{{\rm b}}$ latitudinal
variations are dominated by the opacity of a passive tracer, rather than the
kinetic temperature. Note that evidence for the part of the deep cells
extending from $1.5$ to $6$ bar depends on the flip of the background ammonia
gradient (Fig. 2a), and without it these depths might be part of upper inverse
cells (Fletcher et al.,, 2021). The shape of this vertical profile might be
set by precipitation, diffusion, and small-scale mixing, all of which might
change with latitude and depth (Guillot et al., 2020a, ). Nonetheless, the
remarkable agreement between the model and the data, together with the robust
correlation analysis, provide strong evidence that the observed distribution
of ammonia is governed by the existence, number, position and relative
strength of the Ferrel-like circulation cells in Jupiter.

Figure 5:  (b) A figurative cross section of Jupiter’s meridional circulation
and (a) a magnification of the midlatitude circulation cells. The circulation
cells (blue) are axisymmetric in the zonal direction. The pink shell
represents a deep layer characteristic for all depths within the circulation
cells. The white arrows represent alternating jet streams and are symmetric
around the equator for the purpose of clarity. Each jet between latitudes
$20^{\circ}-60^{\circ}\,{\rm S/N}$ is accompanied by a turbulence-driven
circulation cell (blue arrows) in the meridional plane as illustrated in panel
b. The equatorial upwelling associated with the superrotating jet is drawn at
the equator, as illustrated in panel e, as part of a larger possible
equatorial cell (dark blue).

#### acknowledgments

All the data used in this study is publicly available, see Tollefson et al.,
(2017) for the winds data, Salyk et al., (2006) for the eddies data, Li et
al., (2017) for the ammonia data, Brown et al., (2018) for the lightning data
and Oyafuso et al., (2020) for the brightness temperature data.

### References

  * Aurnou and Olson, (2001) Aurnou, J. M. and Olson, P. L. (2001).  Strong zonal winds from thermal convectionin a rotating spherical shell.  Geophys. Res. Lett., 28(13):2557–2559.
  * Bolton et al., (2017) Bolton, S. J., Adriani, A., Adumitroaie, V., Allison, M., Anderson, J., Atreya, S., Bloxham, J., Brown, S., Connerney, J. E. P., DeJong, E., Folkner, W., Gautier, D., Grassi, D., Gulkis, S., Guillot, T., Hansen, C., Hubbard, W. B., Iess, L., Ingersoll, A., Janssen, M., Jorgensen, J., Kaspi, Y., Levin, S. M., Li, C., Lunine, J., Miguel, Y., Mura, A., Orton, G., Owen, T., Ravine, M., Smith, E., Steffes, P., Stone, E., Stevenson, D., Thorne, R., Waite, J., Durante, D., Ebert, R. W., Greathouse, T. K., Hue, V., Parisi, M., Szalay, J. R., and Wilson, R. (2017).  Jupiter’s interior and deep atmosphere: The initial pole-to-pole passes with the Juno spacecraft.  Science, 356:821–825.
  * Brown et al., (2018) Brown, S., Janssen, M., Adumitroaie, V., Atreya, S., Bolton, S., Gulkis, S., Ingersoll, A., Levin, S., Li, C., Li, L., Lunine, J., Misra, S., Orton, G., Steffes, P., Tabataba-Vakili, F., Kolmasova, I., Imai, M., Santolik, O., Kurth, W., Hospodarsky, G., Gurnett, D., and Connerney, J. (2018).  Prevalent lightning sferics at 600 megahertz near Jupiter’s poles.  Nature, 558(7708):87–90.
  * Busse, (1976) Busse, F. H. (1976).  A simple model of convection in the Jovian atmosphere.  Icarus, 29:255–260.
  * Busse, (2002) Busse, F. H. (2002).  Convective flows in rapidly rotating spheres and their dynamo action.  Phys. of Fluids., 14:1301–1314.
  * Christensen, (2002) Christensen, U. R. (2002).  Zonal flow driven by strongly supercritical convection in rotating spherical shells.  J. Comp. Phys., 470:115–133.
  * Christensen et al., (2020) Christensen, U. R., Wicht, J., and Dietrich, W. (2020).  Mechanisms for limiting the depth of zonal winds in the gas giant planets.  Astrophys. J., 890(1):61.
  * de Pater et al., (2001) de Pater, I., Dunn, D., Romani, P., and Zahnle, K. (2001).  Reconciling Galileo Probe Data and Ground-Based Radio Observations of Ammonia on Jupiter.  Icarus, 149(1):66–78.
  * de Pater et al., (2019) de Pater, I., Sault, R. J., Wong, M. H., Fletcher, L. N., DeBoer, D., and Butler, B. (2019).  Jupiter’s ammonia distribution derived from VLA maps at 3-37 GHz.  Icarus, 322:168–191.
  * Debras and Chabrier, (2019) Debras, F. and Chabrier, G. (2019).  New models of Jupiter in the context of Juno and Galileo.  Astrophys. J., 872(1):100.
  * Dietrich and Jones, (2018) Dietrich, W. and Jones, C. A. (2018).  Anelastic spherical dynamos with radially variable electrical conductivity.  Icarus, 305:15–32.
  * Duer et al., (2020) Duer, K., Galanti, E., and Kaspi, Y. (2020).  The range of Jupiter’s flow structures that fit the Juno asymmetric gravity measurements.  J. Geophys. Res. (Planets), 125(8).
  * Fletcher et al., (2016) Fletcher, L. N., Greathouse, T. K., Orton, G. S., Sinclair, J. A., Giles, R. S., Irwin, P. G., and Encrenaz, T. (2016).  Mid-infrared mapping of Jupiter’s temperatures, aerosol opacity and chemical distributions with IRTF/TEXES.  Icarus, 278:128–161.
  * Fletcher et al., (2020) Fletcher, L. N., Kaspi, Y., Guillot, T., and Showman, A. P. (2020).  How well do we understand the belt/zone circulation of Giant Planet atmospheres?  Space Sci. Rev., 216(2):1–33.
  * Fletcher et al., (2021) Fletcher, L. N., Oyafuso, F. A., Allison, M., Ingersoll, A., Li, L., Kaspi, Y., Galanti, E., Wong, M. H., Orton, G. S., Duer, K., Zhang, Z., Li, C., Guillot, T., Levin, S. M., and Bolton, S. (2021).  Jupiter’s temperate belt/zone contrasts revealed at depth by Juno microwave observations.  Earth and Space Science Open Archive, page 35.
  * Galanti and Kaspi, (2021) Galanti, E. and Kaspi, Y. (2021).  Combined magnetic and gravity measurements probe the deep zonal flows of the gas giants.  Mon. Not. Roy. Astro. Soc., 501(2):2352–2362.
  * Galanti et al., (2021) Galanti, E., Kaspi, Y., Duer, K., Fletcher, L. N., Ingersoll, A., Cheng, L., Orton, G. S., Guillot, T., M., L. S., and J., B. S. (2021).  Constraints on the latitudinal profile of Jupiter’s deep jets.  Geophys. Res. Lett., 48(9):e2021GL092912.
  * Garcıa-Melendo and Sánchez-Lavega, (2001) Garcıa-Melendo, E. and Sánchez-Lavega, A. (2001).  A study of the stability of jovian zonal winds from HST images: 1995–2000.  Icarus, 152(2):316–330.
  * Gastine et al., (2014) Gastine, T., Wicht, J., Duarte, L. D. V., Heimpel, M., and Becker, A. (2014).  Explaining Jupiter’s magnetic field and equatorial jet dynamics.  Geophys. Res. Lett., 41:5410–5419.
  * Gierasch et al., (1986) Gierasch, P. J., Magalhaes, J. A., and Conrath, B. J. (1986).  Zonal mean properties of Jupiter’s upper troposphere from Voyager infrared observations.  Icarus, 67:456–483.
  * Giles et al., (2017) Giles, R. S., Fletcher, L. N., Irwin, P. G., Orton, G. S., and Sinclair, J. A. (2017).  Ammonia in Jupiter’s troposphere from high-resolution 5 $\mu$m spectroscopy.  Geophys. Res. Lett., 44(21):10–838.
  * Grassi et al., (2010) Grassi, D., Adriani, A., Moriconi, M. L., Ignatiev, N. I., D’Aversa, E., Colosimo, F., Negrão, A., Brower, L., Dinelli, B. M., Coradini, A., and Iccioni, G. (2010).  Jupiter’s hot spots: Quantitative assessment of the retrieval capabilities of future IR spectro-imagers.  Planatary and Space Science, 58(10):1265–1278.
  * Grassi et al., (2020) Grassi, D., Adriani, A., Mura, A., Atreya, S. K., Fletcher, L. N., Lunine, J. I., Orton, G. S., Bolton, S., Plainaki, C., Sindoni, G., Altieri, F., Cicchetti, A., Dinelli, B. M., Filacchione, G., Migliorini, A., Moriconi, M. L., Noschese, R., Olivieri, A., Piccioni, G., Sordini, R., Stefani, S., Tosi, F., and Turrini, D. (2020).  On the spatial distribution of minor species in Jupiter’s troposphere as inferred from Juno JIRAM data.  J. Geophys. Res. (Planets), 125(4).
  * (24) Guillot, T., Li, C., Bolton, S. J., Brown, S. T., Ingersoll, A. P., Janssen, M. A., Levin, S. M., Lunine, J. I., Orton, G. S., Steffes, P. G., and Stevenson, D. J. (2020a).  Storms and the depletion of ammonia in Jupiter: Ii. explaining the Juno observations.  J. Geophys. Res. (Planets), 125(8):e2020JE006404.
  * Guillot et al., (2018) Guillot, T., Miguel, Y., Militzer, B., Hubbard, W. B., Kaspi, Y., Galanti, E., Cao, H., Helled, R., Wahl, S. M., Iess, L., Folkner, W. M., Stevenson, D. J., Lunine, J. I., Reese, D. R., Biekman, A., Parisi, M., Durante, D., Connerney, J. E. P., Levin, S. M., and Bolton, S. J. (2018).  A suppression of differential rotation in Jupiter’s deep interior.  Nature, 555:227–230.
  * (26) Guillot, T., Stevenson, D. J., Atreya, S. K., Bolton, S. J., and Becker, H. N. (2020b).  Storms and the depletion of ammonia in Jupiter: I. microphysics of mushballs.  J. Geophys. Res. (Planets), 125(8):e2020JE006403.
  * Heimpel et al., (2005) Heimpel, M., Aurnou, J., and Wicht, J. (2005).  Simulation of equatorial and high-latitude jets on Jupiter in a deep convection model.  Nature, 438:193–196.
  * Iess et al., (2018) Iess, L., Folkner, W. M., Durante, D., Parisi, M., Kaspi, Y., Galanti, E., Guillot, T., Hubbard, W. B., Stevenson, D. J., Anderson, J. D., Buccino, D. R., Casajus, L. G., Milani, A., Park, R., Racioppa, P., Serra, D., Tortora, P., Zannoni, M., Cao, H., Helled, R., Lunine, J. I., Miguel, Y., Militzer, B., Wahl, S., Connerney, J. E. P., Levin, S. M., and Bolton, S. J. (2018).  Measurement of Jupiter’s asymmetric gravity field.  Nature, 555(7695):220–222.
  * Imamura et al., (2020) Imamura, T., Mitchell, J., Lebonnois, S., Kaspi, Y., Showman, A. P., and Korablev, O. (2020).  Superrotation in planetary atmospheres.  Space Sci. Rev., 216(5):1–41.
  * Ingersoll et al., (2017) Ingersoll, A. P., Adumitroaie, V., Allison, M. D., Atreya, S., Bellotti, A. A., Bolton, S. J., Brown, S. T., Gulkis, S., Janssen, M. A., Levin, S. M., Cheng, L., Liming, L., Lunine, J. I., Orton, G. S., Oyafuso, F. A., and Steffes, P. G. (2017).  Implications of the ammonia distribution on Jupiter from 1 to 100 bars as measured by the Juno microwave radiometer.  Geophys. Res. Lett., 44(15):7676–7685.
  * Ingersoll et al., (2021) Ingersoll, A. P., Atreya, S., Bolton, S. J., Brueschaber, S., Fletcher, L. N., Galanti, E., Kaspi, Y., Levin, S. M., Li, C., Li, L., Lunine, J. I., Orton, G. S., and Waite, H. (2021).  Jupiter’s overturning circulation: Breaking waves take the place of solid boundaries.  Geophys. Res. Lett. in review.
  * Ingersoll et al., (2000) Ingersoll, A. P., Gierasch, P. J., Banfield, D., Vasavada, A. R., and Galileo Imaging Team (2000).  Moist convection as an energy source for the large-scale motions in Jupiter’s atmosphere.  Nature, 403:630–632.
  * Janssen et al., (2017) Janssen, M. A., Oswald, J. E., Brown, S. T., Gulkis, S., Levin, S. M., Bolton, S. J., Allison, M. D., Atreya, S. K., Gautier, D., Ingersoll, A. P., Lunine, J. I., Orton, G. S., Owen, T. C., Steffes, P. G., Adumitroaie, V., Bellotti, A., Jewell, L. A., Li, C., Li, L., Misra, S., Oyafuso, F. A., Santos-Costa, D., Sarkissian, E., Williamson, R., Arballo, J. K., Kitiyakara, A., Ulloa-Severino, A., Chen, J. C., Maiwald, F. W., Sahakian, A. S., Pingree, P. J., Lee, K. A., Mazer, A. S., Redick, R., Hodges, R. E., Hughes, R. C., Bedrosian, G., Dawson, D. E., Hatch, W. A., Russell, D. S., Chamberlain, N. F., Zawadski, M. S., Khayatian, B., Franklin, B. R., Conley, H. A., Kempenaar, J. G., Loo, M. S., Sunada, E. T., Vorperion, V., and Wang, C. C. (2017).  MWR: Microwave Radiometer for the Juno Mission to Jupiter.  Space Sci. Rev., 213(1-4):139–185.
  * Juckes, (2001) Juckes, M. (2001).  A generalization of the transformed Eulerian-mean meridional circulation.  Q. J. R. Meteorol. Soc., 127(571):147–160.
  * Kaspi et al., (2009) Kaspi, Y., Flierl, G. R., and Showman, A. P. (2009).  The deep wind structure of the giant planets: Results from an anelastic general circulation model.  Icarus, 202:525–542.
  * Kaspi et al., (2018) Kaspi, Y., Galanti, E., Hubbard, W. B., Stevenson, D. J., Bolton, S. J., Iess, L., Guillot, T., Bloxham, J., Connerney, J. E. P., Cao, H., Durante, D., Folkner, W. M., Helled, R., Ingersoll, A. P., Levin, S. M., Lunine, J. I., Miguel, Y., Militzer, B., Parisi, M., and Wahl, S. M. (2018).  Jupiter’s atmospheric jet-streams extend thousands of kilometres deep.  Nature, 555:223–226.
  * Kaspi et al., (2020) Kaspi, Y., Galanti, E., Showman, A. P., Stevenson, D. J., Guillot, T., Iess, L., and Bolton, S. J. (2020).  Comparison of the deep atmospheric dynamics of Jupiter and Saturn in light of the Juno and Cassini gravity measurements.  Space Sci. Rev., 216(5):1–27.
  * Laraia and Schneider, (2015) Laraia, A. L. and Schneider, T. (2015).  Superrotation in terrestrial atmospheres.  J. Atmos. Sci., 72(11):4281–4296.
  * Lee and Kaspi, (2021) Lee, S. and Kaspi, Y. (2021).  Towards an understanding of the structure of Jupiter’s atmosphere using the ammonia distribution and the Transformed Eulerian Mean theory.  J. Atmos. Sci., 78(7):2047–2056.
  * Lewis et al., (2007) Lewis, S. R., Read, P. L., Conrath, B. J., Pearl, J. C., and Smith, M. D. (2007).  Assimilation of thermal emission spectrometer atmospheric data during the Mars Global Surveyor aerobraking period.  Icarus, 192(2):327–347.
  * Li et al., (2020) Li, C., Ingersoll, A., Bolton, S., Levin, S., Janssen, M., Atreya, S., Lunine, J., Steffes, P., Brown, S., Guillot, T., Allison, M., Arballo, J., Bellotti, A., Adumitroaie, V., Gulkis, S., Hodges, A., Li, L., Misra, S., Orton, G., Oyafuso, F., Santos-Costa, D., Waite, H., and Zhang, Z. (2020).  The water abundance in Jupiter’s equatorial zone.  Nature Astronomy, 4(6):609–616.
  * Li et al., (2017) Li, C., Ingersoll, A., Janssen, M., Levin, S., Bolton, S., Adumitroaie, V., Allison, M., Arballo, J., Bellotti, A., Brown, S., Ewald, S., Jewell, L., Misra, S., Orton, G., Oyafuso, F., Steffes, P., and Williamson, R. (2017).  The distribution of ammonia on Jupiter from a preliminary inversion of Juno microwave radiometer data.  Geophys. Res. Lett., 44(11):5317–5325.
  * Limaye, (2007) Limaye, S. S. (2007).  Venus atmospheric circulation: Known and unknown.  J. Geophys. Res. (Planets), 112(E4).
  * Little et al., (1999) Little, B., Anger, C. D., Ingersoll, A. P., Vasavada, A. R., Senske, D. A., Breneman, H. H., Borucki, W. J., and Team, T. G. S. (1999).  Galileo images of lightning on Jupiter.  Icarus, 142(2):306–323.
  * Liu et al., (2008) Liu, J., Goldreich, P. M., and Stevenson, D. J. (2008).  Constraints on deep-seated zonal winds inside Jupiter and Saturn.  Icarus, 196:653–664.
  * Liu and Schneider, (2010) Liu, J. and Schneider, T. (2010).  Mechanisms of jet formation on the giant planets.  J. Atmos. Sci., 67:3652–3672.
  * Miesch and Hindman, (2011) Miesch, M. S. and Hindman, B. W. (2011).  Gyroscopic pumping in the solar near-surface shear layer.  Astrophys. J., 743(1):79.
  * Oyafuso et al., (2020) Oyafuso, F., Levin, S., Orton, G., Brown, S., Adumitroaie, V., Janssen, M., Wong, M. H., Fletcher, L. N., Steffes, P., Li, C., Gulkis, S., Atreya, S., S., M., and S., B. (2020).  Angular dependence and spatial distribution of Jupiter’s centimeter-wave thermal emission from Juno’s microwave radiometer.  Earth Planet. Sci. Lett., 7(11):e2020EA001254.
  * Porco et al., (2003) Porco, C. C., West, R. A., McEwen, A., Del Genio, A. D., Ingersoll, A. P., Thomas, P., Squyres, S., Dones, L., Murray, C. D., Johnson, T. V., Burns, J. A., Brahic, A., Neukum, G., Veverka, J., Barbara, J. M., Denk, T., Evans, M., Ferrier, J. J., Geissler, P., Helfenstein, P., Roatsch, T., Throop, H., Tiscareno, M., and Vasavada, A. R. (2003).  Cassini imaging of Jupiter’s atmosphere, satellites and rings.  Science, 299:1541–1547.
  * Potter et al., (2014) Potter, S. F., Vallis, G. K., and Mitchell, J. L. (2014).  Spontaneous superrotation and the role of Kelvin waves in an idealized dry GCM.  J. Atmos. Sci., 71(2):596–614.
  * Read et al., (2018) Read, P. L., Lewis, S. R., and Vallis, G. K. (2018).  Atmospheric dynamics of terrestrial planets.  Handbook of Exoplanets, 144:2537–2557.
  * Salyk et al., (2006) Salyk, C., Ingersoll, A. P., Lorre, J., Vasavada, A., and Del Genio, A. D. (2006).  Interaction between eddies and mean flow in Jupiter’s atmosphere: Analysis of Cassini imaging data.  Icarus, 185:430–442.
  * Schneider, (2006) Schneider, T. (2006).  The general circulation of the atmosphere.  Ann. Rev. Earth Plan. Sci., 34:655–688.
  * Schneider and Liu, (2009) Schneider, T. and Liu, J. (2009).  Formation of jets and equatorial superrotation on Jupiter.  J. Atmos. Sci., 66:579–601.
  * Showman and de Pater, (2005) Showman, A. P. and de Pater, I. (2005).  Dynamical implications of Jupiter’s tropospheric ammonia abundance.  Icarus, 174:192–204.
  * Taylor et al., (2004) Taylor, F. W., Atreya, S. K., Encrenaz, T. H., Hunten, D. M., Irwin, P. G., and Owen, T. C. (2004).  Jupiter: the planet, satellites and magnetosphere, chapter The composition of the atmosphere of Jupiter, pages 59–78.  Cambridge University Press.
  * Tollefson et al., (2017) Tollefson, J., Wong, M. H., de Pater, I., Simon, A. A., Orton, G. S., Rogers, J. H., Atreya, S. K., C., R. G., Januszewski, W., Morales-Juberías, R., and S., M. P. (2017).  Changes in Jupiter’s zonal wind profile preceding and during the Juno mission.  Icarus, 296:163–178.
  * Vallis, (2017) Vallis, G. K. (2017).  Atmospheric and Oceanic Fluid Dynamics.  pp. 770. Cambridge University Press., second edition.
  * Vasavada and Showman, (2005) Vasavada, A. R. and Showman, A. P. (2005).  Jovian atmospheric dynamics: An update after Galileo and Cassini.  Reports of Progress in Physics, 68:1935–1996.
  * Wicht and Gastine, (2020) Wicht, J. and Gastine, T. (2020).  Numerical simulations help revealing the dynamics underneath the clouds of Jupiter.  Nature Communications, 11(1):1–4.
  * Wicht et al., (2019) Wicht, J., Gastine, T., Duarte, L. D., and Dietrich, W. (2019).  Dynamo action of the zonal winds in Jupiter.  Astron. and Astrophys., 629:A125.
  * Young and Read, (2017) Young, R. M. and Read, P. L. (2017).  Forward and inverse kinetic energy cascades in Jupiter’s turbulent weather layer.  Nature Physics, 13(11):1135–1140.
  * Young et al., (2019) Young, R. M., Read, P. L., and Wang, Y. (2019).  Simulating Jupiter’s weather layer. Part I: Jet spin-up in a dry atmosphere.  Icarus, 326:225–252.

## Supporting Information

###  S1 Eddy momentum driven Ferrel-like cell

####  S1.1 Standard formulation

Using approximations similar to the commonly used formulation which describe
the terrestrial Ferrel-cell dynamics (Vallis,, 2017), the leading order zonal
mean zonal momentum equation may be written as

$\frac{\partial\bar{u}}{\partial t}+\frac{\partial}{\partial
y}\left(\overline{u^{{}^{\prime}}v^{{}^{\prime}}}\right)-f\bar{v}=-F_{{\rm
sink}},$ (S1)

where $F_{{\rm sink}}$ is a sink term. On Earth, the sink term represents a
surface drag in the Ekman layer and in the Jovian atmosphere, if the cells are
as deep as the jets (Kaspi et al.,, 2020), it might be Ohmic dissipation (Liu
et al.,, 2008; Liu and Schneider,, 2010; Kaspi et al.,, 2020) (i.e., $F_{{\rm
sink}}=\frac{1}{4\pi\bar{\rho}}\overline{\left(\nabla\times\mathbf{B}\right)\times\mathbf{B}},$
where ${\rm{\bf B}}$ is the 3D magnetic field and $\rho$ is density). Eq.
1111Note that equation cross-references without a S refer to equations in the
main text., adequate for the upper branch of the cells, results from applying
time-averaging away from the sink layer. This balance leads to meridional
velocities in the directions illustrated in Fig. 2b222Note that figure cross-
references without a S refer to figures in the main text.. Alternatively,
applying vertical integration cancels the Coriolis terms in the vertical
boundaries of the cell and the relation between the zonal jets, their eddy
source term and the sink term becomes

$\frac{\partial\bar{U}}{\partial t}=-\frac{\partial}{\partial
y}\left(\overline{U^{{}^{\prime}}V^{{}^{\prime}}}\right)-\hat{F}_{{\rm
sink}}.$ (S2)

where $U$ and $V$ are the vertically integrated velocities, and $\hat{F}_{{\rm
sink}}$ is the vertically integrated sink term. Therefore, converging
(diverging) eddy momentum fluxes transfer their momentum to eastward
(westward) jets, as can be seen in the Jovian atmosphere (Fig. 1b,d). In the
descriptions throughout this study, Cartesian approximations are used for the
sake of clarity, but the actual calculations were performed using the
spherical, more accurate, formulations.

####  S1.2 TEM formulation for Jupiter

The Transformed Eulerian Mean (TEM) equations, commonly invoked to quantify
Lagrangian mass-transport in Earth’s Ferrel cells, describe a circulation
driven by the diabatic heating term (Vallis,, 2017). The TEM formulation can
be derived from the momentum and thermodynamic equations, under Boussinesq and
quasi-geostrophic approximations (Vallis, (2017), ch. 10.3), resulting in

$\frac{\partial\bar{u}}{\partial t}-f\bar{v}^{*}=\nabla\cdot{\mathcal{{F}}},$
(S3)

and

$\frac{\partial\bar{b}}{\partial t}-\bar{w}^{*}N^{2}=\bar{S},$ (S4)

where $\bar{v}^{*}=\bar{v}-\frac{\partial}{\partial
z}\left(\frac{1}{N^{2}}\overline{v^{\prime}b^{\prime}}\right)$ and
$\bar{w}^{*}=\bar{w}+\frac{\partial}{\partial
y}\left(\frac{1}{N^{2}}\overline{v^{\prime}b^{\prime}}\right)$ are defined as
the “residual” mean meridional and vertical velocities, respectively, which
approximate mass-transport by both Eulerian mean velocities and eddy fluxes,
$N^{2}=\nicefrac{{\partial\bar{b_{0}}}}{{\partial z}}$ is the Brunt-Vaisala
frequency, $\bar{S}$ is the diabatic heating term, and
${\bf\mathcal{{F}}}=-\left(u^{\prime}v^{\prime}\right){{\rm\hat{j}}}+\left(\frac{f}{N^{2}}\overline{v^{\prime}b^{\prime}}\right){\rm{\hat{k}}}\,\,$
is the Elissan-Plam (EP) flux. $b_{0}$ and $b$ represent the mean (zonally-
and meridionally-averaged) and the deviation from the mean of the buoyancy
force and ${{\rm\hat{j}}}$ (${{\rm\hat{k}}}$) is a unit vectors in the
meridional (vertical) direction. On Earth, the diabatic heating term is
important (Vallis,, 2017), and therefore the residual meridional velocities
accurately represent the total meridional transport of mass in Earth’s
midlatitudes.

The midlatitude atmosphere on Earth is characterized by baroclinicity, and as
a result, the second term of the EP flux is substantial and plays a key role
in the resulting circulation. On Jupiter, the eddy fluxes beneath the cloud
level and the diabatic heating are yet to be measured. However, gravity-
measurement analysis implies that Jupiter’s jets are nearly barotropic in the
depth range relevant to this study (Galanti and Kaspi,, 2021), meaning that
the EP flux is dominated by the first term, and $\bar{w}^{*}$ is comparable to
$\bar{w}$ (Lee and Kaspi,, 2021). Therefore, under the barotropic limit, the
equations describing the Eulerian velocities in a Ferrel-like cell might also
represent the total mass transport in the Jovian cells.

####  S1.3 Correlation analysis

The picture illustrated in Fig. 2, relating the distribution of ammonia and
the zonal winds according to the Ferrel-like cells hypothesis, is tested
quantitatively in a correlation analysis exhibited in Fig. 3. The expected
relations between the zonal jets and the ammonia meridional gradients in the
NH are

$\bar{u}\propto\begin{cases}\begin{array}[]{c}-\partial_{y}m_{a}\\\
\partial_{y}m_{a}\end{array}&\begin{array}[]{c}1.5\leq p<6\,{\rm
bar}\quad\left({\rm channels}\,\,4-5\right),\\\ p\geq 6\,{\rm
bar}\quad\left({\rm channels\,}\,1-3\right),\,\,p<1.5\,{\rm
bar}\quad\left({\rm channel\,}\,6\right).\end{array}\end{cases}$ (S5)

Here, $\partial_{y}m_{a}$ is the latitudinal gradient of the ammonia
concentration ($m_{a}$) and $p$ is pressure. The channels refer to the six
frequencies of Juno’s MWR. In the SH, as the circulation is reversed, Eq. S5
flips signs. Eq. S5 captures also the case of cells with westward jets, as
both $\bar{u}$ and $\partial_{y}m_{a}$ change sign. A Pearson correlation
coefficient (${{S}}(\vartheta,{\rm ch})$) is calculated for each latitude and
MWR channel, and its value is represented by a color between blue,
representing a negative correlation, white, representing no correlation and
red, representing a positive correlation. The $T_{{\rm b}}$ data is measured
in a resolution of $\sim 0.6^{\circ}$ latitude. The data is interpolated such
that the grid size is $0.1^{\circ}$ latitude, and the correlation for each
point $\vartheta_{i}$ is calculated along a span
{$\vartheta_{i}-2^{\circ}$,$\vartheta_{i}+2^{\circ}$}. This choice of a
$4^{\circ}$ latitudinal bin allows having enough data points for the
statistical value of the correlation (more than 6 data pairs), and ensures the
local nature of the results. The correlations on the MWR data are calculated
between the following trends. In Fig. 3a, channels 1-3 and 6 and in Fig. 3c
channels 4-5 the color represents the correlation
$\bar{u}\propto\partial_{y}T_{{\rm b}}$. In Fig. 3a, channels 4-5 and in Fig.
3c channels 1-3 and 6 the color represents the correlation
$\bar{u}\propto-\partial_{y}T_{{\rm b}}$. Note that anomalies of brightness
temperature and ammonia abundance are inversely proportional (Li et al.,,
2017). In Fig. 3b the correlations at all channels are calculated according to
$\bar{u}\propto-T_{{\rm b}}$. Here, $T_{{\rm b}}$ is the Nadir component of
the brightness temperate [${}^{\circ}{\rm K}$] (other emission angles were not
included in the analysis), averaged over nine Juno orbits (PJs
$1,3,4,5,6,7,8,9$ and $12$) (Oyafuso et al.,, 2020). For further discussion
regarding limb-darkening $T_{{\rm b}}$ we point the readers to Fletcher et
al., (2021). $\bar{u}$ is Jupiter’s zonally-averaged zonal wind [${\rm
m\,\,s^{-1}}$] measured by the Hubble space telescope during Juno’s third
perijove (Tollefson et al.,, 2017), projected barotropically along the axis of
rotation (Galanti and Kaspi,, 2021; Galanti et al.,, 2021). The ammonia
distribution by JIRAM is estimated at a depth of $\sim 6$ bar, which is the
depth where a local minimum appears in $M_{{\rm a}}$ (Fig. 2a). Arbitrarily,
the correlation is performed according to $p>6$ bar in Eq. 2 (Fig. 2d) and the
overall positive result points that indeed the depth level of JIRAM
measurements should be deeper than the local minimum of $M_{{\rm a}}$. The
ammonia estimates from JIRAM are measured in a resolution of $1^{\circ}$
latitude. Similar to the $T_{{\rm b}}$ correlation analysis, the data is
interpolated, and the correlation is performed on a $4^{\circ}$ latitude bin
for consistency. In Fig. 3d (f) the color represents the correlation
$\bar{u}\propto-\partial_{y}m_{a}$ ($\bar{u}\propto\partial_{y}m_{a}$) and in
Fig. 3e the color represents the correlation $\bar{u}\propto m_{a}$. Finally,
the eddy momentum flux convergence is measured in a resolution of $1^{\circ}$
latitude, and the correlation is performed using the same latitudinal bin of
$4^{\circ}$. In Fig. 3g (h) the color represents the correlation
$\bar{u}\propto-\partial_{y}\left(\overline{u^{\prime}v^{\prime}}\right)$
($\bar{u}\propto\partial_{y}\left(\overline{u^{\prime}v^{\prime}}\right)$).

We also examine the correlation between the zonal wind and the lightning
gradient. We find a good match in the northern hemisphere and a weak negative
correlation in the southern hemisphere, where Juno is much less sensitive to
lightnings (Fig. S1). Note that the correlation between lightnings and the
Ferrel cells is less indicative, as we should only examine the correlation in
the rising branch of the cells. Therefore, the correlation values away from
the rising branch should be regarded with caution.

Figure S1:  Correlation coefficients (upper panels) between the zonal wind
(middle panels) and the lightning meridional gradient (lower panels) in the
midlatitudes.

####  S1.4 Advection-relaxation model

#####  S1.4.1 $T_{{\rm b}}$ as an indicator for ammonia

Figure S2: (a) Ammonia estimate (Li et al.,, 2017) from Juno MWR PJ1 data. (b)
$m_{a}$ calculated from $T_{{\rm b}}$ (Eq. S8). This field is used as a
benchmark ($m_{a}^{({\rm data})}$) for the model results. (c) the difference
between panel a and panel b.

As ammonia estimates by Juno’s MWR for the high midlatitudes are not yet
available, we express the $T_{{\rm b}}$ measurements as ammonia in order to
examine the model results. For that, we define a reconstructed ammonia
distribution from Juno MWR data ($m_{a}^{({\rm data})}$), used as a benchmark
for the advection-relaxation model, constructed by the mean ammonia calculated
from MWR measurements of PJ1 ($M_{{\rm a}}$, Fig. 2a) (Li et al.,, 2017) and
$T_{{\rm b}}$ measurements averaged over multiple Juno orbits (PJs 1-12)
(Oyafuso et al.,, 2020). The standard deviation between the perijoves is
computed as a function of latitude to validate that the latitudinal variations
appearing in the $T_{{\rm b}}$ data are physical (Fig. S3 for the midlatitudes
and Fig. S4 for the equatorial region). We estimate the $T_{{\rm b}}$
anomalies ($T_{{\rm b}}^{\prime}$) by removing the cross-channel average
($T_{{\rm b},{\rm mean}}({\rm ch})$) from each MWR channel, and then decompose
$T_{{\rm b}}^{\prime}$ into Legendre polynomials and reconstruct the anomalies
without the low polynomials to remove large scale variations. These
variations, representing equator to pole radiation differences, are not
related to ammonia variations by meridional cells (Oyafuso et al.,, 2020). The
reconstructed ammonia is then

$T_{{\rm b}}^{\prime}(\vartheta,{\rm ch})=T_{{\rm b}}(\vartheta,{\rm
ch})-T_{{\rm b},{\rm mean}}({\rm ch})\cong\sum_{i=1}^{N}A_{i}(\vartheta,{\rm
ch})P_{i}\left(\sin\vartheta\right),$ (S6)

$T_{{\rm b,rec}}^{\prime}(\vartheta,{\rm
ch})=\sum_{i=30}^{N}A_{i}(\vartheta,{\rm ch})P_{i}\left(\sin\vartheta\right),$
(S7)

$m_{a}^{({\rm data})}(\vartheta,{\rm ch})=M_{{\rm a}}({\rm ch})-K({\rm
ch})T_{{\rm b,rec}}^{\prime}(\vartheta,{\rm ch}),$ (S8)

where $P_{i}$ are the Legendre polynomials, $A_{i}$ are the associated
coefficients, $N=200$ is the number of polynomials used and $K$ $[{\rm
ppm}\cdot{\rm degrees}^{-1}]$ is a depth-dependent ’key’, optimized at each
depth using Matlab’s ’fmincon’ to best fit the estimated ammonia distribution
from PJ1 (Li et al.,, 2017). Note that $T_{{\rm b}}$ is available between
latitudes $90^{\circ}{\rm S}$ to $90^{\circ}{\rm N}$, therefore $P_{i=30}$ is
equivalent to a wavelength of approximately $12^{\circ}$ latitude. As the jets
widths are not larger the $8^{\circ}$ in latitude, this truncation removes
variations that are not due to the existence of meridional Ferrel-like cells.
The overall structure of $m_{a}^{({\rm data})}$ is very similar to the ammonia
map from PJ1 (Li et al.,, 2017), while the meridional anomalies now represent
well PJs 1-12 (Fig. S2). Finally, the resulting $K$ ranges between $5$ at $1$
bar to $0$ at depth and is used to estimate $m_{a}^{({\rm data})}$ at
latitudes $60^{\circ}{\rm S}$ to $60^{\circ}{\rm N}$. The profile is
interpolated between $1$ and $240$ bar, according to the relevant pressure
levels of each channel (Fig. 4d). These levels are estimated according to the
peak of the contribution function of each MWR channel, to give that channels
$\left\\{{\rm 1,2,3,4,5,6}\right\\}$ correspond to pressure levels of
$\left\\{240,30,9,3,1.5,0.7\right\\}$ bar (Janssen et al.,, 2017; Bolton et
al.,, 2017).

Figure S3:  The standard deviation (STD) of the nadir $T_{{\rm b}}$ in the
midlatitudes computed from nine different perijoves through PJs 1-12 (Oyafuso
et al.,, 2020), in the six MWR channels. The STD is computed after the trend
removal for each perijove, as detailed in eq. S7. For channels 1-5, the STD
values are smaller than the variation seen in the data (Fig. 1h). Note that
although the STD values for channel 6 are higher, at some latitudes, than the
mean latitudinal variation (of the same channel), this channel senses at
altitudes which are above the cells identified in this study. Figure S4:  The
equatorial standard deviation of the nadir $T_{{\rm b}}$ computed from nine
different perijoves through PJs 1-12 (Oyafuso et al.,, 2020), in the six MWR
channels. The STD is computed on the original nadir coefficients.

#####  S1.4.2 Cells construction and parameterization

To describe the meridional cells in the simplest manner, we parameterize each
cell (indexed $k$) with an ellipse, according to a parameter $l_{k}$ as
follows

$l_{k}=\sqrt{\frac{\left(d-a\right)^{2}}{a^{2}}+\frac{\left(\vartheta-\vartheta_{k}\right)^{2}}{b_{k}^{2}}}.$
(S9)

Here, $d$ is defined as downward distance from the cloud level, and $a$ and
$b_{k}$ are the vertical and meridional extents of the cell, respectively.
$\vartheta_{k}$ is the latitude of the center of cell $k$ (black dots in Fig.
S5a,b). $b_{k}$ is set according to half the width of cell $k$ (the distance
between a black line and a black dot in Fig. S5a,b). The outline of cell $k$
(representing the path of the peak tangential velocity along the cell) is
thereby defined by $l_{k}=1$. For simplicity, the velocities in a cell are
defined using a normal distribution according to

$\displaystyle v_{k}=$ $\displaystyle
V_{k}\exp\left[-\frac{1}{2}\left(\frac{l_{k}-1}{\widetilde{\sigma}}\right)^{2}\right]\sin\phi,$
(S10) $\displaystyle w_{k}=$ $\displaystyle
V_{k}\exp\left[-\frac{1}{2}\left(\frac{l_{k}-1}{\widetilde{\sigma}}\right)^{2}\right]\cos\phi,$
$\displaystyle\widetilde{\sigma}=$
$\displaystyle\sigma\left(\cos^{2}\vartheta+\frac{b_{k}}{a}\sin^{2}\vartheta\right),$

where $\phi=\arctan\left[\frac{d-a}{R_{{\rm
J}}\sin\left(\vartheta-\vartheta_{k}\right)}\right]$, $R_{{\rm J}}$ is
Jupiter’s radius and $\sigma$ is a parameter for the broadness of a cell’s
branch. $V_{k}$ represents the relative strength (velocity) of cell $k$,
parameterized according to the square root of the averaged eddy momentum flux
convergence within the cell (Fig. S5b), and its sign represents the cell’s
direction (clockwise/counter-clockwise), set according to the zonal wind sign
at the center of the cell (Fig. S5a).

Figure S5:  (a) Vorticity (purple) and zonal wind (dashed, orange). Blue (red)
triangles represent vorticity peaks
$\left(\nicefrac{{\partial\zeta}}{{\partial y}}=0\right)$, where the vertical
branches of the cells drive upward (downward) motion. The cells’ centers are
positioned between vorticity peaks (black dots), and the cells’ extents
($2b_{k}$) are the distances between pairs of black lines. (b) Eddy momentum
flux convergence (blue) and the cells’ location as in panel a. Cells’ relative
strength is set by the averaged value of eddy flux convergence within the cell
(light blue circles).

#####  S1.4.3 Optimization and numerical solution

To solve for $m_{{\rm a}}$, Eq. 3 is discretized using finite differences as

$\bar{w}_{i,j}\frac{m_{i-1,j}-m_{i+1,j}}{2dr}+\bar{v}_{i,j}\frac{m_{i,j+1}-m_{i,j-1}}{2R_{{\rm
J}}d\vartheta}=-G_{i}\left(m_{i,j}-M_{i}\right),$ (S11)

where the “a” subscript of $m$ and $M$ was removed for clarity. Here $i$, and
$j$ are indices for the grid points in the r and $\vartheta$ directions,
respectively. $dr$ and $d\vartheta$ are the distances between adjacent points
in each direction. Eq. S11 constitutes one of $n^{2}$ equations for $n^{2}$
variables, where $n$ is the resolution of the grid in each direction. Eq. S11
is rearranged in a matrix form as $Ax=b$ such that $m_{{\rm a}}$ can be
calculated from $A^{-1}b$.

The parameters $G_{i}$, $\sigma$ and $a$ are unknowns. For this, the Matlab
optimization function ’fmincon’ is used for deciding $G_{i}$, $\sigma$ and $a$
to best reproduce the $m_{a}^{({\rm data})}$ map. The cost function

$f(G,\sigma,a)=\underset{i,j}{\sum}\left(\left|m_{i,j}^{({\rm
model})}-m_{i,j}^{({\rm data})}\right|\right)^{2}$ (S12)

is the measure used to find the optimal parameters. The resulting value of $G$
is shown in Fig. S6 for the case of $G$ that is varying with depth, and for
comparison, the solution with a constant $G$ is shown as well. The value of
$\sigma$ was found to be $\sim 0.85$ in both cases. The parameter $a$ was
found to be $\sim 1600$ bar.

#####  S1.4.4 Robustness analysis

Figure S6:  Comparison between model results with and without vertical
variation in $G.$ (a) The vertical variation of the normalized source term $G$
used for Fig. 4c. b and e, $m_{a}^{({\rm data})}$ anomalies [ppm] in the SH.
(c) The ammonia anomalies [ppm] map produced by the advection-relaxation model
with the source term from panel a. (d) Constant normalized source term $G$ in
Eq. 3. (f) The ammonia anomalies [ppm] map produced by the advection-
relaxation model with the source tern from panel d. In panels b, c, e and f
the vertical mean profile $M_{{\rm a}}$ is removed from the ammonia map
$m_{{\rm a}}$ and dashed red and green lines are the upward and downward
branches of the cells, respectively. The comparison reveals that although the
solution with varying source term (a) results in a better model solution (c)
compared to the measurements (b and d), the essence of the anomalies (f) is
well captured with a constant source term (d). Figure S7: Model robustness
analysis. An example for a model run without optimization. The optimized
variables in this run are set manually according to physical considerations.
The depth of the cells ($a$) is set to $1,500$ km such that the cells extend
$3,000$ km in accordance with gravity measurements for the depth of the zonal
jets (Kaspi et al.,, 2018). The relaxation constant $G$ is set (without height
dependence) from scaling argument by equating the relaxation term in equation
3 to the vertical advection term, leading to $G={\rm max}(w)/a$. The value of
$\sigma$ is set to $0.5$. (a) The normalized vertical velocity as a function
of latitude and pressure. (b) The ammonia anomalies map that was reconstructed
from Juno’s MWR measurements. (c) The map of ammonia anomalies resultant from
the degenerated model. It can be seen that the optimization procedure doesn’t
change the nature of the results which is robust. The structure of the
anomalies stays largely the same both in this figure and in Fig. 4, and it
stems mostly from the latitudinal structure of the wind and the vertical
stratification of the ammonia, both being derived from observations.

To validate that the model results are robust and not sensitive to the
specific parameters found by the optimization analysis, $m_{i,j}^{({\rm
model})}$ was also solved from equation 3 using a chosen set of parameters
instead of the optimized set. The depth of the cells was chosen to be $3000$
km ($a=1500\,{\rm km}$), in accordance with the depth of the jets that was
estimated from gravity measurements (Kaspi et al.,, 2018). The width of the
cells was chosen such that one standard deviation covers half a cell’s width
($\sigma=0.5$). $G$ is degenerated into a constant (no dependence in $z$) and
is set, from a scaling argument, as $G={\rm max}(w)/a\,{\rm[s^{-1}]}$. As seen
in the results (Fig. S7), the modeled ammonia anomalies map still predicts the
data convincingly. The main difference is the depth where the sign of the
ammonia anomalies flip (between $3$ and $6$ bar), which is now only controlled
by the input $M_{{\rm a}}$ (Li et al.,, 2017), and could not be ’corrected’ by
a depth-dependent relaxation time scale. It can be seen that the essence of
the circulation cells is still very apparent in the results.

###  S2 Equatorial region analysis

####  S2.1 Estimate for the extent of the equatorial region

The tangent cylinder is the projection (along the axis of rotation) of the
planet’s solid-body rotating core on the outer shell. The equatorial latitudes
lie outside of the tangent cylinder, and is thereby characterized by a
different dynamical regime than that of the midlatitudes. To separate
quantitatively the midlatitudes, positioned within the tangent cylinder, from
the equatorial region, it is required to know the depth of the atmosphere
($D_{{\rm atm}}$) and the radius of the planet ($R$). The latitudes of the
cylinder’s edge ($\alpha$) can then be derived from geometrical considerations
as

$\alpha=\arccos\left(1-\frac{D_{{\rm atm}}}{R}\right).$ (S13)

Gravity analysis reveals that Jupiter’s atmosphere is approximately
$3000\,\,{\rm km}$ deep (Guillot et al.,, 2018; Kaspi et al.,, 2018).
Substituting $D_{{\rm atm}}=3000\,\,{\rm km}$ and $R=R_{{\rm
J}}=70,000\,\,{\rm km}$ in Eq. S13 gives $\alpha=\pm 16.8^{\circ}$. This means
that fluid columns parallel to the axis of rotation in the latitude range
$-17^{\circ}\leq\vartheta\leq 17^{\circ}$, can theoretically extend
uninterruptedly between the hemispheres.

####  S2.2 Theory for the leading balance in the Jovian equatorial region

Starting from the primitive equations(Vallis,, 2017), the continuity and zonal
momentum equations in spherical coordinates are

$\frac{\partial\rho}{\partial t}+\rho\left(\frac{1}{R_{{\rm
J}}\cos\vartheta}\frac{\partial u}{\partial\lambda}+\frac{1}{R_{{\rm
J}}\cos\vartheta}\frac{\partial}{\partial\vartheta}\left(v\cos\vartheta\right)+\frac{\partial
w}{\partial r}\right)+\frac{u}{R_{{\rm
J}}\cos\vartheta}\frac{\partial\rho}{\partial\lambda}+\frac{v}{R_{{\rm
J}}}\frac{\partial\rho}{\partial\vartheta}+w\frac{\partial\rho}{\partial
r}=0,$ (S14)

and

$\frac{\partial u}{\partial t}+\frac{u}{R_{{\rm
J}}\cos\vartheta}\frac{\partial u}{\partial\lambda}+\frac{v}{R_{{\rm
J}}}\frac{\partial u}{\partial\vartheta}+w\frac{\partial u}{\partial
r}-2\Omega\left(\sin\vartheta\,v-\cos\vartheta\,w\right)-\frac{uv}{R_{{\rm
J}}}\tan\vartheta=-\frac{1}{R_{{\rm J}}\rho\cos\vartheta}\frac{\partial
p}{\partial\lambda},$ (S15)

respectively, where $\lambda$ is longitude. The mean density ($\rho_{m}$) is
assumed to only change with $r$, and anomalies from it are assumed to be much
smaller than the mean value. In addition, the meridional derivatives and the
meridional velocity are assumed very small near the equator, relative to the
other terms. This assumption is based on the symmetry around the equator.
Assuming also a steady state, Eq. S14 and Eq. S15, evaluated on the equatorial
plane ($\vartheta=0^{\circ}$), are then

$\rho_{m}\left(\frac{1}{R_{{\rm J}}}\frac{\partial
u}{\partial\lambda}+\frac{\partial w}{\partial
r}\right)+w\frac{\partial\rho_{m}}{\partial r}=0,$ (S16)

and

$\frac{u}{R_{{\rm J}}}\frac{\partial u}{\partial\lambda}+w\frac{\partial
u}{\partial r}+2\Omega w=-\frac{1}{R_{{\rm J}}\rho_{m}}\frac{\partial
p}{\partial\lambda}.$ (S17)

Eq. S17 can equivalently be represented as

$\frac{1}{R_{{\rm J}}}\frac{\partial
u^{2}}{\partial\lambda}+\frac{1}{\rho_{m}}\frac{\partial\left(wu\rho_{m}\right)}{\partial
r}-\frac{u}{\rho_{m}}\left(\rho_{m}\frac{1}{R_{{\rm J}}}\frac{\partial
u}{\partial\lambda}+\rho_{m}\frac{\partial w}{\partial
r}+w\frac{\partial\rho_{m}}{\partial r}\right)+2\Omega w=-\frac{1}{R_{{\rm
J}}\rho_{m}}\frac{\partial p}{\partial\lambda},$ (S18)

where the third term vanishes according to Eq. S16. Next, the velocities are
decomposed via Reynolds decomposition
($u=u^{\prime}+\overline{u},\,\,\,\,w=w^{\prime}+\overline{w}$), such that the
zonal mean of Eq. S18 gives

$-\frac{\partial\left(\overline{w^{\prime}u^{\prime}}\rho_{m}\right)}{\partial
r}=\overline{w}\frac{\partial\left(\overline{u}\rho_{m}\right)}{\partial
r}+2\rho_{m}\Omega\overline{w}+\rho_{m}\overline{u}\frac{\partial\overline{w}}{\partial
r}.$ (S19)

This shows that the momentum originating from the eddy momentum flux
convergence
$\left(\frac{\partial\left(\overline{w^{\prime}u^{\prime}}\rho_{m}\right)}{\partial
r}<0\right)$, which drives the equatorial superrotation (Kaspi et al.,, 2009),
is divided between the growing equatorial super-rotating jet
$\left(\frac{\partial\left(\overline{u}\rho_{m}\right)}{\partial
r}>0,\,\overline{u}>0\right)$, the Coriolis force and another residual term.
The growing equatorial jet has been shown in many numerical simulations of
superrotation (Heimpel et al.,, 2005; Kaspi et al.,, 2009; Gastine et al.,,
2014). It is a good assumption that each of the terms on the right side is of
smaller magnitude than the source term on the left side. Finally, rearranging
Eq. S19 gives

$\overline{w}=-\frac{1}{\underset{>0}{\underbrace{\frac{\partial\left(\overline{u}\rho_{m}\right)}{\partial
r}}}+\underset{>0}{\underbrace{2\rho_{m}\Omega}}}\underset{<0}{\underbrace{\left(\frac{\partial\left(\overline{w^{\prime}u^{\prime}}\rho_{m}\right)}{\partial
r}+\rho_{m}\overline{u}\frac{\partial\overline{w}}{\partial r}\right)}}>0.$
(S20)

To further simplify, Eq. S20 can be shown for the case of small Rossby number:

$\overline{w}=-\frac{1}{2\rho_{m}\Omega}\frac{\partial\left(\overline{w^{\prime}u^{\prime}}\rho_{m}\right)}{\partial
r}>0.$ (S21)

This implies that the mean upwelling ($\overline{w}$) correlates with eddy
momentum flux convergence, and therefore with the equatorial superrotating jet
at the equatorial region. Since the superrotating jet is supposed to be driven
by angular momentum fluxes in the direction perpendicular to the rotation axis
(Heimpel et al.,, 2005; Kaspi et al.,, 2009; Schneider and Liu,, 2009),
converging in the equatorial region, Eq. S21 would take the more general form

$\overline{w_{\perp}}=-\frac{1}{2\rho_{m}\Omega}\partial_{\perp}\left(\overline{w_{\perp}^{\prime}u^{\prime}}\rho_{m}\right)>0,$
(S22)

where $w_{\perp}$ and $\partial_{\perp}$ are the velocity and the gradient in
the direction perpendicular to the axis of rotation, i.e.,
$w_{\perp}=w\cos\vartheta+v\sin\vartheta$ and
$\partial_{\perp}=\left(\cos\vartheta\right)\partial_{r}+\left(r^{-1}\sin\vartheta\right)\partial_{\vartheta}$.

Figure S8:  Meridional averaged ammonia values at different latitudinal
regions. The lines are calculated according to the inferred ammonia map from
PJ1 (Li et al.,, 2017). Each line is averaged at the latitudinal range
described in the legend.