# Ammonia and Water Abundance derived from Juno MWR and VLA Observations of
Jupiter

###### Abstract

The vertical distribution and evolution of trace gases in planetary
atmospheres, a key aspect for understanding the global circulation on gas
giants, can be obtained with observation of the atmosphere’s thermal emission.
Inverting radio observations to recover atmospheric structure, however, is
non-trivial, and the solutions are degenerate. We propose a modeling framework
to prescribe a vertical distribution of trace gases that combines a thermo-
chemical equilibrium model with a robust handling of uncertainties. To this
means we retrieve nadir brightness temperatures and limb-darkening parameters,
together with their uncertainties, from the Juno Microwave Radiometer (MWR).
We then apply this framework to MWR observations during Juno’s first year of
operation (Perijove passes 1 - 12) and to longitudinally-averaged scans taken
with the upgraded Very Large Array (VLA) (dePater2016; dePater2019). We use
the model to constrain the atmosphere between -60∘ and 60∘ latitude and down
to 100 bar based on MWR and VLA data and find an excellent agreement between
the two measurement sets. We constrain the deep ammonia abundance to be
$366.5^{+32.8}_{-19.9}$ ppm ($2.47^{+0.22}_{-0.13}\times$ solar abundance),
and find a depletion of ammonia down to a depth of $\sim$ 20 bar, which
supports processes that deplete the atmosphere below the ammonia and water
cloud layers. Lastly, we find a weak evidence for an enhanced water abundance
in the belts compared to the equatorial zone, however, any strong conclusions
about the water abundance on Jupiter are pre-mature and require further
labratory experiments to understand opacities of trace gases at high pressures
and temperatures.

Planetary atmospheres, Planetary climates, Atmospheric composition, Water
vapor

††journal: PSJ

##  1 Introduction

Some of Jupiter’s biggest mysteries remain unsolved despite decades of
observations. Questions about its formation environment, the source of energy
for its atmosphere, and the distribution of trace gases within said atmosphere
remain as open-ended as ever. The vastness of the planet and its atmosphere,
coupled with the short time scales on which the atmosphere evolves, implies
that a global understanding of the planet requires coordinated observations
across many different frequencies. In the early days of observational
astronomy, optical images were used to infer the dynamics of the atmosphere by
tracking the motion of cloud tops (peek_1959; Ingersoll1981). Various
techniques were used to correlate images and find the dominant zonal component
of the global circulation. The meridional component of the atmosphere is much
smaller in magnitude, making it more difficult to observe. Only recent
advances in observational capabilities and computing techniques have allowed
for direct measurement of the north-south component of the wind by tracking
small scale cloud features (Asay-davis2011).

A key feature in expanding this limited two-dimensional view of the atmosphere
is observations of the planet’s thermal emission at infrared (Gierasch1986)
and radio wavelengths. At radio wavelengths the impact of clouds is
sufficiently small such that we can still draw conclusions about the
atmosphere below, allowing radio observers to spectulate as early as the
1960’s about the presence of micowave absorbers in the Jovian atmosphere below
the visible cloud tops (Thornton1963; Welch1966; Wrixton1971; Berge1976).
Given the right combination of frequencies, we can use radio observations to
construct a full three-dimensional view of the planet and reveal global
circulation patterns.

Studying this pattern on the planet requires observations at resolutions high
enough to distinguish between the various climatological zones, i.e., the
visibly white zones and brown belts on Jupiter. Ground-based observatories use
the ability to detect phase at radio wavelengths and have built a complex
interferometers, with baselines on order of 10s of kilometers, such as the
Very Large Array (VLA) and the Atacama Large (sub)Millimeter Array (ALMA),
both providing excellent spatial resolution from Earth. The NASA Juno mission
opted for an alternative solution, building an orbiter that approaches Jupiter
as close as 3000 kilometer. Intriguingly, despite the six orders of magnitude
difference in distance from the target between the two observatories, the
final resolution is comparable, as the baseline of the VLA antennas are six
orders of magnitude larger than the antennas’ diameters on the orbiter. Though
both approaches have orthogonal advantages and disadvantages for atmospheric
observations, the combination provides a very powerful tool for studying the
atmosphere of Jupiter: the ground-based observations have excellent spectral
and spatial (both in latitude and longitude) resolution, but are limited by
the presence of synchrotron radiation at the longest wavelengths, while the
proximity of Juno to the planet allows for observations that avoid synchrotron
radiation; however, data are taken along a limited north-south track so there
is no global context for these observations.

The model we develop in this work is applicable to both sets of measurements,
and provides a reference frame for fitting and interpreting radio
observations. In Section 2, we present our data reduction pipeline developed
for the Juno data, followed by our description of the perturbed thermo-
chemical equilibrium model for the atmospheric structure in Section 3. Based
on a parameterized distribution of trace gases, we can compute the brightness
temperature using our radiative transfer code, as described in Section 4. In
Section 5, we present our results that minimize the difference between the
observations and the model atmosphere for the mean atmosphere based on the
first year of Juno observation, and compare our results to the only other
published result for perijove (PJ) 1. Lastly, we touch upon the most important
interpretations in Section 6, and comment on the state of the water abundance.
We finish the paper with our conclusions and recommendations for how to
advance the state of the art.

##  2 Data

###  2.1 Sources of Data

Since Juno’s insertion into Jupiter’s orbit on July 4th, 2016, the spacecraft
has remained in a 53 day, highly eccentric, polar orbit, in which the
spacecraft spends most of its time outside of the high radiation environment
around the planet. Among the suite of instruments, there are several with the
primary goal of studying the planet’s atmosphere, including a microwave
radiometer (MWR; Janssen2017). This instrument couples six radio receivers to
a set of antennas (an array of patch antennas and corrugated horns) to receive
radio waves with minimal side lobe contributions at six distinct frequencies,
from here on referred to as Channel 1 through 6 (see Table 1 for details).

Table 1: Radiometer and antenna characteristics of Juno’s MWR. Channel | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
Frequency [GHz] | 0.6 | 1.278 | 2.597 | 5.215 | 10.004 | 22
Bandwidth [MHz] | 21 | 43.75 | 84.5 | 169 | 325 | 770
Beamwidth [deg] | 20.6 | 20.6 | 12.1 | 12.1 | 12.0 | 10.8

Jupiter’s emission at radio wavelengths consists of two components: a thermal
component (i.e., the blackbody radiation of the planet’s atmosphere) and non-
thermal emission originating from synchrotron radiation by electrons in the
radiation belts or sferic emission by lightning on the planet (e.g.,
Berge1976; dePater1990; Brown2018). As seen from Earth, part of the
synchrotron radiation is spatially co-located with the atmospheres’ thermal
emission (see, e.g., dePater1986), which makes observations increasingly
difficult below 7 GHz (the frequency at which synchrotron radiation becomes
dominant (carr1983phenomenology)). Juno circumvents this superposition of
emitters by dipping between the planet and the radiation belts during PJs,
which brings the spacecraft as close as $\sim$3000 km above the cloud tops of
the planet.

The spacecraft completes one rotation every 30 seconds, in which time its
beams sweep across the sky (sky coordinates $\theta_{sky},\phi_{sky}$) and
obtain 300 beam-convolved measurements (i.e. the brightness temperature
$T_{b}$) of the sky’s brightness ($T_{sky}$) at positions expressed in
spherical coordinates in the spacecraft frame ($\theta_{sc},\phi_{sc}$)
(Janssen2017):

$T_{b}=\int_{0}^{\pi}\int_{0}^{2\pi}{T_{sky}(\theta_{sky},\phi_{sky})G(\theta_{sky}-\theta_{sc},\phi_{sky}-\phi_{sc})}sin(\theta_{sc})d\theta_{sc}d\phi_{sc}$
(1)

The radiation pattern (or ‘gain’) of the radio antenna $G(\theta,\phi)$ is
designed so that 99% of the emission falls within 2.5 times the half power
beam width (HPBW) of the Gaussian main beam, limiting the amount of radiation
received through the side lobes. The intersection of the main beam with the
planet determines the resolution of the measurements, which is highest at
nadir viewing, and decreases with increasing slant angle, distance from the
planet, and wavelength. For ideal operation of the MWR, the antennas scan the
planet in a north-south direction to obtain multiple measurements per latitude
at varying emission angle (the angle between the antenna’s line of sight and
the surface normal). While measurements at nadir probe deepest into the
planet, observations at decreasing emission angles scan the planet higher up
in the atmosphere. Juno’s other science priorities, such as gravity science,
require communication antennas to be aligned with Earth so that the satellite
does not align with the spin axis of the planet, which deteriorates the
emission angle coverage, even to the degree where we cannot retrieve any
observables (such as PJ 10-11). The most relevant information for the MWR
observation campaign are listed in Table 2.

The radiation pattern for the six antennas, the spacecraft geometry during
PJs, and the corresponding antenna temperature measurements are documented in
the Planetary Data System***Dataversion 3, https://pds.nasa.gov/ and are
utilized in the following analysis.

Table 2: MWR Orbits used in our work

PJ number | Date | Min. dist (km) | PJ [lat,lon] | Lon range
---|---|---|---|---
1 | 2016-08-27 | 4138 | 4.3, 95.8 | -115, -75
3 | 2016-12-16 | 4101 | 6.5, - 5.5 | -18, 5
4 | 2017-02-02 | 4234 | 7.5, 85.1 | 64, 107
5 | 2017-03-27 | 3314 | 8.6, 175 1 | 64, 107
6 | 2017-05-19 | 3383 | 9.7, -140 | 165, -177
7 | 2017-07-11 | 3359 | 10.8, -49.0 | -156, 114
8 | 2017-09-01 | 3332 | 11.8, 41.9 | -73, -21
9 | 2017-10-24 | 3839 | 13., 132. | 30, 73
12 | 2018-05-24 | 3201 | 15.9, 115 | -144, -107

Overview of MRW’s orbits which yield data useful for atmospheric science. The
columns correspond to the PJ number, date of the flyby, minimum distance from
the planet, sub satellite location during closest encounter, and the longitude
range covered during the flyby.

In addition to the Juno data, we also incorporate observations between 3 and
25 GHz from the Very Large Array (VLA) into this analysis. These were
originally published in dePater2016 and dePater2019, and we refer the reader
to these papers for details on the data reduction and calibration.

###  2.2 Data reduction of the Juno observations

The two quantities needed from the Juno data for our atmospheric retrieval
(see Section 5) are the nadir brightness temperature at the observing
frequency, $T_{n}$, and the limb-darkening coefficient $p$, which numerically
represents how the brightness temperature decreases at increasing emission
angle $\theta$ ($\mu\equiv\cos\theta$) as we probe higher up in the atmosphere
compared to nadir observations:

$T_{b}(\mu)=T_{n}{\mu}^{p}$ (2)

with $T_{b}(\mu)$ the brightness temperature of Jupiter at a given frequency
convolved with the gain pattern of the antennas (Equation 1).

For a full treatment of the problem, the observations can be deconvolved based
on an exact knowledge of the antenna pattern and using several overlapping
observations (Oyafuso2020). We choose an alternative approach, and instead
correct the observations for the viewing geometry. Since we are interested in
obtaining the variation in brightness temperature with emission angle, we must
account for the fact that the beam captures a wide range of emission angles.
Based on the shape of the planet and the relative position of the antenna
beam, we derive a map of emission angle for each observation (see inset b) in
Figure 1, and Appendix A). Knowing the detailed shape of the antenna pattern,
we can convolve the two quantities to obtain a beam-averaged temperature
${T_{b}(\bar{\mu})}$. Assuming thatbrightness temperature variations the
observable quantities vary across the beam are small compared to variations in
brightness due to changes in emission angle, we obtain the final brightness
temperature of the observation as the mean brightness temperature across the
beam ${T_{b}(\bar{\mu})}$ at the mean emission angle across the beam
$\bar{\mu}$:

$T_{b}(\bar{\mu})=T_{n}\bar{\mu}^{p}$ (3)

$\bar{\mu}=\frac{\iint_{0}^{\theta_{sc}}G(\theta_{sc},\phi_{sc})\mu(\theta_{sc},\phi_{sc})}{\iint_{0}^{\theta_{sc}}G(\theta_{sc},\phi_{sc})}$
(4)

We note here that other parameterizations for the limb-darkening can be used
(such as the two term expansion by Li2017)); however, the exact form does not
influence the atmospheric retrieval. All we care about for the retrieval is by
how much the brightness temperature has dropped at an emission angle of 45∘
for our model.

We set the limits of the beam convolution equal to 2.5 times the HPBW of the
beam, corresponding to the angle within 99% of the total emission is received
(Janssen2017). In case a beam element falls off the planet, we chose to
discard that measurement due to drastic brightness contrast between the planet
and either the background sky or synchrotron radiation.

We use the boresight emission intercept with the planet to bin the
measurements in 1 degree latitude bins. Inset a) in Figure 1 shows an example
of this, where we plot all 60 observations where the beam-center falls within
the selected region. The resolution of the actual observations is set by the
size of the footprints, which varies as the spacecraft is moving constantly
while it scans the planet. The 1-degree bin size therefore does not correspond
to the resolution of the data, but this binning scheme still captures
variations with latitude. Each individual footprints corresponds to a
brightness temperature at a given emission angle (Figure 3). From these
individual observations we obtain the desired quantities $T_{n}$ and $p$ at
the six MWR frequencies and their uncertainties by using a weighted linear
least squares (Equation 5) to fit the data to a parametric relationship
between emission angle and brightness temperature (Equation 3). Using the
logarithm form, we split the parametric equation into a summation of two
parts. The cosine of the emission angle functions as a weight (${W}$) to
downweigh observations (${Y}$) at large emission angles where we suspect
contamination by synchrotron radiation to be largest. Since synchrotron
contamination becomes less problematic at higher frequencies, we choose to
apply the weights only to Channel 1 - Channel 3.

Figure 1: PJ3 flyby during December 2016. The map shows X-band (8-10 GHz) VLA
observations taken concurrently with the flyby to show the longitudinal
structure of the planet (Moeckel2019, These VLA data will be fully published
in a subsequent paper). The red footprints correspond to MWR’s Channel 5
observation, where we have plotted the highest resolution (i.e. smallest
emission angle), while the blue footprints correspond to the VLA resolution in
its most extended A-configuration. Inset a) is zoom in for the indicated
region and shows all the 60 individual observations colorized by their beam-
convolved emission angle that are centered around the equator. These 60
observations are the basis for fitting the nadir brightness temperature and
the limb-darkening at 0∘. Inset b) is zoom in to illustrates how the emission
angle varies across the beam, where the red central dot corresponds to the
boresight center and the black dot indicate the position of the spacecraft.
The dashed ellipses correspond to the beams’ extent projected onto the planet
with the contours showing how the emission angle varies across the beam. The
final beam-convolved emission angle is 59.2 deg when integrating over the 2.5
$\theta_{HPBW}$ contour that corresponds to 99% of the emission.

$\frac{\delta}{\delta X}\sum\limits_{i=1}^{n}\vec{W}({Y}-{A}{X})=0$ (5)

$\displaystyle W=\begin{bmatrix}\mu_{1}\\\ \vdots\\\ \mu_{i}\\\ \vdots\\\
\mu_{n}\end{bmatrix}$ $\displaystyle Y=\begin{bmatrix}logT_{b}(\mu_{1})\\\
\vdots\\\ logT_{b}(\mu_{i})\\\ \vdots\\\ logT_{b}(\mu_{n})\end{bmatrix}$
$\displaystyle A=\begin{bmatrix}1&log\mu_{1}\\\ \vdots&\vdots\\\
1&log\mu_{i}\\\ \vdots&\vdots\\\ 1&log\mu_{n}\end{bmatrix}$ $\displaystyle
X=\begin{bmatrix}logT_{n}&p\end{bmatrix}$

Based on best fit parameters we compute the error $\epsilon$ by subtracting
the data from the parametric fit, and compute the reduced chi-squared value
$\chi^{2}$ for $n_{d}$ data points and $n_{p}$ parameters:

$\displaystyle\epsilon$ $\displaystyle={Y}-{A}{X}$ $\displaystyle\chi^{2}$
$\displaystyle=\frac{\epsilon^{T}W\epsilon}{$n_{d}$-$n_{p}$}$
$\displaystyle(A^{T}WA)^{-1}\chi^{2}$
$\displaystyle=\begin{bmatrix}\left(\sigma^{fit}_{\bm{log(T)}}\right)^{2}&\bm{\rho}\\\
\bm{\rho}&\left(\sigma^{fit}_{p}\right)^{2}\\\ \end{bmatrix}$

Figure 2: Verification of the fit uncertainty derived from the weighted least
squares approach by comparing the results to an Markov-Chain Monte-Carlo
(MCMC) method using the same underlying model to fit the data Equation 3. The
solid black lines corresponds to the 16%, 50% and 84% quantiles, while the
dotted red line indicate the errors retrieved from the weighted least square
approach. The best-fit quantities for both methods agree well within their
uncertainties. The largest difference is found in the uncertainty on the fit
limb-darkening parameter, where the weighted least square approach returns
smaller errorbars than the MCMC results.

The 1-$\sigma$ uncertainty corresponds to the diagonal element of the
covariance-variance matrix, while the off-diagonal elements are the
correlation between the factors. From there, we obtain the uncertainties for
the brightness temperature $\sigma^{fit}_{T}$ and for the limb-darkening
parameter $\sigma^{fit}_{p}$. In Figure 2 we compare the fit uncertainties
derived from the weighted least square approach to a more computational
expensive MCMC approach. The best-fit quantities derived from each other agree
within the derived uncertainties, with the uncertainties derived from the
weighted least squares are smaller when compared to the spread seen in the
MCMC results. We therefore will carefully assess the impact on the
uncertainties on the derived ammonia distributions. The MCMC results also
indicate the correlation between the two quantities
The final uncertainty for the brightness temperature is a combination of the
fit uncertainty $\sigma_{T}^{fit}$ and the absolute calibration uncertainty
$\sigma_{T}^{cal}$ (Janssen2017), which we add in quadrature. We adopt an
absolute calibration uncertainty $\sigma^{cal}_{T}$ of 10% for Channel 1 due
to the unknown opacitiy of absorbers (e.g., ammonia, water, and alkalis) at
such high pressures (Hanley2009; Bellotti2016; Bellotti2017; Bellotti2018),
and 1.5% for all other channels (Janssen2017). The final uncertainty for the
limb-darkening coefficient is a relative measurement and therefore corresponds
only to the fit uncertainty.

$\displaystyle\sigma_{T}$
$\displaystyle=\sqrt{(\mathbf{n}\;\sigma_{T}^{fit})^{2}+(\sigma_{T}^{cal})^{2}}$
(6) $\displaystyle\sigma_{p}$ $\displaystyle=\;\;\mathbf{n}\;\sigma_{p}^{fit}$
(7)

Spatial variation of the observables within the beam and remaining synchrotron
contribution both deteriorate the reliability of the fit, making the retrieval
prone to errors in case of an overreliance on these measurements. To fully
understand the impact of the uncertainties, we allow $\mathbf{n}$ to vary
between 1 and 10 (inset in Figure 3 shows the various uncertainty levels), and
test the robustness of our results to the uncertainty fidelity $\mathbf{n}$.
Thus, as function of latitude we have reduced the MWR measurements to two
parameters: nadir brightness temperature $T_{n}$ and limb-darkening
coefficient $p$ (with a 1 deg step-size in latitude per PJ) and their
corresponding uncertainties.

Figure 3:  Brightness temperature at 12∘ deg north latitude for the 6
different channels, plotted against the beam-convolved emission angle. The
closed circles correspond to the observations considered in the fit, and the
open circles to the observations that were removed. Open circles inwards of
the cut-off are observations with high susceptibility to synchrotron radiation
(spatial filter – see Section 2.3), while open circle outwards are
measurements where part of the beam is off the planet. The cut-off is
determined such that all beam elements interior 2.5 HPBW are on the planet.
The impact of the synchrotron radiation leaking in through the side lobes can
be seen in Channel 1 and 2, where the brightness temperature increases beyond
60 deg, while the drop in brightness temperature at higher frequencies is due
the cold sky leaking in. The solid line with shading corresponds to the fitted
model according to Equation 3 and three uncertainty fidelities (n=1,3,10). The
intersect at an emission angle of zero is the nadir brightness temperature,
and the limb-darkening parameter defines the shape of the curve.

In addition to the individual PJs, our final goal is to produce a global
(i.e., longitude, or orbit-averaged) latitudinal distribution of trace gases.
We therefore combine the individual PJs and compute the mean for each channel.
We caution the reader that the global mean is based on nine narrow slices of
the atmosphere separated both in time and longitude, and is thus not a perfect
representation of a global atmosphere. For the uncertainties, we elect not to
use the standard deviation between the individual PJs since they (i) reflect
spatial variations instead of measurement uncertainties and (ii) do not
reflect the absolute uncertainty for the temperature measurements. Instead we
compute the average fit uncertainty $\overline{\sigma}^{fit}$ at a given
latitude among the PJs and add them to the absolute uncertainty:

$\displaystyle\sigma_{T}^{mean}$
$\displaystyle=\sqrt{\left(\mathbf{n}\;\overline{\sigma}_{T}^{fit}\right)^{2}+\left(\sigma_{T}^{cal}\right)^{2}}$
(8) $\displaystyle\sigma_{p}^{mean}$
$\displaystyle=\;\;\mathbf{n}\;\overline{\sigma}_{p}^{fit}$ (9)

We use $\mathbf{n}$ = 3 (i.e., 3-$\sigma$) for the fit uncertainty when
plotting the trace gas distribution, but we also show the retrieved
atmospheric quantities for different fidelities (n=1,5,10) in Section 5.

###  2.3 Data selection

Despite limiting observations to viewing geometries with 99% of emissions
coming from the planet, the low frequency observations show signs of
synchrotron contamination. The tell tale sign of synchrotron contamination is
a breakdown of the typical relationship between emission angle and brightness
temperature. Observations in which synchrotron radiation off the planet is
stronger than the atmospheric signal show a non-physical limb-brightening
effect, unlike the limb-darkening that we are see for atmospheric
observations. In the absence of accurate synchrotron radiation models, and
given the complexity of the electron distribution and the magnetic field close
to the planet (e.g., dePater1981; dePater1997; SantosCosta2017; Moore2019), we
use two different filters to remove observations that are suspect to
contamination.

#### Spatial Filter

The main synchrotron radiation belt of Jupiter is constrained to low latitudes
around the magnetic equator (e.g., roberts1976pitch; sault1997first;
dePater1998). Despite the low sensitivity of the antenna beam in the side-
lobes, the combination of high intensity radiation and large emission area
causes problems in observations where synchrotron emission enters through the
side lobes. We therefore restrict measurements where the sensitivity to the
synchrotron radiation in the side lobes is minimal. In practice, when the
spacecraft is at higher latitudes (10-60∘), we remove all observations that
are looking towards the equator (where the equatorial synchrotron radiation is
leaking in) and instead use only observations that look polewards. See Figure
3 and Appendix B for examples of this filter.

#### Limb-Darkening Filter

The fitted limb-darkening coefficient and its deviation from the mean is
another constraint that we use to assess the synchrotron contribution to the
observations. Observations that are bright at the limb (i.e., a lower than
expected value for the limb-darkening) are a telltale sign of contribution by
synchrotron radiation. Latitudes where the limb-darkening coefficient deviates
more than 5 standard deviations from the mean for the given orbit are replaced
by the orbit mean at that latitude.

###  2.4 Orbit average

The geometry during each flyby is dictated by the type of science that is
prioritized. For MWR observations, we therefore restrict ourselves to a subset
of nine orbits from the first 12 to compute the orbit average (see Table 2).
We use the retrieved uncertainties in the parameters as weights to compute the
mean profiles for brightness temperature and limb-darkening after applying the
spatial and limb-darkening filters discussed above.
The results of the data reduction for these Juno data are shown in Figures 4
(nadir brightness temperature) and 5 (limb-darkening coefficient), where the
solid lines correspond to the mean and the thin dashed lines correspond to the
various PJ passes. For verification purposes, we also compare our results to
published orbit-averaged data of Oyafuso2020, indicated by the dotted line.
The variations with latitude of radio-cold and radio-warm regions indicate the
presence of enriched and depleted trace gases, which modify the opacity for
the escaping blackbody radiation. Some of the latitudinal variations seen in
the mean profiles can be traced down across all channels, such as the radio
cold region around the equator, indicating the presence of atmospheric
variations spanning several orders of magnitude of pressure. At higher
latitudes, the profiles show that the variations decay rapidly with depth,
highlighting more shallow atmospheric features. The majority of atmospheric
variability, as seen by the dispersion between individual PJ profiles, is seen
primarily at the higher frequencies, which correspond to the upper layers of
the atmosphere where temporal and spatial variability dominate (dePater2016).
The structure of the limb-darkening is more chaotic, but to first degree it
mirrors a structure similar to the brightness temperature, with the biggest
diversions seen in the equatorial zone and the neighboring north equatorial
belt. We can see smaller features at higher latitudes corresponding to regions
where the white zones and reddish belts alternate. The limb-darkening
coefficient is very sensitive to variations in spatial structure within the
binned regions, so that it becomes difficult to interpret the signal at higher
frequencies, and we caution reading too much into the fine scale structure at
the highest frequencies.

Figure 4: Nadir brightness temperature for the 6 channels between -60 and 60
deg, where the thick lines indicate the orbit averaged values, and the dashed
thin lines correspond to the individual PJs. The dotted line represents the
orbit averaged value based upon a deconvolution algorithm (Oyafuso2020) and
shows a good agreement between the two methods, validating both our
approaches.  Figure 5: Same as Figure 4 but for the limb-darkening parameter.

##  3 Methods

The presence of radio-absorbing trace gases in a planet’s atmosphere modulates
the thermal emission from the atmosphere as it propagates upwards. By
prescribing a trace gas distribution with a temperature profile for the
atmosphere of Jupiter, we can simulate its brightness temperature using a
radiative transfer (RT) model, such as that described by e.g., Berge1976;
dePater1985. Such simulated data (see Section 4 for the RT code) can then be
compared to Juno MWR and VLA observations: by adjusting the abundances and
altitude profiles of radio-absorbing trace gases (and/or the temperature
profile), one can find a best fit to the data. In Section 5 we focus on the
$NH_{3}$ abundance, which is the dominant absorber at the wavelengths Juno MWR
and the VLA are operating (e.g., dePater1986). At the longest wavelengths the
presence of water (vapor and cloud particles) starts to have an effect on the
observations as well (dePater2005). This signal however is difficult to
interpret and we dedicate a separate Section 6 to the water abundance.

But what is the best way to prescribe the distribution of trace gases? In the
past, researchers determined “arbitrary” jumps in the ammonia profile (e.g.,
dePater1986) to find a best fit to their data. dePater2001 pointed out that
profiles which gave a best fit to the available data did not follow the
profile expected from thermo-chemical equilibrium calculations
(Weidenschilling1973; Atreya1985). Li2017 followed a similar route, opting to
split the atmosphere in several layers and then use a Markov-Chain Monte-Carlo
(MCMC) algorithm to find the stochastic best fit to Juno PJ1 observations.
Considering the few Juno observations that we have (i.e., data at 6
frequencies), the problem is under-constrained, that is, we have more
variables to fit than observations constraining them. Additionally, the
received radiation is an integrated signal over many layers of the atmosphere,
and as such the inversion is not unique. Further information can be obtained
The problem can be further constrained through regularization of the
atmosphere, limiting how much the atmosphere can change between latitudes
pressure levels (Li2017), or through the use of thermo-chemical equilibrium
models as mentioned above (dePater2016; dePater2019). We know from
observations of the troposphere on Earth that the atmosphere is more complex
than a simple thermo-chemical equilibrium model prescribes. Nevertheless, it
adds information a-priori information, and we opt to modify the thermo-
chemical equilibrium models by introducing factors that quantify the departure
from equilibrium in these models.

###  3.1 Model 1: Perturbed thermo-chemical equilibrium model

In the absence of other observations, the simplest description of an
atmosphere is a well-mixed atmosphere, where the trace gas (†††We express
trace gases in both parts per million (ppm) and in proto-solar abundances:
$\frac{N}{H^{2}}$ = 1.48e-4, $\frac{S}{H^{2}}$ = 2.89e-5, $\frac{O}{H^{2}}$ =
1.07e-3 (Asplund2009)) abundance ($x_{i}$) is constant with pressure until the
saturation pressure, $P_{sat}$, at temperature $T_{c}$ for the specific trace
gas is reached, at which point the gas starts to condense out to form clouds.
We modify the relationship by introducing the relative humidity $r_{h}$, where
a value greater than unity indicates super-saturation, delaying the
condensation of cloud particles. On Earth, super-saturation occurs due to a
variety of reasons, most commonly the absence of cloud condensation nuclei, or
the presence of strong updrafts. Sub-saturation at these levels, corresponding
to values less than unity, are most likely a reflection of averaging areas
inside the radio beam that are saturated and areas that are not saturated,
showing small-scale regions of down and upwelling.

$v=r_{h}\frac{p_{sat}(T)}{P}$ (10)

Below the ammonia ice clouds, the absence (or lack of detection) of hydrogen
sulfide (H2S) gas on Jupiter has been explained by the formation of ammonia
hydrosulfide clouds (NH4SH), where essentially all the H2S is trapped in the
NH4SH as it combines with the more abundant NH3. The formation of the cloud is
controlled by the partial pressure of the two gases, and the reaction rate has
been parameterized by Weidenschilling1973 and Atreya1985.

$NH_{3}+H_{2}S\rightleftharpoons NH_{4}SH$ (11)

The only measurement of sulfur on the planet so far comes from the Galileo
probe measurements (Niemann1996; Wong2004) in a so-called 5-$\mu$m hot spot
(i.e., a region devoid of clouds), which showed a solar enrichment of H2S
similar to that of NH3. Since our results for the ammonia abundance at high
pressures ($\gtrsim$20 bar) are lower than the estimates from the Galileo
experiment (see Figure 8), we opt to not use the Galileo probe measurement.
Instead, we retrieve the nitrogen enrichment factor compared to solar at
pressures below any cloud formation (Section 5), and apply the same enrichment
factor to compute the abundance of H2S. In order to Instead we retrieve the
H2S abundance indirectly by observing its effect on the NH4SH cloud formation.
Additionally, we account for uncertainties in the NH4SH cloud formation, we
introduce a factor that allows cloud formation to be delayed,
$\eta_{NH_{4}SH}$, the efficiency of the reaction. At a value of one, the
reaction occurs at the predicted pressure, and for values less than one, the
clouds form higher up in the atmosphere than predicted by thermo-chemical
equilibrium models:

$P_{NH_{3}}P_{H2S}\geq\eta_{NH_{4}SH}P_{NH_{4}SH}$ (12)

We did not find a constant ammonia abundance below the NH4SH layer that could
explain both the VLA and Juno observations; therefore, we introduced a
parameterization for the region between an assumed deep, globally uniform
atmosphere and the NH4SH cloud layers. Processes such as dynamics
(Showman2005) or microphysics (Guillot2020) can locally enrich or deplete the
atmosphere, leading to gradients in abundance with pressure. We started with
the simplest model by fitting a single mixing parameter between a top and a
bottom pressure, and a deep ammonia abundance. The mixing gradient is defined
such that a positive gradient constitutes a decrease in ammonia with
decreasing pressure. We also tested a second model where we prescribed a deep
ammonia abundance and fit two different scale heights to distinguish between
different mixing regimes in the atmosphere.

$\displaystyle NH_{3}(P)$ $\displaystyle=NH_{3}^{deep}e^{H_{mix}(P_{b}-P)}$
(13) $\displaystyle H_{mix}$ $\displaystyle=\frac{\delta NH_{3}}{\delta P}$
(14)

We postulate that at some depth, the atmosphere is well mixed across all
latitudes. The most well-suited observations for constraining a global deep
abundance come from MWR Channel 1, which probes deepest into the atmosphere
(see Figure 7). Unfortunately, the inversion of the longest wavelength
observations is also the most uncertain, due to scarcity of high pressure and
temperature experiments simulating the conditions in the gas giants
(Hanley2009; Bellotti2016; Bellotti2017). Different published ammonia
absorption coefficients lead to differences as large as 100K for models at
Channel 1 frequencies. We therefore opt to reduce the influence of the longest
wavelength, by applying large absolute uncertainties to the Channel 1
measurements, so that we instead rely on the shorter wavelength channels.
Lastly, for the atmospheric inversion based on the individual PJs, we opt to
prescribe the deep ammonia abundance and derive the pressure below which the
atmosphere is well mixed. This is in line with our assumption that there
exists a deep, well-mixed reservoir of trace gases on the planet.

While H2S gas on Jupiter does not contribute noticeably to the absorption at
radio wavelengths, water has a non-negligible effect on the observations at
the longest wavelengths (see, e.g., dePater1985; dePater2005; Janssen2005).
For the main analysis, we enriched the water abundance by the same solar
enrichment factor that we retrieved for the ammonia gas, and let it condense
out at the water condensation point. This is the simplest assumption and is in
line with the results for the Equatorial Zone (EZ) (Li2020).We test this
assumption in Section 6 where we try to retrieve the water abundance and
ammonia abundance simultaneously. For this part of the analysis, we included
only an indirect effect of the water cloud on the ammonia: about 3% of gaseous
ammonia can be absorbed into the liquid water rain drops
(Weidenschilling1973). This effect causes a small local decrease in ammonia at
the water cloud level, $\sim$6 bar.

###  3.2 Model 2: Stochastic model for the atmosphere

Additionally to the perturbed thermo-equilibrium model we developed a model
similar to the one used by Li2017, that fits gradients between pressure
levels. Starting from 100 bar we introduce pressure nodes between which we
assume an exponential gradient capturing that mimic processes locally
enriching or depleting the atmosphere in trace gases. We set the pressure
nodes to 40,20,15,10,6,4,2 and 1 bar above which the atmosphere follows the
saturation-vapor pressure curve as defined in Equation 10.

##  4 Radiative Transfer Model

We use our open-source radiative transfer model radioBEAR (Radio BErkeley
Atmospheric Radiative-transfer‡‡‡https://github.com/david-deboer/radiobear) to
simulate the brightness temperature for a given distribution of trace gases
and a prescribed temperature structure in the atmosphere by splitting the
atmosphere into multiple layers and calculating the contribution from each
individual layer. This model has been described by dePater2005; dePater2014;
dePater2019. The opacities for ammonia and water were derived from lab
experiments that simulate the high-pressure environment for a given hydrogen-
helium atmosphere by Devaraj2014; Hanley2009; Bellotti2016; Bellotti2017. At
the longest wavelength, the Bellotti2016; Bellotti2017 modeled brightness
temperatures were 100K warmer than the Hanley2009 calculations and are not
reconcilable with the Juno observations. Hence, we favored the earlier lab
measurements.

###  4.1 Temperature Structure

The temperature structure in the atmosphere is a crucial element for the
retrieval, since both trace gas and temperature variations can alter the
observed brightness temperature, and we can only observe the combined effect.
While observers have used mid-infrared wavelengths to derive the temperature
at pressures $P\lesssim 0.7$ bar (Fletcher2009), the temperature structure
below the radiative-convective boundary cannot be obtained with these
techniques.

The temperature at the depth which radio waves originate from are controlled
by atmospheric dynamics, microphysics, and wave patterns in the atmosphere,
and to date, there are very few temperature measurements below the radiative-
convective boundary ($\sim$0.5 bar). The first measurements come from the
Voyager flyby, which used radio-occultation measurements to probe the
temperature structure down to the 1 bar level on ingress and egress to measure
the temperature of the North Equatorial Belt (NEB, 165K $\pm$ 5) and the
Equatorial Zone (EZ, 170.4K $\pm$ 5), respectively (lindal1981atmosphere;
Guillot2004). The latent heat release of condensation in the EZ decreases the
temperature gradient in the atmosphere, leading to higher temperatures at the
same pressure levels compared to an atmosphere that follows a dry adiabat.
Hence, the observed temperatures in the EZ and NEB are consistent with a local
enrichment of trace gases in the EZ compared to the NEB. The only true in-situ
measurements come from the Galileo probe experiment, which entered a highly
trace-gas depleted region of the atmosphere, somewhat similar to the
conditions found in the NEB, and measured a temperature structure resembling a
dry adiabat down to 20 bars.

We have built two temperature profiles—a wet and a dry adiabat—that are
anchored to the Voyager and Galileo measurements (Figure 6). Starting from
deep within the planet, the absence of condensible species implies that the
temperature follows a dry adiabat, until the water condensation level is
reached, at which point the two temperature profiles start to diverge. While
the dry adiabat reaches 166.1 K at the 1 bar level, the latent heat from water
decreases the temperature gradient for the wet adiabat, so that in the EZ a
temperature of 170.4 K is reached at the 1 bar level. We note that the latent
heat release from the more abundant water content on Jupiter is about an order
of magnitude larger than the latent heat release of ammonia. Higher up in the
atmosphere (P $\lesssim$ 0.7 bar), we use profiles obtained from the mid-
infrared (Fletcher2009) and merge the two temperature profiles to cover the
full depth which radio waves probe. In the radiative zone the profiles
reverse, that is, the NEB becomes warmer than the EZ, as the heat can more
efficiently escape from the NEB and heat up the stratosphere.

We tested our results for a dry and a wet adiabat, however, found that the dry
adiabat produces consistently a better results across all latitudes, including
the equatorial zone. All following results are based on the dry atmospheric
profile.

Figure 6: Compilation of the relevant temperature profiles which we assume for
our retrieval. The temperature profiles are anchored at the 1 bar level to the
radio occultation measurements of Voyager (lindal1981atmosphere) and to the
Galileo probe measurements (Seiff1996). Above the 0.7 bar level we use
profiles obtained via infrared measurements (Fletcher2009). Below the 1.0 bar
level we use the expected wet and dry adiabatic lapse rate until the two
profiles merge at the height of the water cloud. The little cutout shows the
anchor points at the 1 bar level, with the corresponding uncertainties.

###  4.2 Weighting functions

Radio emissions originate over a range of atmospheric pressures, and the final
signal is the integrated emission over these layers. The weighting functions
illustrate this effect numerically, by showing the depth of and range from
where we receive the emission. We plot the weighting functions in Figure 7 for
the EZ and the NEB, exemplifying the end members in terms of ammonia gas
distribution covering the most enriched and most depleted region of the
planet. At the longest wavelengths, the radio waves originate over a wide
range of pressures going down to a few hundred bars, while at shorter
wavelengths the weighting functions become more peaked, indicating that they
probe a smaller pressure range. The center of the absorption feature is at 22
GHz, so that observations at this frequency probe highest in the atmosphere,
while observations both at smaller and longer wavelengths probe deeper into
the atmosphere. Lastly, the depletion of trace gases in the NEB allows for
probing deeper into the atmosphere, enabling retrieval of ammonia abundance at
higher pressures.

Figure 7: Weigthing functions visualize the pressure levels from which
radiation is received. We plot the functions for both end members on the
planet, a highly enriched part in the atmosphere (solid line corresponding to
the EZ), and a highly depleted part (dashed line corresponding to the NEB)
showing the depth and range from which we receive the signal. The abundance
profile for which we compute the weighting functions are taken from our fit to
the orbit-averaged MWR data (see Figure 9). The left hand side shows the
relevant functions for the Juno spacecraft, and the right hand side shows the
function for the VLA.

###  4.3 Optimization

For a given distribution of trace gases, we can simulate the brightness
temperature and the limb-darkening coefficient; these modeled parameters are
referred to as $T_{m}$ and $p_{m}$. By matching our simulated atmosphere to
the observations, we retrieve the best-fit distribution of ammonia gas. We
opted to retrieve the distribution per 1-deg bins, independent of neighboring
latitudes, for a robust retrieval. We first determined the optimum using the
scipy §§§https://www.scipy.org/ Nelder-Mead optimizer to find an atmospheric
composition that minimized the difference between the model ($T_{m}$, $p_{m}$)
and the observations ($T_{b}$, $p$), weighted by the uncertainties
($\sigma_{T_{b}}$, $\sigma_{p}$). The uncertainties are essential for
combining the various observational quantities and weighing them
appropriately.

$C=\sum_{1}^{6}\frac{(T_{m}-T_{b})^{2}}{\sigma_{T}^{2}}+\sum_{2}^{5}\frac{(p_{m}-p)^{2}}{\sigma_{p}^{2}}$
(15)

We only use the limb-darkening information from Channel 2 to Channel 5. We
discarded Channel 1 limb-darkening measurements as we cannot match the opacity
of Channel 1 with only ammonia. Other absorbers such as water (dePater2005,
(both vapor and droplets)) and alkalis (Bellotti2018) are known to add opacity
at these low frequencies. We also discarded Channel 6 because it is very
insensitive to our variations in the ammonia abundance and mostly driven by
spatial variations (see the VLA maps published in dePater2016; dePater2019).
This established a combination of parameters that minimize the difference, but
does not indicate the range of possible solutions within the constraint. We
therefore coupled our retrieval with a Markov-Chain Monte-Carlo (MCMC)
(Foreman-Mackey2013) for the set of parameters $x$ to find the uncertainties
in the fit parameters and the correlations between the parameters. The MCMC
maximizes the log-likelihood function $lnp$:

$lnp(T,p|x,\sigma)=-0.5\left(\sum_{1}^{6}ln(2\pi\sigma_{T})+\frac{(T_{m}-T_{b})^{2}}{\sigma_{T}^{2}}+\sum_{2}^{5}ln(2\pi\sigma_{p})\frac{(p_{m}-p)^{2}}{\sigma_{p}^{2}}\right)$
(16)

##  5 Results

In this section we present our best-fit ammonia distribution for the average
of the first 12 orbits when MWR was able to observe the atmosphere (see Table
2) and compare our results for the first PJ to the inversion by Li2017.

###  5.1 Well-mixed reservoir

The first quantity we constrain is the amount of ammonia in the region that is
globally well-mixed. Based on Figure 7, we can see that even at the long
wavelengths of Juno, the signal from depths over 50 bar is poorly constrained,
due to the wide range of pressures the signal originates from. There are four
observed quantities from the orbit-averaged quantities that constrain the deep
ammonia abundance: the brightness temperatures of Channel 1 to Channel 3, and
the limb-darkening coefficient of Channel 2. Since our model cannot find an
ammonia distribution that fits all four quantities, we instead search for a
solution that minimizes the overall error. The key parameter we vary is the
factor $n$ for the fit uncertainties ($T_{b}\pm n*\sigma^{fit}_{T}$,$p\pm
n*\sigma^{fit}_{p}$) and we test the sensitivity of our solution to the chosen
fidelity factor $n$. As a reminder, a value of 1 corresponds to the 1-sigma
results when fitting the nadir brightness temperature and limb-darkening
(Equation 8). We consider four different combinations:

  * 1-$\sigma$

This option uses the mean error in the fit uncertainty based on the nine limb-
darkening profiles (Note: this is not the standard deviation of the fitted
limb-darkening coefficient). This combination maximizes the fit to the limb-
darkening coefficient of Channel 2, at the cost of fitting the brightness
temperature in Channel 1 and Channel 2. We interpret this number as a lower
limit to the deep abundance.

  * 3-$\sigma$

This combination compromises the fit between all quantities, and we use this
for the final result.

  * 5-$\sigma$

Similar to 3-$\sigma$, where the limb-darkening coefficient still influences
the retrieval.

  * 10-$\sigma$

With these large uncertainties, the brightness temperature is the main
quantity driving the retrieval. While the limb-darkening uncertainty changes
by a factor of 10, the total uncertainty in the temperature changes by less
than 2% due to the large calibration errors. This retrieval we interpret as an
upper limit to the retrieved ammonia abundance.

For our final results, we first fit the ammonia abundance between -45 and 45
degrees independently at 1 deg steps, and then average the result over the
chosen latitude range. We adopt the value at 100 bar as the deep NH3
abundance. We compile previous results for the averaged deep ammonia abundance
as a function of pressure in Figure 8 and compare that to the four recovered
profiles. The 3-$\sigma$ profile converges at $366.5^{+32.8}_{-19.9}$ ppm,
where the uncertainties correspond to the 1-$\sigma$ and 10-$\sigma$ profiles.
The final value corresponds to an enhancement over solar of
$2.47^{+0.22}_{-0.13}$. We can see that our results fall in between the VLA
results (dePater2019), which is very insensitive to abundances at these
depths, and the results from the preliminary inversion of the temperature only
measurement of PJ 1 (Li2017).

Figure 8: Globally averaged [-45∘,45∘] ammonia profile as a function of depth.
The various shades of red correspond to our results for different factors of
$n$. The blue curve is from the preliminary inversion of Li2017, and the green
curve is based on the VLA data from dePater2019. The radio attenuation
(Folkner1998) and Galileo Probe Mass Spectrometer (GPMS)(Wong2004) are
superimposed.

###  5.2 Orbit mean

With the deep ammonia constrained to a constant global value, we can turn our
attention to obtaining a profile averaged over Juno’s primary mission by
fitting parameters to the average of the first year of the mission between
-60∘ and 60∘ ( Table 2). Poleward of 60 degrees the spacecraft’s larger
distance from the planet reduces the resolution and increases contamination
from the sidelobes, reducing the reliability of the data inversion. The nine
orbits are widely separated in longitude reflecting a global latitudinal
structure, but we caution that the observations were taken over the course of
a year, during which time the atmosphere underwent some large scale changes,
such as the storm outbreak between the North Temperate Belt (NTB) and the
North Tropical Zone (NTrZ) (Sanchez-Lavega2017; dePater2019_2), and the
outbreak in the South Equatorial Belt (SEB) (e.g., dePater2019_2). In Figure 9
we compare the average atmosphere based on the 3-$\sigma$ option to the
results obtained from the VLA, which were taken between Dec. 2013 and Dec.
2014 dePater2016; dePater2019. The VLA results are based on the brightness
temperature at 21 different frequencies between 3 and 25 GHz, i.e., at a
frequency resolution of 1 GHz; these data are sensitive to pressures
$\lesssim$ 8-10 bar (and down to $\sim$20 bar in the NEB). In contrast, the
Juno observations cover a much wider frequency range, from 0.6 to 22 GHz, at 6
distinct frequencies (Table 2), probing down to 100’s of bars (see Figure 7).
Although these data have a lower frequency resolution, there is further
information in the Juno data by adding a second quantity to fit in models: the
variation of the limb-darkening coefficient with latitude.
The Juno observations confirm the overall structure of the planet that the VLA
has seen at $P\lesssim$ 8 bar, while adding more depth and more details to the
observations. The overall structure is characterized by alternating columns of
enriched and depleted ammonia gas. The equatorial zone is highly enriched in
ammonia, corresponding to the cold radio-temperatures around the equator. The
Juno observations add fine scale structure, such as the increase of ammonia
with altitude in the EZ at P $\lesssim$ 50 bar, the only location on the
planet with such peculiar structure. We can also observe the expansion of the
equatorial column into the NEB at pressures between 2 and 20 bars. The ammonia
abundance in the EZ peaks just below the NH4SH clouds, constituting an
enhancement in the upper troposphere compared to the deep ammonia abundance.
This feature (’ammonia enhancement factor’) was instrumental to derive the
water abundance in the EZ based on Channel 2 - 6 observations by Li2020. The
largest depletion of ammonia on the planet is in the adjacent North Equatorial
Belt, where we fit a depletion of trace gases down to several tens of bars,
quite similar to that derived from the VLA data, though slightly shifted
northwards in latitude at larger pressures. At pressures of 40 bars, the
depletion tapers off to the deep ammonia abundance. At higher latitudes the
variation in brightness temperature with latitude become smaller resulting in
a more homogeneous atmosphere, with most changes seen in the upper few bars,
similar to the VLA results (e.g., dePater2016).

Figure 9: Ammonia abundance map based on the Juno’s orbit average brightness
temperature and limb-darkening with a 3-$\sigma$ fit uncertainty level. The
map shows familiar features such as the column of enriched ammonia in the EZ
and the depletion of ammonia in NEB, with smaller scale features at the mid-
latitudes. Plotted on the same color scale is the original VLA map based on
observations between 2013 and 2014 (dePater2016; dePater2019). The overall
structure is very similar between the two maps, with the Juno data allowing to
resolve more structure at depth due to its low frequency observations. The
dotted lines indicate the latitude at which the zonal wind speed peaks, and
with that denotes the difference between the zones and belts. The dashed lines
corresponds to the vertical cut as shown in Figure 13.

We assess the goodness of our fit by computing the brightness temperature and
limb-darkening coefficient for the given ammonia distribution in Figure 9.
Next to the best-fit results that is the basis for our interpretation (i.e.,
the 3-$\sigma$ option), we show forward model results for the other
combinations of uncertainty (n=1 and n=10) in Figures 10 and 11. Note that we
constrain the value for the deep ammonia corresponding to the values as shown
in Figure 8, so that the deep ammonia for n = 1 corresponds to 346.6 ppm (2.34
solar abundance), and the n=10 to 399.3 ppm (2.69 solar abundance). In Figure
10 we show that the factor $n$ for the fidelity of the fit uncertainty matters
little for the higher frequencies in fitting the observations, and our best-
fit distribution can explain the observed brightness distribution within the
uncertainties. The biggest impact of our confidence in the uncertainties
appears in Channel 2, where the various combinations trend to a more depleted
deep atmosphere (n=1, 1-$\sigma$) or a more enriched deep atmosphere (n=10,
10-$\sigma$). Comparing the brightness temperature fit in Channel 2 to the
limb-darkening fit in Channel 2, we can see that they anti-correlate. The
temperature fits best when using a large scaling factor $n$ (such as 10) for
the fit uncertainty, in which case the brightness temperature drives the
retrieval. A scaling factor of 1, on the other hand, means we fit the limb-
darkening at the cost of brightness temperature. We also show the brightness
temperature (green line) for the distribution determined by the VLA
(dePater2019). We can see that the VLA observations, while not being in
perfect agreement reflect the overall structure of the atmosphere pretty well.
The VLA model (based on the lower end of the Galileo Probe data (Wong2004)) is
generally too cold for Channel 1 to Channel 3, due to its lack of sensitivity
to radiation from the deeper layers in the planet (see Figure 7).
While latitudinal structure in the limb-darkening coefficient is more chaotic,
we still manage to fit the overall structure of the planet to first degree.
The large features around the equator and the neighboring belt are mostly
captured by the fit, only the smaller diversion at higher latitudes tend to
not be captured well by our fit. Similar to the brightness temperature, the
various linestyles for the red shaded lines correspond to the different
scaling factors for the fit uncertainties. This figure shows the additional
value of the limb-darkening coefficient in constraining the atmosphere. While
the brightness temperature for the VLA distribution was comparable to the Juno
observations (Figure 10), the limb-darkening coefficient for Channels 2 and 3
does not match, and indicates that the atmosphere is well-mixed at altitudes
much deeper than the 10 bar level that was inferred from the VLA data. Higher
up in the atmosphere the VLA observations agree well with the Juno
observations; near the ammonia ice cloud, the VLA shows larger latitudinal
variations than the Juno data, which likely results from the higher frequency
resolution. While we do not fit Channel 1, it is still interesting to see that
all models produce similar results, with much larger coefficients than
observed.

Figure 10: Forward modeling results based on the best fit distribution of
ammonia trace gases compared to the Juno observations, plotted against
latitude. The Juno nadir brightness temperatures are the solid lines and their
uncertainties represented by the shading. The various other lines correspond
to forward models, based on this paper(red colors) and the green line is the
forward model based on the VLA data (dePater2019). The various line-styles for
the red lines indicated the uncertainty limits as described in Section 5.1.
Overall, our forward model fits the brightness temperature well within the
uncertainties. Figure 11: Same as Figure 10 but for the limb-darkening
coefficient. We can see how the model has difficulties fitting both the
brightness temperature in Channel 2 and the limb-darkening. The results that
best fit the brightness temperature (10-$\sigma$) have the largest misfit for
the Channel 2 limb-darkening, and the 1-$\sigma$ results fit the limb-
darkening best but not the brightness temperature.

###  5.3 Perijove 1

To date only results from PJ1 are published (Li2017, between 40∘S and 40∘N),
and can be compared with our results. Figure 12 shows our map
(3-$\sigma^{fit}$) for PJ 1 together with the results from Li2017; the latter
were based solely on the brightness temperature Li2017. We can reproduce the
enhancement of ammonia at the equator and the depletion in the North
Equatorial Belt. The Li2017 map shows that latitudinal variations are
constrained to the top few bars of the atmosphere, above an ammonia inversion,
i.e., minimum in the ammonia abundance at all latitudes, with a well-mixed
deeper layer. In contrast, our model converges to a solution that
monotonically increases the ammonia abundance with increasing pressure. The
depth that we can trace this variation to, however, is deeper than the Li2017
results. This can also be seen in Figure 13 where we show three slices through
the representative regions of the planet: the equatorial zone (3∘), the north
equatorial belt (12∘), and the mid-latitudes (30∘). We also added a slice from
the VLA distribution to highlight how and where Juno observations have added
information.

Figure 12: Comparison of ammonia abundance map for PJ1, based on the
preliminary retrieval (Li2017) and the best-fit model based on this paper. The
overall structure looks very similar, with the biggest difference that we do
not find evidence for the ammonia inversion centered around 6 bars in the
preliminary retrieval.  Figure 13: Vertical cut through the ammonia abundance
maps corresponding to the locations indicated in Figures 9 and 12. The three
colors indicate the three different regions, where the line styles
differentiate between different sources. The thick solid lines are our models
with the thin lines corresponding to 100 random draws from the MCMC
posteriori. The dotted line is the atmosphere based on the same data (Li2017),
while the dashed lines are based on the zonally averaged VLA data
(dePater2019).

In Figures 14 and 15 we forward model the ammonia distribution in PJ1 using
our radiative transfer code and compare it to the observations for PJ1 (our
model: red dashed line; Li2017: blue dashed line). The four highest
frequencies corresponding to the upper layers of the atmosphere are well fit
within the uncertainties, and our results in all channels match the data
better than the results in Li2017. This highlights how degenerate inversions
based on radio observations can be, where quite different ammonia profiles can
approximate the observations well. We assessed the uncertainties in our
solution by performing a MCMC retrieval on the two most extreme end members in
the ammonia distribution: the equatorial zone and the north equatorial belt.

Figure 14: Same as Figure 10 but for the PJ1 observations. The individual
brightness temperature observations are given by the solid lines, their
uncertainty is given by the shaded region. The dotted line corresponds to the
brightness temperature as published by Li2017. The dashed lines are the
brightness temperatures for the abundances as plotted in fig. 12 and a dry
adiabat. The red line is our best fit, while the blue line corresponds to the
preliminary inversion (Li2017). We note that while the structure of the
atmosphere is very different, both models are close to the uncertainty, with
our model performing better at all frequencies. Figure 15: Same as Figure 14
but for the limb-darkening coefficient.

###  5.4 Uncertainties for the retrieved parameters of PJ1

In order to assess the uncertainties in our results for PJ1, we performed a
MCMC retrieval on the data for the EZ and NEB. MCMC results add two more
layers of information: uncertainties of the and potential correlations between
parameters. The diagonal elements of the corner plots in Figures 16 and 17
show the distribution of solutions for the given parameters with the error
bars corresponding to the mean, and the 0.16 and 0.84 quantiles. The off-
diagonal elements indicate the correlation between parameters. Starting at the
top of the atmosphere, we can see that the relative humidity in the ammonia
ice cloud is well constrained. We find a high relative humidity in the EZ, and
low humidity in the NEB, which is in line with the thick white cloud deck seen
at visible wavelengths around the equator and the lack of the white ammonia
clouds in the NEB. Such differences in relative humidity were first suggested
by Klein1978, and later verified from VLA observations (e.g., dePater1986;
dePater2016; dePater2019).

Sitting below the ammonia clouds is the postulated $NH_{4}SH$ cloud, that
forms around 2-3 bar. From our retrievals we find that while we retrieve high
values for the cloud forming efficiency (close to 1) based on our optimizer
results, the value itself is poorly constrained due to the low frequency
coverage of MWR that probe the atmospheric layers where the cloud forms. We
need the higher frequency coverage of VLA observations simultaneous with the
MWR observations to study this cloud formation process.

We define the region where mixing occurs by three parameters: the gradient
($H_{mix}$) between a bottom pressure ($P_{b}$) and a top pressure ($P_{t}$).
The standing column of ammonia around the EZ leads to poor constraints on the
mixing region, however all solutions in the MCMC run support an increase of
ammonia with altitude (see more in Section 6.4). In contrast, the mixing range
in the NEB occurs deeper in the planet, and both $P_{t}$ and $H_{mix}$ are
better constrained.

Figure 16: Uncertainties for and correlations between the retrieved parameters
in the Equatorial zone for PJ1 observations.  Figure 17: Same as Figure 16 but
for the North Equatorial Belt.

###  5.5 VLA 2013-2014 results

With the updated deep ammonia abundance based on the Juno data, we will now
re-analyze the latitudinal profiles of the VLA data taken $\sim$2 years prior
to the Juno orbit insertion dePater2016; dePater2019. The results are shown in
Figure 18, where the top panel is the ammonia map retrieved from the PJ1-12
averaged Juno data (copied from Figure 9, which used both the temperature and
limb-darkening information) and the bottom panel is based on the VLA
brightness temperature data alone. The only piece of information we added to
the VLA retrieval is the bottom mixing pressure at which the value of the deep
ammonia abundance $366.5^{+32.8}_{-19.9}$ ppm ($2.47^{+0.22}_{-0.13}\times$
solar abundance) is reached. The original model was already in pretty good
agreement at altitudes above the 10 bar level, showing regions of alternating
enriched and depleted columns of ammonia. However, applying this new
constraint, we recover a remarkable similar atmosphere, with many of the main
features retained, such as the ammonia enhancement at the equator, the
structure in the NEB, and the general low ammonia abundance down to a few 10s
of bars in the NEB and at mid latitudes. Not only do they agree well, but the
better frequency resolution allows to better constrain the upper troposphere
of Jupiter, the so-called weather layer. The biggest difference can be found
in the northern hemisphere, where abundances below the $NH_{4}SH$ cloud trend
to higher values (also shown in Figure 19). This is particularly interesting
since a storm outbreak in the northern hemisphere (Sanchez-Lavega2017)
happened around the time of perijove 2 (October 2016), and resulted in a
relative warming at northern latitudes in ALMA observations taken in early
2017 (see Figure 2 in dePater2019_2). The difference between the VLA and Juno
data suggest the storm to have had a long-lasting effect.

Figure 18: Ammonia abundance map based on the Juno data (top row) and based
purely on the VLA data (bottom row) using only the fact that atmosphere
becomes well at 40 bar. Many of the main features of the distribution can be
recovered with the VLA.  Figure 19: Slices through the atmosphere based on the
ammonia abundance maps shown in Figure 18. We can see a good agreement between
the profiles derived from the VLA data alone (dot-dashed lines) and the
profiles derived from MWR (solid line). The biggest difference is that the VLA
found an enhancement at mid-latitudes that MWR did not see, most likely a
consequence of storm outbreak during the early Juno mission (Sanchez-
Lavega2017).

##  6 Discussion

In the following we use our retrieved parameters for the ammonia distributions
as a proxy for atmospheric dynamics. As we have shown in the MCMC analysis
above, some of these parameters are poorly constrained. Our approach to bin
measurements into 1 deg and retrieve the parameters per bin results in a large
spread between adjacent latitudes, while the spatial resolution is much
coarser than 1 deg (Figure 1). To improve the signal we average the parameters
over 5 degrees of a latitude-rolling-average to show trends across the zones
and belts. Figure 20 shows the retrieved best-fit parameters for the set of
uncertainty combinations.

Figure 20: Retrieved atmospheric parameters as function of latitude (with the
individual zones and belts denoted by the vertical dashed lines) and for
various combinations of uncertainty (different colors and linestyles). Left:
a) Relative humidity of ammonia gas, b) Cloud formation efficiency, c) The
${\chi}^{2}$ value for the retrieved atmosphere. Right: d) Bottom pressure of
mixing region, e) Pressure range over which mixing occurs, f) Absolute change
in ammonia abundance in the mixing region.

###  6.1 Relative humidity

We obtain an increased relative humidity in the equatorial zone along with a
decreased relative humidity in the adjacent NEB and SEB compared to other
latitudes. Moving further towards the poles, the signal of the South Tropical
Zone and North Tropical Zone can be seen by the increase in relative humidity
until the relative humidity settles to a stable mean below 50% (see panel a)
in Figure 20). The Juno beam is large in comparison to cloud features, so that
we interpret a value of less than 100% saturation as a mixture of clouds and
cloud-free regions. This effect could explain the relative humidity results we
find at the mid latitudes.

###  6.2 NH4SH cloud formation

The MCMC runs revealed that the actual pressure where $NH_{4}SH$ cloud
formation ($\eta_{NH_{4}SH}$) starts is poorly constrained, however, our
3-$\sigma$ models return a consistent value close to unity. A value close to 1
indicates that the onset of cloud formation happens as predicted by the
thermo-chemical equilibrium models, while a value less than one indicates that
the reaction is happening at lower pressures than predicted. The nominal
pressure at which $NH_{4}SH$ is the more favorable state of the trace gases is
nominally around 2.3 bars for a value of 1, whereas a value of 1e-3 pushes
that cloud formation to a pressure of 1.4 bar. This can either mean that the
process is a lot less efficient than postulated (Weidenschilling1973), for
example if there are updrafts that ”supercool” the gas before NH4SH formation
takes place. Alternatively, this could also reflect a lower abundance of H2S,
which we only constrain indirectly. The disagreement between the various
linestyles in panel b) of Figure 20 however indicates that this result depends
much on the choice of uncertainties, and a higher spectral resolution than MWR
is needed to infer the cloud formation efficiency.

###  6.3 Mixing regime

Juno’s MWR covers a wide frequency range at the expense of spectral coverage,
resulting in very few constraints on the deeper atmosphere. We developed two
different models to describe the atmosphere below the ammonia ice cloud. The
simplest model prescribes a single gradient ($H_{mix}\equiv\frac{\delta
NH_{3}}{\delta P}$) starting from a bottom pressure ($P_{b}$), over a certain
pressure range ($\Delta P$). For clarity we report on the absolute change in
the ammonia abundance ($\Delta NH_{3}$) over the mixing regime, as the
gradient itself is a lot less constrained. The more complex model allows for a
secondary mixing regime between the two pressure levels; however, because we
found no significant improvement in our fits, we only present the simplest
model here. The mean bottom pressure at which we constrain the mixing to start
is around 20 bars, but with significant difference between the tropical region
and the mid latitudes (see panel d) in Figure 20). While the signal in the EZ
is poorly constrained by this model (see Section 6.4), we find that the mixing
in the adjacent belts (NEB, SEB) starts higher up in the atmosphere than the
planetary mean.

Interestingly, the mixing occurs at much greater depths than the water clouds
($\sim$7 bar), which was long thought to be one of the main drivers for the
tropospheric variations (Showman2005; dePater2019_2). These results,
nevertheless, confirm that the atmosphere is not well mixed below the cloud
levels (NH3, NH4SH, solution cloud), as inferred as early as the 1980s when
radio measurements showed the presence of a radio-hot NEB, which was
interpreted as a depletion in the ammonia gas abundance (dePater1986). VLA
measurements showed this depletion to extend down to the $\sim$20 bar level
(dePater2019), and we show that it might extend even deeper.
Despite its large size, and the inherent computational problems with modeling
planets this size, Young2019; Young2019p2 have managed to simulate the
dynamics of Jupiter, and found an ammonia depletion down to $\sim$8 bars, with
the exception of the EZ, implying that dynamics may play a large role in
shaping the atmosphere of Jupiter. Recently, microphysics has been invoked to
explain the depletion of trace gases by trapping ammonia in water ice and
raining it out (Guillot2020). It remains to be seen if this process can
deplete an atmosphere down to tens of bars, i.e., well below the melting point
of water ice.

While the mid-latitudes seem to have a very consistent signal for the mixing
region, with $\Delta NH_{3}\sim$ 100 ppm (panel f) in Figure 20), the EZ and
the NEB show the largest variations from that number. Most notable the EZ
shows sign of the gradient reversing, resulting in an increase of ammonia with
altitude, while the NEB shows the biggest gradient, due to the highly depleted
upper atmosphere.

###  6.4 Equatorial Zone structure

The equatorial zone in Figures 13 and 20 shows a peculiar structure with
ammonia increasing with altitude from $\sim$50 to $\sim$20 bar, in stark
contrast with the monotonically decreasing ammonia at all other latitudes. The
retrieval favors a solution where the ammonia abundance peaks in the upper
troposphere, just below the NH4SH cloud. While a similar structure was found
by Li2020, we want to show the fidelity of this solution. In our nominal
retrieval we fix the deep ammonia abundance and only allow the depth to vary
at which the atmosphere is uniformly mixed. We ran an MCMC retrieval at 3 deg
North, the region where the retrieval found the strongest signal consistent
with an inversion, and ran a full atmospheric retrieval. For the uncertainties
we chose the combination of 3-$\sigma$, which compromises between the limb-
darkening measurements and the brightness temperature measurements. In Figure
21 we show the most relevant MCMC posteriori distribution: the mixing gradient
for three different prescribed deep ammonia abundances.

(a) 366.5 ppm/2.47 solar abundance

(b) 400 ppm/2.69 solar abundance

(c) 428 ppm/2.89 solar abundance

Figure 21: Vertical mixing gradient (change in ammonia per pressure scale
height) for three prescribed deep ammonia abundances in the EZ (3∘N). Negative
values indicates an increase in ammonia with decreasing pressure, or in other
words an increase in ammonia with increasing altitude.

Based on our nominal deep ammonia abundance, the ammonia abundance increases
by 47 ppm per pressure scale height with a large spread as indicated by the
14% and 86% quantiles. Nevertheless, we can exclude a decrease with altitude,
as the distribution does not include positive values. Even if we relax the
deep ammonia abundance to our 10-$\sigma$ combination, that is increasing the
deep ammonia abundance by 41ppm with respect to the nominal model, the
negative gradient remains, however, the distribution now also includes
positive values. If we wanted to force a positive mixing gradient in the
equatorial zone, we have to increase the deep ammonia abundance to 428ppm. To
understand how the gradient correlates with the deep ammonia abundance, we
expand our retrieval to also determine the deep ammonia abundance. We show in
Figure 22 two cases, the scenario where we allow both ammonia and water to
vary, and one scenario where we force the mixing gradient to be positive.

(a) Mixing gradient with unconstrained deep ammonia abundance

(b) Forced positive mixing ratio

Figure 22: Vertical mixing gradient when allowing the deep ammonia abundance
to vary for 3 deg N.

Even when we let the deep ammonia vary (see Figure 22(a)) the retrieval favors
solutions where the ammonia increases with altitude at the 2-$\sigma$ level.
The retrieved ammonia abundance is in line with the results that we showed in
Section 5.1, confirming our previous analysis. When we force the retrieval to
only explore positive mixing gradients (Figure 22(b)) the solutions tend
towards a uniform mixed column (gradient 0) with the deep ammonia abundance at
2.9$\times$ solar, outside of our limits for ammonia. Lastly, we further
relaxed the model and allow for two mixing regimes, with a break point between
the two. We compute the fit of each model by comparing the $\chi^{2}$ for the
MCMC solutions in Table 3. We can see that no other parameterization of the EZ
can outperform our simple model of having the deep ammonia fixed with an
increase in ammonia with altitude.

Table 3: $\chi^{2}$ values for the various models in the EZ Model description | $\chi^{2}$ |  | $\chi^{2}$
---|---|---|---
Fixed deep ammonia |  | Other models
2.47 | 11.5 | Free deep ammonia | 11.9
2.69 | 11.5 | Two mixing regimes | 11.9
2.89 | 13.2 |  |
2.47 + positive mixing | 14.3 |  |

###  6.5 Water abundance

The value for the water abundance in Jupiter’s troposphere is poorly
constrained, yet it is a key parameter to distinguish different formation
scenarios (Hersant2004) and to understand the dynamics of the atmosphere
(Ingersoll1998). One of the main goals of the Juno mission is therefore to
determine, or at least to better constrain, the water abundance (Janssen2017).
This is not an easy task, as the opacity in Jupiter’s atmosphere is dominated
by ammonia gas; only observations below $\sim$1.5 GHz are sensitive to opacity
by water vapor. Furthermore, opacity by the clouds themselves is another
factor that can affect the low frequency observations (e.g., dePater2005).
Previous observations to constrain the water abundance have both focused on
small sections of the planet (Li2020; Wong2004) or only found indirect
evidence through the use of lower limits (i.e., measurements above the base of
the water cloud, Bjoraker2018). We collected all previous observations for
water and display them in Table 4 to provide reference values. In an effort to
improve our knowledge on the water abundance, we ran a suite of MCMC
retrievals, where we simultaneously constrain the ammonia and the water
abundance profiles. While the various parameters for the ammonia abundance are
well constrained by the various frequencies, the water abundance only
corresponds to the signal at the longest wavelength (see Figure 25). Due to
the low information content in the data, we opt to try and fit the simplest
model for the water abundance, that is a constant mixing ratio until the
condensation point of water vapor, from where on the water profile follows the
saturated vapor curve. Due to the importance of this measurement, we define a
few criteria for a robust measurement of the water abundance. First, the
measurement must be consistent within a climatological zone. We define the
locations of the climate zones by the locations of the jets as identified by
the motion of cloud features (Ingersoll1981; Tollefson2017), and we expect
that the measurement is consistent across these regions. This assumption is
based on the ammonia abundance profiles (see Figure 9), where ammonia varies
only slightly within a zone or belt. Secondly, the uncertainties from the
retrieved parameters influence the results through the $\chi^{2}$
minimization, that weighs each data point by the fitted uncertainty. We have
shown how the results for the various atmospheric quantities depend on the
uncertainties (Figure 20). We, therefore, require the water abundance to be
robust to the choice of uncertainties. Lastly, the water abundance should be
robust to adding further information, such as VLA observations.

Table 4: Measurements of Water on Jupiter Location on Jupiter | Pressure | Derived water abundance (solar) | Reference
---|---|---|---
Radio hot spot | 17.6-20.9 bar | 0.3 $\pm$ 0.15 | (Wong2004)
Great Red Spot | 5 bar | >1.1 | (Bjoraker2018)
Equatorial zone | >6 bar | 2.5 (+2.2, -1.6) | (Li2020)
Figure 23: Flowchart showing the various combination of data and uncertainties
used for the retrieval of the water abundance at the low latitudes on Jupiter.

We have selected 6 different combinations of data and uncertainties (Figure
23) to retrieve the water abundance. We start out by retrieving the water
abundance based on the orbit average of the Juno mission (PJ1-12). We use the
3-$\sigma^{fit}$ values which formed the basis for our ammonia results.
Ignoring large scale global variations, we can reduce the fitting
uncertainties (Equation 8) by a factor of $\sqrt{n_{PJ}}$, the number of PJs
($n_{PJ}$ = 9) with sufficient coverage to fit the nadir brightness
temperature and the limb-darkening coefficient. This reduces the nominal fit-
uncertainties to 1-$\sigma^{fit}$. Lastly, taking full advantage of the data
set, we can average over 5 degrees of latitude to further reduce the
uncertainties by a factor $\sqrt{5}$ to $\sim$ 0.5-$\sigma^{fit}$. Averaging
greatly reduces the uncertainty for the limb-darkening coefficient
$\sigma^{fit}_{p}$, while the absolute calibration errors dominate the
temperature measurements ($\sigma^{cal}_{T}\gg\sigma^{fit}_{T}$). As a
consequence, averaging not only reduces uncertainties, but it also forces the
retrieval to favor the limb-darkening measurement over the brightness
temperature measurement, as suggested at the time by Janssen2005. Since the
limb-darkening measurement in Channel 2 has such a large influence, we limit
the ammonia abundance to stay within our limits $2.47^{+0.22}_{-0.13}\times$
solar defined in Section 5 to avoid unrealistic solutions. In addition to the
orbit average Juno data, we include the VLA longitude-smeared brightness
temperature observations (dePater2016). While only the longest frequency
observations of the VLA probe below the water condensation level, the
observations are relevant to constrain the top part of the atmosphere. The VLA
data help constrain the top mixing pressure ($P_{t}$) and the mixing gradient
($H_{mix}$), which affect the structure even below the VLA sensitivity. We
also test if water alone could explain the extra source of opacity at the low
frequencies by adding the limb-darkening coefficient of Channel 1 to our
retrieval.

Before we show the results of our ammonia and water retrievals, it is
important to understand the effect of the water abundance on the observables:
nadir brightness temperature and the limb-darkening coefficient. dePater2005
showed the sensitivity of radio spectra to the water vapor and clouds, while
Janssen2005 examined the possibility of detecting water through measurements
of the limb-darkening parameter. We combine those two to understand the effect
of water vapor on the radio spectra and the limb-darkening. To this purpose,
we have forward modeled the atmosphere (ignoring the effects of cloud
particles) based on a range of vertical ammonia and water vapor profiles that
are shown in Figure 24. We have selected 6 profiles from the last 1000 walkers
of our MCMC-simulations (see Figure 23 \- VLA & Juno 1$\sigma$). Profiles
1,2,5, and 6 were chosen to be close to region within which 68% of the MCMC-
results fall, hence are all consistent with the observations. Profile 3
corresponds to the peak of the posteriori distribution for the deep water and
deep ammonia abundance, and profile 4 indicates the maximum ammonia in the
last 1000 walkers. Based on these profiles we can compute the brightness
temperature as a function of emission angle, and compare these results to the
Juno observations and their uncertainties. In Figure 25 we plot the brightness
temperature corresponding to the selected profiles. The largest effect on the
nadir brightness temperature is seen in Channel 1, but the difference between
the model results becomes smaller with increasing frequency. The offset
between the forward models and the observations in Channel 2 is consistent
with our previous analysis (Figure 10). Despite the large range in water
abundance, the effect on the brightness temperature is smaller than we can
distinguish based on the calibration uncertainty shown by the light shading.
These results are consistent with the expectations laid out by dePater2005 and
Janssen2005. The right-hand-side shows the effect of the water abundance on
the limb-darkening parameter, where we normalized all observations to the Juno
nadir brightness temperature (shown by the solid gray line with the fit
uncertainty indicated by the shading). The largest effect can be seen in the
Channel 1 observations, where we can clearly distinguish between the high
water abundance (profile 1-3) and the low water abundance (profile 4 and 5).
In Channel 2, the differences between the various profiles can no longer be
distinguished within the 3$\sigma^{fit}$ uncertainty, while in Channel 3 most
solutions are degenerate within the 1$\sigma^{fit}$ uncertainty level. Due to
the unknown sources of opacity in Channel 1 limb-darkening results, our best
chance of retrieving the water abundance is based on the Channel 2 limb-
darkening parameter, especially if we can reduce the uncertainties.

Figure 24: A subsample of vertical abundance profiles for ammonia and water
consistent with the observations based on our MCMC results (see Figure 23 \-
VLA Juno 1$\sigma$). We have taken 6 representative profiles and use these as
inputs for the forward model (see Figure 25). There are two classes of
solutions we select from: high water abundance (profiles 1-3) and low water
abundance (profiles 4-6). The colorbar assigns a unique color for each
combination of the deep ammonia and the deep water abundance. The contour
lines correspond to the posteriori distribution, denoting the region within
which 32% and 68% of the MCMC results fall. We note that the water abundance
spans an order of magnitude, while the variation in ammonia is much smaller.

Figure 25: The forward modeling results of the various profiles as shown in
Figure 24 are plotted as dashed lines compared to the observations (gray
lines). The left-hand-side shows the effect on the brightness temperature as a
function of emission angle for the various profiles, compared to the MWR
observations and the absolute calibration uncertainty indicated by the
shading. The right-hand-side shows the effect on the limb-darkening, where the
individual simulations are normalized to agree to the Juno nadir observations.
The shaded region corresponds to the fit uncertainty $\sigma^{fit}$ (n=1,3).
The large absolute calibration uncertainty $\sigma^{cal}$ does not allow
distinguishing between the various profiles based on the nadir brightness
temperature. The limb-darkening is more sensitive to the water abundance with
Channel 1 showing signals that could be distinguished given the uncertainty.
The signal in Channel 2 can only be seen when we chose n=1 (i.e.
1$\sigma^{fit}$) for the our uncertainty value.

Due to computational limitations, we focus our analysis on the SEB, the EZ,
and the NEB, the regions with the greatest variations in ammonia, and thus the
region where we expect the largest differences in the water abundance. We
retrieve the water abundance (as well as all other parameters in our model)
for 7 different regions within these zones and belts. For the 3-$\sigma$ and
1-$\sigma$ results the distribution corresponds to the center of the box
(-15∘,-11∘,-3∘,0∘,3∘,11∘,15∘), and the 0.5-$\sigma$ combinations correspond to
the 5 degree average centered on that latitude.
Figure 26 compiles all the results for the water abundance, where we display
the a-posteriori distribution based on the MCMC results for the various
combinations of data and uncertainties. Revisiting our assumptions for a
robust water signal, we see that none of the three criteria are met. The
retrieved water abundance changes drastically between adjacent slices in the
atmosphere, especially relevant in the EZ where it varies from less than 1
$\times$ solar at -3∘ to a value around 9 $\times$ solar at +3∘. In the NEB
(11∘), where we probe the deepest levels and expect the most reliable results
for the water abundance, MCMC cannot fit a single abundance, instead it is
fully unconstrained regardless of the choice of data and uncertainty. The most
promising results are in the SEB, where we find weak evidence for an elevated
water abundance compared to the ammonia results, nevertheless the analysis
across all the different latitude slices does not inspire a lot of confidence
in this signal. Interestingly, we can match the extra opacity in Channel 1
with a well-constrained water abundance of 17$\times$ solar, and the results
are consistent between the NEB and EZ. Nevertheless, we know that there are
other absorbers (dePater2005; Bellotti2018) at these frequencies, and further
experiments are required to understand this signal.

For completeness we also show the results for the ammonia abundance that we
retrieve simultaneously with the water abundance in Figure 27. The ammonia
abundance is consistent with the expectations we laid out for the water
abundance. Results look consistent across latitudes, the impact of the
uncertainties is notable, however results are consistent for various
combinations of data.

Based on these results, the water abundance on Jupiter remains elusive. The
biggest factor to improve these results are understanding how water and
ammonia absorption behave at large pressures, allowing both Channel 1 and
Channel 2 data to be more reliable. This also includes understanding the extra
absorption that we see in Channel 1, where we can only match the data with an
order of magnitude larger water abundance than was seen in the past. At these
wavelengths the contribution of Alkalis also becomes significant
(Bellotti2018), and shortwards of $\sim$5 GHz scattering/absorption by water
droplets starts to affect radio waves (dePater2005).

Figure 26: A-posteriori distributions of retrieved water abundances based on
various combinations of data and uncertainties as shown in Figure 23, with
each color corresponding to a different data combination. The 7 boxes
correspond to the 7 regions we selected to compute the abundance, with at
least two boxes per zone or belt. Even the combinations with the strongest
signal (yellow and purple) cannot constrain a consistent water signal across a
climatological zone. The differences within the EZ are too large to ascribe to
atmospheric variations, and are mostly due to the low signal strengths given
the current data available. Figure 27: Same as Figure 26 but for the ammonia
abundance retrieved simultaneously. We can see that the ammonia signal is
reproducible with various data combinations, and consistent between different
slices in the EZ, SEB and NEB.

##  7 Conclusions

Radio observations allow to disentangle some of the mysteries that Jupiter
still holds for us. A combined analysis of VLA and Juno observations allow us
to derive the vertical structure of the atmosphere; however, this is not
unique due to the fact that radio waves originate from a range of pressures.
We introduce disequilibrium factors into a thermo-chemical equilibrium model
(Weidenschilling1973; Atreya1985) that describe the distribution of trace
gases and use these factors in our inversion technique. Fitting the required
quantities, however, also requires a good understanding of the fidelity of the
data product. For this reason, we developed a data reduction pipeline, that
obtains the brightness temperature and limb-darkening coefficient along with
their uncertainties. Coupling these two, results in a distribution of trace
gases (see Figure 9), the sensitivity of our solutions to the fit
uncertainties (see Figure 20) and our certainty in the derived atmospheric
distribution (see Figures 16 and 17).

Our approach recovers an atmosphere that confirms the structure that was seen
by the VLA in the past (dePater1986; dePater2019) and by Juno (Li2017),
however, we also find some remarkable differences. We find that the deep
ammonia abundance of $366.5^{+32.8}_{-19.9}$ ppm ($2.47^{+0.22}_{-0.13}\times$
solar abundance) is lower than the Galileo (Wong2004) and VLA results
(dePater2016; dePater2019, which were based on the Galileo results), however
slightly higher than the Juno results (Li2017). The atmosphere only approaches
a state of being well-mixed at pressures greater than 40 bars, sitting under a
layer between 20 and 40 bar in which the ammonia abundance gradually decreases
with increasing altitude. The atmospheric structure requires processes that
deplete the atmosphere down to such high pressures, or alternatively require a
mechanisms that locally enrich the atmosphere, such as in the equatorial zone.
The formation of the NH4SH and ammonia cloud layers dominates the vertical
structure in the ammonia abundance in the upper troposphere. When we use the
Juno results for the deep ammonia abundance and apply it to the VLA
observations (dePater2016; dePater2019), we can recover a remarkably similar
atmosphere from the ground, with the largest difference in the northern
hemisphere, that was subject to a large storm outbreak early on in the Juno
mission (Sanchez-Lavega2017).

The uncertainties that we retrieve from the raw Juno observations allow us to
study the sensitivity of our results to these observational uncertainties.
Especially the low frequency observations are hard to fit, where the limb-
darkening coefficient and the brightness temperature favor different solutions
for the structure of the atmosphere for pressures below 10 bar. Higher up in
the atmosphere the uncertainties influence the retrieval less (see Figure 8).
The equatorial zone represents an interesting diversion from the general
atmosphere. It shows the highest ammonia abundance on the planet, with an
indication that the ammonia is increasing with altitude, in stark contrast to
the rest of the planet. This signal is robust to the chosen uncertainties (see
Figure 20) and can be reproduced even when we loosen the assumption and let
the deep ammonia vary. This structure can also be reproduced with the VLA
observations (see Figure 18).
The highly sought after water abundance is a key feature of Jupiter science
and the Juno mission. The longest wavelength observations of Juno MWR,
especially the limb-darkening coefficient, are sensitive to the water
abundance. We set out to retrieve the water abundance in the tropical regions
of Jupiter using an MCMC algorithm that combined different sets of data and
their corresponding uncertainties. While looking at an isolated case, such as
a single combination of data (only Juno data) or a single slice of latitude
(only the EZ), we can retrieve a water abundance. Unfortunately, the signal is
not robust to different combinations of data, or different levels of
uncertainty, i.e., the water abundance is only marginally constrained by the
current set of measurements.
So how can we advance our understanding of the deep atmosphere and the
corresponding water abundance? First and foremost, we need to investigate the
signal in Channel 1 which so far has been mostly ignored. With the differences
between the various ammonia absorption coefficients at high pressures
(Hanley2009; Bellotti2016; Bellotti2017) alone, Jupiter’s deep atmosphere
remains poorly understood. Furthermore, at these long wavelengths the
absorption by water clouds (dePater2005) and other species such as alkalis
(Bellotti2018) becomes significant, both of which are poorly understood and
often ignored in the RT models. Both lab measurements at high pressures and
temperature and ab-initio calculations are essential to fully model these
contributions to the opacity. The other big unknown is the temperature
structure below the clouds. While infrared observations can probe down to 0.8
bar, the structure below is largely unconstrained, and usually assumed to be
adiabatic (wet or dry). It is quite impossible to determine the temperature
structure from remote observations; understanding the dynamics can help
construct temperature profiles that are more representative of an atmosphere
controlled by intermittent convection.

## Appendix A - Beam-Convolved Emission Angle

As explained in Section 2, for each observation we compute the beam-convolved
emission angle to account for the large range of emission angles within the
beam. Since we are interesting in the relationship between the emission angle
and the brightness temperature we must assure that the emission angle is
representative of the observations. Especially at low emission angles, the
extended beam covers a wide range of emission angles, and the boresight
emission angle is not an adequate indicator for the observation. Figure 28
shows the relationship between the two quantities for a full rotation of the
spacecraft. The large beam-size of Channel 2 makes the effect more pronounced
than the observations in Channel 5.

(a) Channel 2

(b) Channel 5

Figure 28: Relationship between the beam-convolved emission angle and
boresight emission angle. The different colors indicate the integration region
in units of the HPBW. Further away from the beam center, the beam sensitivity
is dropping rapidly, and the contribution to the averaging becomes small.
Therefore, the results converge for sufficiently large integration regions.
Our final results are based on 2.5 HPBW, which contains 99% of the emission.

## Appendix B - Synchrotron filter

To best understand the spatial synchrotron filter, in Figure 29 we show a
region of the planet where the radiation from the main synchrotron belt is
leaking in when looking towards the equator.

Figure 29:  Same as Figure 3 but for a latitude where the synchrotron
contribution is noticeable. Note that the there are two branches in Channel 1
and 2 that have a large difference in the temperature for the same emission
angle. The increase temperature is due to the synchrotron leaking in through
the side lobes. Our spatial restriction can handle the synchrotron
contribution well, ignoring the upper, contaminated branch of the
observations.