Jupiter's Multi‐Year Cycles of Temperature and Aerosol Variability From Ground‐Based Mid‐Infrared Imaging

We use a long‐term record of ground‐based mid‐infrared (7.9–24.5 μm) observations, captured between 1984 and late 2019 from 3‐m and 8‐m class observatories (mainly NASA's Infrared Telescope Facility, ESO's Very Large Telescope, and the Subaru Telescope), to characterize the long‐term, multi‐decade variability of the thermal and aerosol structure in Jupiter's atmosphere. In this study, spectral cubes assembled from images in multiple filters are inverted to provide estimations of stratospheric and tropospheric temperatures and tropospheric aerosol opacity. We find evidence of non‐seasonal and quasi‐seasonal variations of the stratospheric temperatures at 10 mbar, with a permanent hemispherical asymmetry at mid‐latitudes, where the northern mid‐latitudes are overall warmer than southern mid‐latitudes. A correlation analysis between stratospheric and tropospheric temperature variations reveals a moderate anticorrelation between the 10‐mbar and 330‐mbar temperatures at the equator, revealing that upper‐tropospheric equatorial temperatures are coupled to Jupiter’s Equatorial Stratospheric Oscillation. The North and South Equatorial Belts show temporal variability in their aerosol opacity and tropospheric temperatures that are in approximate antiphase with one another, with moderate negative correlations in the North Equatorial Belt and South Equatorial Belt changes between conjugate latitudes at 10°–16°. This long‐term anticorrelation between belts separated by ∼15° is still not understood. Finally we characterize the lag between thermal and aerosol opacity changes at a number of latitudes, finding that aerosol variations tend to lag after thermal variations by around 6 months at multiple latitudes.

Plain Language Summary Jupiter's atmosphere displays a wide variety of perturbations in its temperatures, clouds and aerosols.In this study, we use a large set of ground-based observations captured in the mid-infrared between 1984 and 2019 to characterize long-term changes in the temperatures and aerosols.This long-term analysis show a number of cyclic disturbances, and allows us to distinguish between seasonal and non-seasonal changes in Jupiter's atmosphere.In particular, we observe that the northern mid-latitudes above 30° are continuously warmer than their counterpart latitudes in the south at 10 mbar pressure level (the stratosphere), potentially due to differences in the polar haze in Jupiter, which extends to lower latitudes in the north compared to the south.Additionally, our study reveals for the first time that the thermal oscillation present in Jupiter's equatorial stratosphere at the 10-mbar pressure level (known as Jupiter's Equatorial Stratospheric Oscillation) is also observed to descend to higher pressures (330 mbar), meaning that it is not confined to the stratosphere.Finally, we also discuss the lag between temperature and aerosol changes at diverse latitudes to try to identify the mechanisms responsible for the different atmospheric disturbances observed on Jupiter.
• Ground-based multi-wavelength images are used to compute stratospheric and tropospheric temperature and tropospheric aerosol opacity maps • Results reveal that upper-tropospheric equatorial temperatures are coupled to Jupiter's Equatorial Stratospheric Oscillation • Stratospheric temperatures at 10 mbar show a permanent hemispherical asymmetry at mid-latitudes, with northern mid-latitudes overall warmer than their counterparts

Supporting Information:
Supporting Information may be found in the online version of this article.
Ground-based 5-μm time-series spanning 34 years, sensing thermal emission from the 2-7 bar region modulated by overlying clouds (Bjoraker et al., 2015;Giles et al., 2015), enabled the discovery of a rare cyclic disturbance in Jupiter's Equatorial Zone (EZ,±7° latitude).A comparison between the periodic brightening found in the EZ at 5 μm and long-term changes at visible wavelengths (Antuñano et al., 2018), showed that these disturbances occur contemporaneously with (a) visible reddening events of the EZ (Rogers, 1995), (b) a decrease of the ultraviolet reflectivity at 410 nm (Simon-Miller & Gierasch, 2010), and (c) a decreased aerosol opacity at 600-800 mbar (Antuñano et al., 2020), suggesting a clearing of the NH 3 clouds, and potentially hinting at changes in the ammonia upwelling from the deeper tropospheric levels.A long-term analysis of Jupiter's zonal winds by Simon-Miller and Gierasch (2010), using observations captured by the HST between 1994 and 2008, and later extended to 2017 and 2019 by Tollefson et al. (2017) and Wong et al. (2020) respectively, revealed similar periodic changes in the equatorial zonal-winds at cloud level, with faster winds observed near-contemporaneously to the 5-μm brightness change (Antuñano et al., 2018), potentially due to tracking faster clouds located deeper in the troposphere during the disturbances.
These long-term studies of Jupiter's atmosphere also reported non-periodic changes in the zonal winds (Wong et al., 2020) related to large convective outbreaks at the SEB and NTB, helping researchers to predict future events, and revealed an anticorrelation of the 5-μm brightness changes between the NEB and the SEB, suggesting a potential interhemispheric connection between Jupiter's prominent equatorial belts (Antuñano et al., 2019).
In this study, we aim to extend the investigations of Orton et al. (1991Orton et al. ( , 1994)), Antuñano et al. (2019Antuñano et al. ( , 2021)), and Orton et al. (2023) to a larger number of mid-infrared wavelengths to explore Jupiter's temperature and aerosol variability over long spans of times, from the NH 3 -ice clouds to the stable mid-stratosphere, enabling 10.1029/2022JE007693 3 of 34 us to separate seasonal from non-seasonal variability.We aim to show that, although Jupiter's troposphere and stratosphere exhibit myriad dynamical phenomena, there exist quasi-periodic patterns (e.g., Rogers, 1995) that may aid in future predictions of planetary-scale changes to the banded structure.We will use a large dataset of ground-based mid-infrared observations spanning 36 years  at eight wavelengths between 7.9 and 24.5 μm, allowing us to derive latitude-and time-resolved temperature profiles and aerosol distributions via spectral inversion.This extends the analysis of Orton et al. (1991Orton et al. ( , 1994) ) and Simon-Miller et al. (2006), which investigated tropospheric and stratospheric thermal variations between 1978 and 1992, and was the original motivation for our longer-term study.It also uses a larger filter set than the study by Orton et al. (2023) over a similar time frame.
The paper is organized as follows.In Section 2, we describe the observations, image reduction process, and data processing techniques applied in this study.Long-term variations in brightness temperature at each wavelength are explored independently in Section 3 and are compared to previously reported meteorological events in Jupiter.
Optimal estimation inversions, performed to characterize the zonally averaged stratospheric and tropospheric temperatures and aerosol opacity between 1983 and 2019, are described in Section 4. The results are shown in Section 5 and discussed in Section 6.Finally, we summarize our conclusions in Section 7.

Observations and Image Reduction
This study uses ground-based mid-infrared images captured between 1983 and 2019 at wavelengths between 7.9 and 24.5 μm using a number of different instruments (Table 1): BOLO-1 (1983BOLO-1 ( -1993)), the Mid InfraRed Array Camera (MIRAC, 1993(MIRAC, -1999)), MIRLIN (1996MIRLIN ( -2003)), and the Mid-Infrared Imager and Spectrometer (MIRSI, 2003(MIRSI, -2011) instruments mounted at the 3-m NASA Infrared Telescope Facility (IRTF) in Hawai'i; the Very Large Telescope (VLT) Imager and Spectrometer for mid-InfraRed (VISIR (2006(VISIR ( -2011(VISIR ( , 2016(VISIR ( -2018) ) instrument mounted at the 8-m VLT in Chile; and the Cooled Mid Infrared Camera and Spectrometer (COMICS, 2005(COMICS, -2019) ) instrument on the 8-m Subaru Telescope in Hawai'i.A summary of the wavelength range, detector, and plate-scale of these instruments is given in Table 1, and more detailed information on these instruments can be found in the references given in the table.Figure 1 summarizes all the observations used in this study as a function of wavelength and full information on the exact dates, instruments, wavelengths and central meridian longitude coverage is given in Table S1 in Supporting Information S1.As described later in Section 2.2.1, wavelengths shown in Figure 1 are approximate due to differences between instrument filters, and correspond to the wavelengths used during the calibration process to enable coverage over longer time spans.We note that mid-infrared observations at a wider range of wavelengths (e.g., 9.8 μm, 11-12 μm) and captured by other instruments (e.g., MICHELLE and T-ReCS mounted on the Gemini telescope in Hawai'i) are also available.However, these cover shorter time-frames (from months to less than 20 years) or were captured at epochs that were very scattered, impeding a long-term analysis.All 18.72 and 20.5 μm data used in this study are the same as those used in (Orton et al., 2023).
All images used in this study, except those captured by the BOLO-1 instrument, are reduced using the Data Reduction Manager pipeline described by Fletcher et al. (2009).This process consists of subtracting the sky and telescope background from each image to be able to detect Jupiter's weak emission via a chopping-nodding  technique, correcting detector non-uniformities and bad-pixels by the application of flat-fields and bad-pixel masks to each image, geometrically calibrating (i.e., navigating) the data by limb-fitting, and projecting them into cylindrical maps at the desired spatial resolutions.As BOLO-1 images were captured using raster-scan techniques sampling the planet in a 1″ regular grid, the previous steps cannot be followed to process these data.Here, we use the reduced and projected BOLO-1 images as described by Orton et al. (1991Orton et al. ( , 1994)).
Jupiter images captured by the VISIR instrument are partially obscured by the negative beam when chopping the telescope for sky subtraction.This is due to the 38 × 38″ field of view of the VISIR instrument and the maximum chopping amplitude of 25″ of the VLT, which is smaller than Jupiter's typical angular diameter of 40″.In this study, no attempt is made to correct the obscured regions and instead, we remove them from each image by cropping the affected latitudes before the navigation step.Additionally, images captured by the VISIR instrument from 2016 onward display a pattern of vertical and horizontal stripes across the detector that cannot be entirely removed by a flat-field.Here, we remove the striping pattern following the steps in Donnelly (2021), which corrects independently the horizontal and vertical stripes by using a Gaussian smoothing for the former, and low-pass and high-pass filtering for the latter.This correction, however, does not entirely remove the central stripe where the two halves of the VISIR detector meet.For this reason, we perform an additional step, which consists of removing the central burst by fitting it to a Gaussian function and subtracting it from the brightness profile.Finally, due to the detector sensitivity and the brightness of the telluric atmosphere at 17.65 and 18.72 μm, the field of view of the VISIR observations between 2016 and 2018 at these wavelengths is smaller than at other wavelengths, meaning that Jupiter images are windowed, cutting off the high-latitudes.
Examples of reduced 8.6-μm Jupiter images captured during each of the four decades analyzed in this study are shown in Figure 2a, showing the evolution of the quality of ground-based mid-infrared observations.Figures 2b-2i shows examples of fully reduced and mapped images of Jupiter captured in May 2018 in the eight filters shown in Figure 1 by the VISIR (left column, 7.9-13.04μm) and COMICS (right column, 17.6-24.5μm) instruments, and correspond to some of the highest quality ground-based mid-infrared observations of Jupiter obtained to date.The filters shown in this figure probe the following pressure levels: 7.9 μm probes the stratosphere at around 10 mbar; 8.6 and 10.77 μm probe temperatures, ammonia and aerosol near the NH 3 cloud level at around 600-800 mbar; and the upper troposphere and tropopause are probed by the 13.04-24.5μm filters.VISIR images shown in this figure have been corrected to remove the artificially obscured regions mentioned above, by adding back in the subtracted flux (Antuñano et al., 2020).This correction, however, leaves a "residual arc" (shown in Note that the wavelength of the observations might not exactly match the wavelength annotated here.See Section 2.2.1 for information on the wavelength shifting performed in this study.A summary of these data is given in Table 1 and a full description is available in Table S1 of Supporting Information S1. 10.1029/2022JE007693 6 of 34 black in Figure 2) at the edge of the obscured regions that is more apparent at wavelengths with limb-brightening or weak limb-darkening (Donnelly, 2021).As mentioned above, this correction is only performed to show the quality of the observations and it is not used in our subsequent quantitative analysis.

Radiometric Calibration
Each image and each BOLO-1 map of the disk is radiometrically calibrated following Antuñano et al. (2020) and Orton et al. (2023), with Q-band (i.e., 17.2-24.5μm) and N-band (7.9-13.04μm) images scaled in radiance to match Voyager IRIS and Cassini CIRS profiles, respectively, within ±50° latitude of the equator.CIRS and IRIS were both Fourier transform spectrometers calibrated in space by viewing a warm black body and cold deep space, and using the difference to calibrate the resulting interferograms.Although separated in time by decades, the CIRS and IRIS spectra of Jupiter are self-consistent in spectral regions of overlap.They remain our most reliable reference sources for in-space calibration.Scaling our long-term observations to Cassini CIRS and Voyager IRIS assumes that Jupiter's overall brightness has not varied significantly with time, resulting in this technique being insensitive to any global-scale change in the mid-infrared brightness over these long spans of time.However, relative changes in brightness from latitude to latitude are robust.The selected scaling region differs from Antuñano et al. (2020), where data were scaled to the average radiance over ± 20°-60° latitude to avoid the equatorial and off-equatorial latitudes, which are more variable than the higher latitudes.However, as mentioned above, VISIR Jupiter images between 2016 and 2018 at 17.6 and 18.7 μm are windowed, cutting off the high-latitudes (i.e., above ∼50° latitude) and preventing the use of the same scaling region as in Antuñano et al. (2020).The choice of the scaling region used in this study was made after several experiments, and was found to be the best method to calibrate the data without compromising temporal changes at any latitude.
As the instruments used in this study have slightly different filters, BOLO-1, MIRAC, MIRLIN, MIRSI and COMICS filters are shifted in wavelength to match those of the VISIR instrument and COMICS 19.5-μm filter (all in Figure 1).This is done before scaling the data to match IRIS and CIRS, to treat all filters equally over the entire time series (e.g., shifting the MIRLIN 18.67-μm filter to treat it as the VISIR 18.72-μm filter).Additionally, due to the higher spatial resolutions of VISIR and COMICS observations compared to the observations acquired with the 3-m IRTF telescope, VISIR and COMICS images are smoothed before the calibration to match the spatial resolution of MIRSI observations and are used throughout the subsequent analysis.We follow Orton et al. (1994Orton et al. ( , 2023) ) in scaling the 17.80-μm images captured by the BOLO-1 instrument in 1980BOLO-1 instrument in -1982BOLO-1 instrument in and 1984BOLO-1 instrument in -1985 as if it was a 18.72 μm filter, as the difference in the peak contribution functions at these wavelengths are relatively small and the Q-band BOLO-1 filter was extremely broad (16-26 μm according to the IRTF photometry manual from April 1988).This process ignores potential small variations of the radiance due to the slightly different filters.However, these variations are expected to be within the absolute-calibration uncertainty, which is as large as 14% in the Q-band data and around 4%-5% in the N-band data (Antuñano et al., 2020).

Zonal-Mean Radiance and Smoothed Radiance Profiles
Latitudinal radiance profiles as a function of date and wavelength are obtained by longitudinally averaging the radiance of each image within 30° longitude of the minimum emission angle in 1° latitudinal bins.For the chop-obscured regions of VISIR images, only latitudes where the obscured region does not fall within 15° longitude east and west of the central meridian were considered.This choice is introduced to ensure that an adequate region around the central meridian is used for the zonal-mean radiance.Images where the Great Red Spot was centered on the central meridian were avoided to exclude anomalous regions that do not represent the zonal radiance.We note that other anomalous regions such as NEB plumes, 5-μm hot-spots are relatively well distributed in longitude and therefore, could still be representative of the zonal mean.In addition, in the CIRS or IRIS data, as well as in ground-based observations in the N-band, 5-μm hot spots do not appear any more prominent than warm areas elsewhere.
Potential longitudinal variability at each latitude and wavelength is given by the standard deviation of the zonalmean radiance.This is usually an order of magnitude smaller than the absolute calibration uncertainty in the Q-band, but can be as large as the absolute calibration error in most of the N-band filters.Furthermore, the longitudinal variability at 8.59 μm (see Figure 2), sensing tropospheric aerosol contrasts, can be up to an order of magnitude larger than the calibration uncertainty.In this study, we include the longitudinal variability in the measurement errors used in our spectral retrievals (see Section 4), adding the standard deviation of the zonal average in quadrature to the absolute calibration uncertainty.Zonal-mean brightness temperature maps, converted from the zonal-mean radiance, are shown in Figure 3 as a function of time and latitude, showing the temporal and latitudinal variability of Jupiter at different wavelengths between 1983 and late 2019.This variability is fully described in Section 3.
To obtain radiance profiles that better represent the full dataset, and to fill in gaps when instruments were unavailable, zonally averaged radiance profiles are smoothed as a function of date for each latitude and wavelength.This is performed using a Savitzky-Golay (SG) smoothing filter (Savitzky & Golay, 1964), which is defined as a moving average weighted by a polynomial of a certain order.This technique allowed us to fit the data with lower uncertainty compared to simple smoothing and interpolation of the data.In this study, a 24-point wide (covering 1,440 days) and fourth-order polynomial SG filter is used.This technique requires data to be on a regular grid, so zonal-mean profiles are first linearly interpolated using a 60-day window in 1° latitudinal bins.Larger window sizes resulted in excessively smoothed profiles, while smaller windows showed an artificially wavy profile.We tested the sensitivity of our results to the chosen smoothing function and found that the Lomb-Scargle periodograms remain invariant for window sizes that better reproduce the observed trend (e.g., ranges between 1,000 and 1,800 days, see Figure S1 in Supporting Information S1).The SG smoothing is repeated 200 times for each latitude and wavelength, using randomized values of the zonal-mean radiance within the estimated measurement errors to consider the longitudinal variability and the uncertainties introduced during the radiometric calibration.
Finally, the 200 smoothed profiles for each latitude and wavelength are averaged together at each date to obtain averaged and smoothed radiance profiles in Figure 4 (with the corresponding standard deviation).This technique is identical to that used in Orton et al. (2023).
These smoothed time series are used to analyze the temporal variability of the brightness temperature at each wavelength, and also to retrieve temperature and aerosol opacity maps on a regular temporal grid (see Sections 3 and 4).Examples of average smoothed brightness temperature profiles at 8.6 μm and two different latitudes (red solid lines), compared to the measured zonal-mean brightness temperatures (black dots) as a function of time, are shown in Figure 4.This demonstrates that the average smoothed profiles adequately reproduce the trends seen in the observations, enabling us to represent the data on a regular temporal grid.These smoothed profiles also reduce our sensitivity to outliers, where discrete meteorological features may influence our approximate zonal mean on each date.
In particular, we use the "Scargle" function written in IDL, setting a false-alarm probability of 0.02 (or 98% significance) and a minimum of 1 year periodicity.The uncertainty of the obtained periods is assumed to be equal to the FWHM of the power spectrum peak.Results are shown in Sections 3 and 5, where only periods with a significance higher than 98% are represented.

Cycles of Variability in the Mid-Infrared
Brightness temperature anomaly maps between ±48° latitude, computed by subtracting the average brightness temperature over all latitudes and dates from the smooth brightness temperature profiles at each latitude, are shown in Figures 5 and 6 as a function of date.This technique enables investigation of potential temporal changes with respect to an average state of Jupiter, and provides a clearer view of the changes than Figure 3.The temporal variance of the brightness temperatures and Lomb-Scargle periodograms (Scargle, 1982) are also presented in Figures 5 and 6, showing the regions of highest variability and the timescales of these changes, respectively, at each of the eight wavelengths analyzed in this study.Brightness temperature anomaly maps show remarkable cyclic activity, which will be discussed in the following subsections.

Jupiter's Equatorial Stratospheric Oscillation
Jupiter's Equatorial Stratospheric Oscillation (JESO, also known as the Quasi-Quadrennial Oscillation or QQO, Leovy et al., 1991;Orton et al., 1991) can be observed at 7.9 μm in Figure 5a as a warm and cold temperature pattern at the equator, with the off-equatorial latitudes at ∼12°-14° changing in anti-correlation to the equator (Cosentino et al., 2017).Our 7.9-μm brightness temperatures also show the JESO disruptions -a change in the phase of JESO at the equator by 180°-from 1992 to 2008 reported by Antuñano et al. (2021), where the brightness temperatures display early brightness temperature maxima compared to the 4-4.5-year period inferred in previous studies (e.g., Simon-Miller et al., 2006).Our brightness temperature data, however, do not show the early temperature increase at the equatorial latitudes in mid 2017 suggested in Giles et al. (2020).This difference is further discussed in Section 5.
At wavelengths between 17.6 and 24.5 μm, sounding upper tropospheric temperatures at 100-300 mbar pressures (Fletcher et al., 2009), the brightness temperature at the EZ also varies in time.The Lomb-Scargle periodograms in Figures 6a-6d hint at two potential periodicities at the equatorial latitudes -a 4-year period mainly present in 17.6, 18.7, and 24.6 μm (indicated by green squares), and a 8-9-year period (indicated by blue squares) observed in all Q-band wavelengths -in agreement with Orton et al. (2023).The 4-year period observed in this tropospheric emission does coincide with the periodicity of the stratospheric JESO seen in Figure 5a, suggesting that this phenomenon could extend downwards into the upper troposphere, as predicted by Leovy et al. (1991).We discuss the potential coupling between stratospheric and upper-tropospheric temperatures in Section 6.

Equatorial Zone Disturbance in the Mid-Infrared
At 8.6 μm, sensing tropospheric temperatures and aerosol opacity at around 600-800 mbar pressure level (Fletcher et al., 2009), Figure 7 shows a periodic change in brightness temperature (i.e., a decrease in the aerosol opacity or increase in the kinetic temperature) of 7 years at the equatorial latitudes (indicated by pink square in Figure 5) related to the EZ disturbances.These disturbances have been previously observed with a ∼7-year periodicity at the NH 3 cloud tops (Antuñano et al., 2020;Rogers, 1995); below the NH 3 clouds sounded by 5-μm observations (Antuñano et al., 2018); in the upper tropospheric haze as brightness changes at blue wavelengths (Simon-Miller & Gierasch, 2010); and in the zonal wind field (Simon-Miller & Gierasch, 2010;Tollefson et al., 2017;Wong et al., 2020).The increase in the brightness temperature in 1992, 2000, 2006-2007, and 2018-2019 match all the coloration events at visible wavelengths since 1983 (e.g., Rogers, 1995) in agreement with a thinning of the NH 3 clouds suggested by Antuñano et al. (2018).
At 10.7 and 13.0 μm (the former sensing tropospheric ammonia gas and temperature near 500 mbar, and the latter sounding tropospheric temperatures at similar pressure levels, Fletcher et al., 2009) Figure 7 shows subtle increases of the brightness temperatures of <1 K at the EZ nearly contemporaneous to the 8.6-μm increases (shown by black arrows), suggesting either small warming of the equatorial latitudes during the EZ disturbances and/or changes on the ammonia gas distribution at the ∼500-mbar pressure level (see Figures 5c and 5d).This is in agreement with the small changes observed in the brightness temperature during the 2006-2007 EZ disturbance analyzed in Antuñano et al. (2020) and the 7-year periodicity seen in Figures 5c and 5d (shown by pink squares).A correlation analysis between brightness temperature variations at 8.6 and 10.7 μm at the equator shows that changes at 8.6 μm lag 360 days behind those at 10.7 μm.Similarly, 8.6-μm changes lag 180 days behind those at 13 μm (as indicated in Figure 7, where the dashed black lines indicate the 13.0 μm brightness temperature minima).
The 8-year period found at the equator in the Q-band wavelengths is not related to the EZ disturbances, as a comparison of the variability at the EZ in Figures 6a-6d with the EZ disturbances in Figures 5b-5d shows that the changes observed in the upper troposphere are not influenced by the EZ disturbances, confirming previous findings of Antuñano et al. (2020).

Off-Equatorial and Mid-Latitudes
The NEB and SEB display the largest temporal variance in the 8.6-μm brightness temperature (see Figure 5b), indicating a highly variable aerosol distribution in these two belts.The NEB displays overall higher brightness temperatures than the SEB mainly at 8.6 μm, 10.7 and 13.0 μm, in agreement with the stronger depletion in ammonia and aerosols found at the NEB compared to the SEB (de Pater et al., 2019;Fletcher et al., 2016, Fletcher, Kaspi, et al., 2020;Fletcher, Orton, et al., 2020;Li et al., 2017).At longer wavelengths (i.e., 17.6-18.72μm), the NEB displays higher temporal variance compared to the SEB, while at 20.5 and 24.5 μm no significant differences are observed between the NEB and the SEB.The Lomb-Scargle periodogram shows a subtle ∼4.5-year period at 17°-20°N latitude at 8.6 and 10.7 μm (indicated by yellow squares), coinciding with the average time intervals between the NEB expansions (Fletcher, Orton, Sinclair, et al., 2017;Rogers, 1995) and confirming the change in tropospheric aerosols and temperatures during these events.
Finally, at latitudes higher than 20°, all wavelengths show a periodic 12 ± 2-year oscillation of their brightness temperatures (indicated by white squares in Figures 5 and 6).These variations are more prominent in the Q-band wavelengths, where the brightness temperatures are observed to change in anti-phase between the northern and southern hemispheres.The similarity between the observed periodicity and the jovian year (11.8 Earth years), together with the anti-correlation between hemispheres, suggests that these observations are dominated by changes to our viewing geometry and are showing the changing emission angle (difference between the sub-observer point and the local normal at the latitude/longitude of interest) as Jupiter moves around the Sun.A comparison between variation of the minimum emission angle between 1983 and late 2019 to 20.5-μm brightness temperature changes at 40°N, show a clear anticorrelation (see Figure 8), suggesting that the observed variations at mid-latitudes are not real and should disappear when accounting for Jupiter's weak seasons.Spectral inversions, which properly account for this effect by calculating radiance for a path at the specific emission angle, will be presented in Section 5 and discussed in Section 6.
In the next section, we describe the retrievals performed to quantify the changes in stratospheric and tropospheric temperatures and aerosol opacity during the cyclic atmospheric events described above.
In this study, spectra represented by the 5-point and 8-point image cubes are inverted independently using the radiative-transfer and retrieval code called NEMESIS (Irwin et al., 2008) to provide estimations of the stratospheric and upper tropospheric temperatures, tropospheric ammonia, and aerosol opacity between 1983 and 2019 and 1996 and 2019, respectively.In both cases, NEMESIS calculates synthetic spectra for a given atmospheric profile using pre-tabulated k-distributions for each gases, along with additional "continuum" opacity sources, such as aerosol cross sections and collision-induced absorption.NEMESIS then fits the crude 5-or 8-point spectra with a non-linear Levenberg-Marquardt method, iteratively changing the free atmospheric parameters to obtain a final optimal fit that does not deviate significantly from a prior (see Sections 4.1 and 4.2).The retrievals do not account for reflected sunlight and we assume pure thermal emission at these mid-infrared wavelengths.A comparison between the retrieved results from the 5-point and 8-point spectral cubes provides an understanding of potential degeneracies inherent in the retrievals due to the low number of wavelengths used in both cases.
Additionally, although the 8-point spectra provide more information to better constrain the retrievals, the 5-point spectral image cubes allow us to extend our study back to 1983.For both retrievals, Jupiter's atmosphere is divided in 80 regular layers in log(p) between pressures of 10 and 10 −6 bar.Descriptions of the reference atmospheric models used in each case are given in the following subsections.

Retrievals With the 5-Point Spectra
In this case, we use a reference atmospheric model where tropospheric ammonia, phosphine, the vertical temperature profile (T(p)) and stratospheric hydrocarbons (methane, ethane, and acetylene) come from a low-latitude Figure 9. Functional derivative profiles dR/dx, where R is the spectral radiance and x is the atmospheric temperature, for each of the 8 filters used in this study, showing that our dataset provides sensitivity between ∼10 and ∼600 mbar.These were calculated at nadir emission angle; the peaks would move to lower pressures with higher emission angles.
The sources of spectral line data used are listed in Fletcher, Orton, et al. (2018).The reference aerosol profile corresponds to a thick cloud of 10±5-μm radius NH 3 ice particles with base at 800 mbar, top at 400 mbar and with a fractional scale height of 0.2 times the gas scale height.However, assuming particles with a 1 ± 0.5-μm radius does not change the retrieved temperatures and gas abundances.
With only five spectral points to define the mid-IR spectrum, the retrieval suffers from significant degeneracy between tropospheric temperatures and the ammonia distribution.In particular, low brightness temperatures at 10.7 μm could result from an increased NH 3 abundance, or a decreased temperature at 500 mbar, and without the added constraint of the 8-filter retrieval we are unable to distinguish between these cases.For this reason, we model our observations using a two-step approach: (a) first, we retrieve temperature, ammonia, and aerosol opacity simultaneously by allowing the vertical temperature profile to vary while scaling the reference aerosol and ammonia distributions, holding para-H 2 fixed at equilibrium; second (b) we use the average of the previously retrieved latitudinal ammonia profile over all dates as the prior profile and retrieve temperatures and aerosol opacity simultaneously whilst assuming that the ammonia gas does not vary with time.Comparisons of the quality of the fits between method (a) and (b) (χ 2 , shown in Figure 10a) and the retrieved temperatures and aerosol opacity profiles for the Equator (Figure 10b) and for 10°N (Figure 10c) obtained from these two approaches show improvements of the fit at 10°N in 1993, 2007, and 2011, but do not show any notable differences at the equator, suggesting that either the ammonia gas distribution has not changed significantly at this latitude in the last 36 years or that these five filters alone are insufficient to be able to characterize its latitudinal and temporal variability.We believe it is the latter, as changes in ammonia distributions are expected to occur at the equatorial and off-equatorial latitudes during the EZ disturbances, NEB expansions, and SEB fading and revival events (Antuñano et al., 2018(Antuñano et al., , 2020;;Fletcher et al., 2011;Fletcher, Orton, Rogers, et al., 2017;Fletcher, Orton, Sinclair, et al., 2017), although no such change has been yet reported at the equatorial latitudes.In this study, therefore, we show and discuss the results obtained from the second approach.

Retrievals With the 8-Point Spectra
A similar reference atmospheric profile to that described above is used for the 8-point spectral retrievals, with the difference that the a priori temperature profile comes now from an average of the retrieved temperatures from the 5-point spectra between 1983 and 2019.This is done to ensure that the reference atmospheric model used represents the average thermal state of Jupiter's atmosphere during the period analyzed in this study.The a priori aerosol profile is the same as in Section 4.1 -a cloud of 10±5-μm radius NH 3 ice particles in the 400-800 mbar pressure level with a fractional scale-height 0.2 times the gas scale-height.However, unlike the 5-point spectra, using a reference aerosol profile with 1 ± 0.5-μm-radius NH 3 particles does not fit the observations at 10.77 and 13.04 μm, particularly at latitudes with high aerosol opacity.
As in the 5-point spectral analysis, we investigate potential degeneracies inherent in the retrievals by comparing the retrieved atmospheric parameters from three different retrieval approaches: Model (1) retrieves stratospheric and tropospheric temperatures and tropospheric ammonia gas distribution and aerosol opacity simultaneously by allowing temperatures to vary freely and scaling ammonia and aerosol opacity profiles; Model (2) is the same as (1) but assumes that ammonia remains unchanged with time, retrieving only stratospheric and tropospheric temperatures and tropospheric aerosol opacity; Model (3) is the same as (2) but also allows para-H 2 to vary as well (due to changes of the opacity of the S(0) and S(1) lines potentially controlling the thermal infrared continuum sensed by the Q-band filters, that is, 17.6-24.5μm Fletcher et al., 2009).
Residual maps of χ 2 , computed by subtracting the goodness-of-fit χ 2 of model 1 from those of model 2 (Figure 10d), and model 3 from to those of model 2 (Figure 10g), show that there is no significant difference on the fits performed in model 1 and 2, except at a narrow region around 10°N (Figure 10f), where ammonia is required to deviate substantially from the assumption of a spatially and temporally uniform prior.The poor fits of the observations at the southern edge of the NEB compared to other latitudes, even when our model enables ammonia to vary during the fitting process (see Figures 10e and 10f), suggests a lack of ammonia gas at these latitudes, in agreement with the large depletion of ammonia observed at the NEB in the microwave, radio and thermal infrared (e.g., Fletcher, Kaspi, et al., 2020;Fletcher, Orton, et al., 2020;Li et al., 2017).The small differences in the goodness of fits between model 1 and 2 indicates that, as with the 5-point spectral analysis, the low number of spectral points prevents us from confidently retrieving the ammonia gas distributions with this multi-filter method.Finally, a comparison of the goodness of fits between models 2 and 3 is presented in Figure 10g, showing that due to the low number of filters used to cover the Q-band, we cannot reliably identify any variability in para-H 2 (i.e., NEMESIS can adequately fit the 8-point spectra by assuming a para-H 2 fraction that remains static with time, see Figures 10h and 10i).Therefore, we neither retrieve ammonia gas nor the para-H 2 fraction in this study, focusing the remainder of this discussion on changes to aerosols and temperatures.Long-term spectroscopic observations, rather than filtered imaging, are likely to be needed to assess variability in NH 3 and para-H 2 .

Results: Retrieved Temperatures and Aerosol Opacity
Retrieved aerosol opacity at 600-800 mbar (constrained by the 8.6-μm filter) and temperatures in the stratosphere (i.e., 10 mbar, constrained by the 7.9-μm filter), upper troposphere (330 mbar, constrained by 18.7-24.5μm filters) and mid-troposphere (i.e., 500 mbar, constrained by the 10.77 and 13.04 μm filters) are shown in Figure 11.These are accompanied by the temporal variance of these temperatures and the corresponding Lomb-Scargle periodograms, showing potential cyclic behaviors of the temperature and aerosol opacity variations.These temperature and aerosol opacity maps are built by combining the 5-filter retrieval results between 1983 and 1996 and the 8-filter retrievals between 1996 and 2020, to obtain a map covering the full 36 years of our dataset.
Due to the different number of filters used in both retrievals, the retrieved temperatures and aerosol opacity using 5 or 8 filters differ in absolute values, although similar temporal tendencies are observed in both cases (see Figure S2 in Supporting Information S1).Differences in the absolute values between the 5-filter and the 8-filters retrievals are corrected in Figure S2 of Supporting Information S1 by shifting the results between 1983 and 1995 up or down depending on the latitude so that the absolute values of the 5-filter and the 8-filter retrievals are consistent in 1996 (when Jupiter was regularly observed with all 8 filters).The applied corrections fall within the retrieval uncertainties (see Figure S3 in Supporting Information S1 where the differences between the 8-filter and 5-filter retrievals are shown as a function of date and latitude).We therefore note that absolute values in Figure 11 are highly uncertain, due to the degeneracies inherent in imaging retrievals, and we only focus on the relative temporal and meridional variability of the temperatures and aerosol opacity, which can be considered to be robust.

Hemispherical Asymmetry
The northern hemisphere stratosphere at 10 mbar, poleward of 30°N, appears to be warmer than the corresponding southern latitudes, irrespective of the solar longitude (see Figure 12), in agreement with the warmer northern mid-latitudes reported in 2000 and 2014 by Fletcher et al. (2016), the preliminary results from brightness-temperature measurements (e.g., Antuñano et al., 2021;Orton et al., 2023) and the 1-D seasonal radiative seasonal model by Guerlet et al. (2020).At 40° latitude, for example, the northern region is on average 3.1 ± 2.6 K warmer than 40°S, while this temperature asymmetry weakens toward the equator, reaching differences of 2.7 ± 2.5 K on average between 30°N and 30°S.This apparently permanent asymmetry, with the warmer northern mid-latitudes, may correspond to the higher number density of aerosols found poleward 30°N at pressures smaller than 50 mbar compared to the southern latitudes, which could provide additional radiative heating (Zhang et al., 2013).This asymmetry could also arise from the combination of the former and the differences in hemispheric solar insolation (Guerlet et al., 2020) because the northern hemisphere receives 21% greater solar energy inflow due to Jupiter's northern summer solstice (solar longitude L s = 90°) occurring close to Jupiter's perihelion (L s = 57°).The northern hemisphere often exhibits enhanced stratospheric wave activity at northern mid-latitudes that may be connected to tropospheric disturbances (Fletcher, Orton, Sinclair, et al., 2017), which could provide additional mechanical forcing to drive this hemispheric asymmetry.
Figure 10.Maps of residual χ 2 values (a, d, g), showing differences in the quality of the spectral fits achieved by NEMESIS between model ii (i.e., retrieving stratospheric and tropospheric temperatures and tropospheric aerosol opacity simultaneously using the 5-point spectra) and model i (retrieving stratospheric and tropospheric temperatures and tropospheric ammonia gas and aerosol opacity simultaneously using the 5-point spectra) in panel a; model 2 (i.e., retrieving stratospheric and tropospheric temperatures and tropospheric aerosol opacity simultaneously using the 8-point spectra) and model 1 (retrieving stratospheric and tropospheric temperatures and tropospheric ammonia gas and aerosol opacity simultaneously using the 8-point spectra) in panel d; and comparing model 2 and model 3 (retrieving stratospheric and tropospheric temperatures, tropospheric aerosol opacity and para-H 2 abundances simultaneously) in panel (g).Actual χ 2 values of the fits using models i and ii (green and black dots, respectively) and models 1-3 (blue, black and red dots, respectively) for the equator (b, e, and h) and 10° north (c, f, and i), showing small variations between models in all panels except (f), where ammonia gas deviates substantially from the prior.The equator and 10° north are highlighted in (a), (d), and (g) by white dashed lines for reference.
Stratospheric temperatures at mid-latitudes are also observed to oscillate with a ∼12-year period in anti-phase between the northern and southern hemisphere, as reported by Simon-Miller et al. (2006).A Pearson correlation analysis of associate conjugate latitudes between 26° and 40° returns large negative correlations, with a mean correlation coefficient of −0.67, confirming the anticorrelated nature of the stratospheric temperature variations at temperate and mid-latitudes (Pearson coefficients as a function of latitude are displayed in Figure S4 in Supporting Information S1).
The top panel of Figure 12 shows a very complex variability of the temperatures with the solar longitude, suggesting that stratospheric temperatures are not solely modulated by radiative forcing.This is reinforced by the lag of around 1 year observed between the temperature maxima and the solar forcing between 1983 and 2008 (in agreement with Simon-Miller et al. (2006)), which is around half than the expected value from radiative response models (Li et al., 2018).Here we speculate that stratospheric temperatures at mid-latitudes are modulated by a combination of radiative and mechanical forcing from meteorological activity at the deeper levels (Fletcher et al., 2016;Guerlet et al., 2020;Simon-Miller et al., 2006).

NTB Disturbances
At lower latitudes -between 20° and 25°-similar stratospheric temperatures are observed in both hemispheres, with average differences of 1.3 ± 2.5 K and 0.2 ± 1.9 K at 24° and 20° latitude, respectively (see Figure 12).However, unlike in the southern region, a large temporal variance is found at the north temperate domain, suggesting the northern stratosphere to be more dynamic than the south.
At 20° and 25° latitude, stratospheric temperatures are observed to oscillate in time with a 5.5-6-year period (indicated by blue squares in Figure 11) in both hemispheres (in agreement with the brightness temperature maps in Figure 5a), displaying a maximum variation of ∼5-6 K at the north hemisphere and ∼3-4 K in the south (see Figures 11 and 12).This period is slightly longer than that of the NTB disturbances occurring deeper in the troposphere (i.e., ∼5-year period, Barrado-Izagirre et al., 2009;Rogers, 1992Rogers, , 1995;;Sánchez-Lavega et al., 2008, 2017), where the eruption of one or more strong convective plumes rising from the deep water cloud layer interact with the background flow forming wakes and vortices (Hueso et al., 2002), to finally completely disturb the ∼22°-27° latitude band at visible wavelengths (Sánchez-Lavega et al., 2017).These convective plumes have been previously observed to reach the lower stratosphere at ∼60 mbar (Sánchez-Lavega et al., 2008).However, no sign of the plumes have been observed at higher altitudes so far.
Figure 12c shows a stratospheric warming around a year after the NTB outbreaks marked in Figure 12c with stars.However, understanding whether the NTB disturbances could be affecting the stratospheric temperatures is challenging.NTB outbreaks could spawn waves that propagate to higher altitudes, warming the stratospheric temperatures at 10 mbar.This is in agreement with the large longitudinal variance (indicating the presence of waves) observed in 7.9-μm images at 20°-25°N in 1992 (Antuñano et al., 2021).However, we note that large temporal longitudinal variance is also found at epochs where no NTB outbreaks were observed (e.g., 1996 and 2002), and little variance is found after the 2008 NTB outbreak, although the latter could be due to the lack of global coverage in ground-based images (Antuñano et al., 2021).Monitoring and characterizing future NTB outbreaks at visible and mid-infrared wavelengths, as well as new numerical model studies will be crucial to understand whether stratospheric temperature variations are influenced by NTB outbreaks.

Jupiter's Equatorial Stratospheric Oscillation
Figure 11a shows the distinct signature of the JESO, with the equatorial and off-equatorial (i.e., 12°-14°) latitudes displaying the largest temporal variance of the stratospheric temperatures.The off-equatorial and equatorial latitudes present warm and cool temperature patterns in anti-phase.As in the 7.9-μm brightness temperature maps shown in Figure 5a, the JESO disruptions from 1992 to 2008 reported by Antuñano et al. (2021) are also observed in the retrieved stratospheric temperatures.The retrieved temperature profile at 13.5 mbar shown in Cosentino et al. (2020) and (Giles et al., 2020) from TEXES spectroscopic mapping (around the pressure level sensed by the 7.9-μm filter) suggests an early temperature increase at the equatorial latitudes in early to mid 2017, although our retrieved temperatures and brightness temperature profile do not show this increase (see Figures 5a  and 11a).However, the retrieved temperatures do show that the latest temperature maximum happened in early 2019, ∼0.5 year earlier than expected from a 4-year period, as the previous maximum happened in late 2015 (see Figure 11a).This is potentially related to the phase shift found at higher altitudes in Giles et al. ( 2020) and the 10.1029/2022JE007693 20 of 34 differences between these studies are potentially due to the 7.9-μm imaging filter blending together the flux from a range of altitudes, compared to the superior vertical resolution provided by the TEXES observations.The Lomb-Scargle periodogram shows a 4-year period (indicated by green squares in Figure 11) at the equatorial and off-equatorial latitudes, in agreement with Figure 5a and previous studies (Cosentino et al., 2017;Giles et al., 2020;Orton et al., 1991;Simon-Miller et al., 2006).Antuñano et al. (2021) showed that, although a ∼4-year period is found when fitting all the 7.9-μm data simultaneously, this does not adequately fit the observations prior to 1992 (where a 5.7-year period is found instead) and after 2008 (where the observations do not show any particular periodicity).The Lomb-Scargle periodogram does not hint at the longer period previously reported in the observations prior to 1992, obtained using a wavelet-transform analysis, but it can be clearly observed in Figures 5a and 11a, where an almost 6-year interval is observed between the brightness temperature maxima in 1984 and 1990.A wavelet transform analysis of the retrieved stratospheric temperatures similar to that performed by Antuñano et al. (2021) returns a statistically significant 4-year periodicity between ∼1994 and 1996, but does not show either the longer period reported in Antuñano et al. (2021) between 1980 and 1988 (see Figure S4 in Supporting Information S1) or the 4-year period between 2012 and 2019 (Giles et al., 2020).
Differences in the wavelet transform periodogram analysis between this study and Antuñano et al. (2021) come from the fact that this study does not include 7.9-μm observations from 1980 to 1982, unlike in Antuñano et al. ( 2021), as we could not create 5-point spectral cubes due to the lack of data in some of the filters used in this study during those years.This results in only one oscillation being identified between 1983 and 1990, instead of the 1.5 oscillations seen in Antuñano et al. (2021).Similarly, differences between the periodogram results in Giles et al. (2020) and ours are potentially due to the lack of 7.9-μm observations between early 2013 and early 2016, where our smoothing technique (see Section 2.2.3) might not be representing the real state of Jupiter's atmosphere.Retrieved stratospheric temperatures from 1984 to 1990 and 2016-2019 shown in Figure 11, however, do show the ∼6-year and 4-year time intervals, respectively, between temperature maxima.These confirm the change in the JESO's periodicity after the 1992 disruption, and reveal that the 2017 disruption of the JESO did not alter its periodicity.

Tropospheric Temperatures and Aerosol Opacity
Retrieved tropospheric temperature profiles as a function of time at 330 mbar (solid lines) and 500 mbar (dashed lines) are shown in Figure 13 for four different conjugate latitudes, representative of the temperature variations of the diverse regions described in this paper.The simultaneous representation of the temperatures at both tropospheric altitudes and conjugate latitudes enables us to perform a correlation analysis of changes, both in altitude and latitude at the same time (discussed in Section 6).
Jupiter's upper troposphere at 330 mbar displays temperature differences at mid-latitudes of ∼1 K at around 40°, although this falls within the formal uncertainties on our retrievals (see Figure 13a).This asymmetry is opposite and much weaker than that observed at the stratospheric mid-latitudes, and it is not observed at 500 mbar, where temperatures are largely symmetric at latitudes >20°, as shown in Figure 13b.At both altitudes tropospheric temperatures at latitudes >20° are less variable than those in the stratosphere.The off-equatorial tropospheric temperatures at 500 mbar exhibit the largest temporal variance, mostly due to the large variability observed in the NEB and SEB compared to other belts in Jupiter.
Tropospheric aerosol opacity in Jupiter's belts and zones, shown in Figure 11, is correlated with the albedo at visible wavelengths, and anti-correlated with the overall tropospheric temperatures.High aerosol opacity (see Figure 11) is found in the cold and typically white zones, whereas lower values are found in the usually warm and low-albedo belts, in agreement with volatiles condensing at cool temperatures and forming cloudy regions.The aerosol opacity is highest at the EZ, followed by the north temperate latitudes between ∼20°-40° and the south temperate regions at ∼20°-40°S.The NEB displays the lowest aerosol opacity and the highest tropospheric temperatures, with a difference in aerosol optical depth between the NEB and the SEB of −0.3 (at 10 μm), and maximum temperature differences between these two belts that could reach 1.8 K at 330 mbar and 5 K at 500 mbar (see Figures 11 and 13).The north and south temperate latitudes display approximately twice the aerosol opacity of the NEB.
Multiple periods are found in Figures 11b-11d, both in the temperature and aerosol opacity variations.The mid-latitude tropospheric temperatures poleward of 30° show 4-4.5-year periods (shown by red squares), while the north temperate latitudes (i.e., 20°-30°N) display a 7-year period for the temperature variations and 5-year cycle for changes in the aerosol opacity (indicated by orange squares).These periodicities differ from those found in the stratosphere (see Section 5.1).
Additionally, the ∼12-year periods that we originally found in the brightness temperatures at all wavelengths between 10.7 and 24.5 μm (see Figure 4) are not observed in the retrieved tropospheric temperatures at 500 mbar nor in the aerosol opacity, indicating that these were largely due to the change of the minimum emission angle over the 36 years (Orton et al., 2023), which are now properly accounted for in the temperature inversion.However, at 330 mbar retrieved temperatures continue to show a 10-14-year period for latitudes larger than 40° (yellow squares), in agreement with (Orton et al., 2023) at 20-30°.
The absence of a 5-5.5-year periodicity in the tropospheric temperatures at 20°-30°N suggests that although the NTB outbreaks darken the southern edge of the north temperate latitudes at visible wavelengths during the NTB disturbances (e.g., Rogers, 1995;Sánchez-Lavega et al., 2017), they do not affect the tropospheric temperatures, as expected from the mid-infrared images, where the NTB disruptions do not appear to have significant signatures.The lack of temperature variations during the NTB disturbances contrasts with the SEB revivals (Fletcher et al., 2011;Fletcher, Orton, Rogers, et al., 2017;Pérez-Hoyos et al., 2012;Rogers, 1992;Sanchez La Vega, 1989) and NEB expansions (described in the following sections), where the troposphere warms, accompanied by a removal of aerosols, allowing the 5-μm emission to escape from the deeper levels.This difference between the NTB disturbances and the SEB revivals, both triggered by convective storms, could come from various sources, such as the underlying environmental differences between the SEB and NTB, with the SEB notably depleted in NH 3 compared to the NTB, or differences in the accumulated convective available potential energy at these latitudes, which would result in differences between the magnitudes of the convective storms observed at the NTB and the SEB, where the former seem to reach higher altitudes (Sánchez-Lavega et al., 2017).The faster wind velocities at the NTB (Barrado-Izagirre et al., 2009) could also produce differences in the interaction of the storms with the background at the NTB and the SEB.

Thermal and Aerosol Opacity Variations During NEB Expansions
The periodograms of the aerosol opacity and tropospheric temperatures at 500 and 330 mbar peak at a 4-5-year period (see white squares in Figure 11), coinciding with the periodicity of the NEB expansions (shown as black lines in Figures 13c, 14d, and 14e and white dashed lines in Figures 14a-14c), where the low opacity of the NEB seems to expand poleward reaching the NTrZ(S) (e.g., Fletcher, Orton, Sinclair, et al., 2017;Rogers, 1995).This is in agreement with the tropospheric warming at 650 mbar reported during these events in previous studies (e.g., Figure 2c in Fletcher, Orton, Sinclair, et al. (2017)) and the brightness temperature maps at 10.7 and 13.0 μm shown in Figure 5.
The NTrZ (17°-22°N) displays no evidence of thermal changes related to NEB expansions at 330 mbar.However, temperatures at 500 mbar in Figures 13c and 14d, suggests increases of <1.5 K accompanied by a subtle decrease in the aerosol opacity (Figure 14c) during some (1988, 1994, 2009-2011, 2012, and 2016), but not all (1996 and 2004-2007) NEB expansions.Although there exists a degeneracy between tropospheric temperature and aerosol opacity due to the low number of filters used in this study, the reduction of the aerosol opacity during the NEB expansions is in agreement with the lower reflectivity at 890 nm shown in Fletcher, Orton, Sinclair, et al. (2017).
A correlation analysis between temperature and aerosol opacity variations at 16°N, suggests a temporal lag of 240 ± 60 days between these two, with aerosol opacity changes happening before the variation in the tropospheric temperature (see also Figures 14d and 14e).This suggests removal of the white aerosols at the NTrZ to reveal the darker chromophore at deeper levels during the NEB expansions.However, we note that this study cannot constrain the altitude of the chromophore, which is still an open question (e.g., Braude et al., 2020;Dahl et al., 2021;Pérez-Hoyos et al., 2020;Sromovsky et al., 2017).

Thermal and Aerosol Opacity Variations During SEB Fading and Revival Events
Since 1983, the SEB has undergone three fade and revival events (1989-1990, 1992-1993, and 2009-2011) and a partial fading (in 2007), where the SEB transforms from being the darkest and broadest belt on Jupiter to a faded and white zone, and restores its normal dark brown coloration after several years (e.g., Rogers, 1995).These fading and revival events are represented in Figure 15 by horizontal black lines and vertical dashed white lines.
The retrieved aerosol opacity at 400-600 mbar shows significant variability related to the aforementioned SEB fading and revival events at all latitudes between 8° and 18°S.The aerosol opacity displays large increases during the SEB fading of 1989-1990 and 2009-2010, with optical depth increases of ∼1 on average at 12°-18°S, and around 0.85 on average at 8°-10°S.During these two whitening events our aerosol opacity retrievals show that, although the aerosol opacity seems to start increasing almost contemporaneously at all latitudes, the most equatorward latitudes of the SEB reach their peak in aerosol opacity months before the poleward latitudes of the SEB.This is not in complete agreement with the observed evolution of the SEB fading at visible wavelengths, where the whitening of the SEB seems to start at mid-SEB latitudes and rapidly expand into lower latitudes (e.g., Rogers, 2017).We note that our study deals with zonally averaged profiles and that at some dates we do not have complete global coverage, leading to a potential loss of the start of the SEB fading.
In 2007, the SEB underwent a partial fading in the visible and at 5-μm wavelength (the later sensing the radiance below the ammonia cloud) that lasted around half a year and where the most equatorward latitudes of the SEB remained undisturbed (e.g., Rogers, 2007).Our retrievals show that the aerosol opacity at 12°-18°S started to increase in early 2007, reaching its maxima in early-to mid-2007 and displaying its lowest aerosol opacity value in early 2009, right before it started to increase again for the 2009-2010 fading event.At the northern latitudes of the SEB between 8° and 10°S, the aerosol opacity seems to stay almost invariant between 2006 and mid-2007 and then decreases following the tendency seen at 12°-18°S.
Considering the anticorrelation between aerosol opacity and tropospheric temperature changes, one would also expect to observe significant changes in the tropospheric temperatures during the SEB fading and revival events.In fact, Fletcher, Orton, Rogers, et al. (2017) reported a 2-4 K localized warming of the troposphere at 500 mbar during the revival of the faded SEB in 2010-2011, although no significant changes were reported during the fading episode (Fletcher et al., 2011).Our zonally averaged tropospheric temperatures at 330 and 500 mbar show a very complex temporal variability, with some epochs (2009 to mid-2010) hinting at a 0.5-1 cooling of the SEB during the fading, while others showing the complete opposite, a cooling of the temperatures during SEB revivals.These perplexing results demonstrate that further analysis using spectroscopy, fully accounting for longitudinal variability, is needed to be able to disentangle temperature and aerosol changes during SEB fading and revival cycles.

The Equatorial Zone
At the equatorial latitudes (i.e., ±7° of the equator), tropospheric temperatures vary with time, with maximum temperature contrasts of ∼1.5 K, both at 330 and 500 mbar (see Figure 13d).The Lomb-Scargle periodograms in Figure 11 display both 4 and 8-year periods for the oscillations of the upper-tropospheric temperatures at 330 mbar (marked by green squares), while a 8-to 9-year periodicity is found for the temperature variations at 500 mbar.The periodicities at 330 mbar are in agreement with the brightness temperature changes at 17.6-24.5μm (Figure 4) and (Orton et al., 2023).As mentioned in Section 3, the 4-year periodicity observed in the upper-tropospheric temperatures at 330 mbar coincides with the period of JESO, suggesting a potential coupling between Jupiter's equatorial upper troposphere and stratosphere.This coupling in further discussed in Section 6.
Brightness temperatures at 10.77-μm and 13.0-μm shown in Section 3, revealed a ∼7-year period at the equatorial latitudes, coinciding with the EZ disturbances (see Figure 4).This period, however, is not observed in the retrieved temperatures at 500 mbar, indicating that the equatorial troposphere might not sufficiently warm during the remarkable EZ disturbances to show up as clearly in our retrievals.Although hard to confirm due to the lack of data during some of the EZ disturbances (e.g., 1992 or after mid-2019), the 500-mbar tropospheric temperatures seem to increase before the EZ disturbances appear at 5 μm, mainly in 2006 and 2018 (see Figure 13d), and decreases again right after the EZ disturbances finish at 5 μm.However, this temperature pattern is not observed in all the EZ disturbances known so far; we note that a further analysis using spectroscopy rather than photometry will surely improve our understanding of the temperature variations during these rare disturbances.Nevertheless, if temperatures do increase before the EZ disturbances, this increase would be smaller than 1.5 K. Whether such a variation could clear the ammonia clouds by sublimation is still unclear.
At these latitudes, the aerosol opacity in Figure 16 is observed to start decreasing around 1-1.5 years before the brightening of the EZ at 5 μm.It usually takes around 4 years to restore the typical equatorial aerosol opacity.This is in agreement with the observed duration of the coloration events at visible wavelengths (Antuñano et al., 2018;Rogers, 1995) and the 8.6-μm brightness temperature changes in Figure 5b, confirming that EZ disturbances start at the ammonia cloud level and expand downward with time to finally brighten the EZ at 5 μm.The contemporaneous decrease of the albedo and aerosol opacity also agrees with the hypothesis of the EZ disturbances being the result of a decrease in the ammonia upwelling due to a potential change in the deep ammonia column beneath the EZ (Bolton et al., 2017;Ingersoll et al., 2017;Li et al., 2017), or due to strong precipitation surrounding the large plumes observed at the northern edge of the NEB (Antuñano et al., 2020).Our study cannot distinguish between these two cases.A correlation analysis between aerosol opacity and tropospheric temperature changes returns the maximum correlation factor (0.41) for lag of 180 ± 60 days between tropospheric temperature variations and aerosol opacity changes.This means that a moderate correlation is found between these two magnitudes, with temperature variations leading to the changes we observe in the aerosol opacity.This clarifies that the 8.6-μm changes that we observe are first due to temperature changes at around 500-600 mbar and later due to aerosol opacity changes.Finally, the aerosol opacity decrease observed in 2012-2013 did not lead to 5-μm brightening of the equatorial latitudes nor to strong coloration events as did the decreases in 1992, 2000, and 2007.Figure 16 also shows a secondary decrease of the aerosol opacity, usually observed halfway between EZ disturbances.This secondary decrease of the aerosol opacity at the equator is more subtle than during the EZ disturbances, with decreases of around half of those found during the EZ disturbances.These shorter-period (3.5-4 years, shown with a light blue square in Figure 11) of aerosol opacity changes at the EZ are somewhat surprising, in that they have no counterpart in the tropospheric temperatures at 500 mbar, nor are they related to notable coloration changes or tropospheric wind speed changes (Tollefson et al., 2017;Wong et al., 2020).This suggests that the aerosol opacity at the equator is modulated not only by the EZ disturbances, but also by other dynamical forces.
These shorter-period changes in the equatorial aerosol opacity do not follow the anticorrelation observed in the tropospheric temperatures at the NEB and SEB, which present a longer period (see Section 6.3), confirming that the potential connection between the NEB and SEB originates deeper in the atmosphere and does not appear to influence the equatorial latitudes.Understanding how ammonia varies with time, both at the ammonia cloud tops sensed via mid-infrared spectroscopy and the deeper levels sensed by Juno's MWR, will be essential to investigate the nature of the subtle aerosol opacity decreases found between the EZ disturbances.
Finally, an asymmetry in the aerosol opacity is also observed at the EZ, with higher aerosol opacity found at 0°-5°S (EZ(S)) compared to 0°-5°N (EZ(N)) overall (see Figure 11d).This asymmetry is present at visible wavelengths too, where the EZ(S) sometimes displays higher albedo than the EZ(N) (e.g., Rogers, 1995).However, the asymmetry, both at the visible wavelengths and aerosol opacity, is opposite to the observed NH 3 distribution Bolton et al. (2017), and could be an effect of the richer cloud morphology observed at the EZ(N) related to the NEB hot spots and plumes (Fletcher, Orton, et al., 2020;Orton et al., 1998;Westphal, 1969) compared to the more quiescent and calm EZ(S).

Discussion
The previous sections have described significant variability in Jupiter's stratospheric and tropospheric temperatures and aerosol opacity, with strong correlation/anticorrelation to previously known cyclic activity and hemispherical changes.A summary of the observed variability is given in Table 2.In this section, we further discuss the most remarkable temperature and aerosol opacity variations, focusing on understanding their nature by characterizing the potential connections between tropospheric and stratospheric changes.

Stratospheric Temperature Asymmetry at Mid-Latitudes
The consistent calibration and analysis approaches to the entire mid-infrared imaging time series have revealed that (a) the northern mid-latitude stratospheric temperatures at 10-20 mbar pressure level are overall 2.7-3.1 K warmer on average than the southern mid-latitude stratospheric temperatures; and (b) temperature variations are highly anticorrelated between hemispheres, with periodicities close to a Jovian year.The consistency of higher brightness temperatures at 7.9 μm and retrieved stratospheric temperatures at the northern mid-latitudes independent of the solar longitude, provides compelling evidence of a hemispherical thermal asymmetry that is weakly modulated by the jovian seasons.This is in agreement with the warmer northern mid-to high-latitudes predicted from the 1-D radiative-convective models of Guerlet et al. (2020), where they attribute the north-south asymmetry to the strong asymmetry in the polar haze abundance as constrained by Zhang et al. (2013).The hemispherical asymmetry in the stratospheric temperatures, aerosol number density, C 2 H 2 and C 2 H 6 could be related to the auroral oval, which is larger and extends further south in the northern hemisphere compared to the more compact auroral oval in the south.However, to date, there has been no systematic study of the time variability of polar stratospheric aerosols to confirm this.
Predicted temperature variations from the 1-D radiative-convective model of Guerlet et al. (2020) are shown in Figure 17 as dashed lines for 40°N and 40°S.The southern mid-latitude stratospheric temperatures do not follow the radiative climate predictions of Guerlet et al. (2020), which expected a smaller amplitude of temperature variations in the south, and that the warmest temperatures would be encountered during mid-northern (rather than southern) summer, because Jupiter is closer to perihelion during southern winter.The temperature time series in our southern dataset in Figure 17 is apparently more random than in the north.Nevertheless, the data  reveal the opposite trend to the model (warmer southern temperatures are experienced during southern spring and summer), suggesting that Jupiter's orbital eccentricity plays a smaller role than that accounted for in the model, in agreement with Orton et al. (2023).Furthermore, our data generally have temperature contrasts twice as large as predicted by the model, indicating that radiative heating of the southern hazes may continue to dominate the temperature trends.Future studies of stratospheric aerosol variability with latitude and time would greatly improve the ability of radiative climate models to predict Jupiter's stratospheric temperatures.
Finally, the strong anticorrelation of the temperature variations between hemispheres is still a mystery.Quasi-seasonal variations are observed also at other latitudes, such as JESO and at other planets like Saturn and Earth.However, these are usually observed at low latitudes.Additionally, the anticorrelation found at changes in the NEB and SEB could be potentially explained by deep connections via cylinders parallel to the rotation axis (see Section 6.3).However, these are most effective at latitudes lower than 16°.Further analysis and models are needed to understand the mechanisms responsible of such anticorrelation.

Stratosphere-Troposphere Coupling
Many studies from the past decades have shown that, in the Gas Giants, deep tropospheric meteorological and/or convective activity can affect the temperature and aerosols in the upper troposphere at ∼100-330 mbar and lower stratosphere at 5-60 mbar (e.g., Li & Ingersoll, 2015;Sugiyama et al., 2014).In most cases, this is thought to be related to wave activity that carries internal heat flux from the deeper levels into the stratosphere, either due to deep convective activity or generated by wind-shear instabilities at the cloud level.
The power spectrum analysis of the retrieved stratospheric and upper-tropospheric temperatures at 10 and 330 mbar, respectively, shows a 4-year period for changes in their temperature profiles at the equator.As mentioned in previous sections, this hints at a potential connection between temperature changes at these two pressure levels.However, the power spectrum alone cannot determine the kind of connection between temperature changes at 10 and 330 mbar.For this reason, in Figure 18 we compare the temporal variability of the retrieved temperatures at 10 and 330 mbar, which enables us to investigate whether such a connection is really seen in the temperature profiles.
Figure 18 shows that the retrieved temperatures at 330 mbar are observed to be anticorrelated to those at 10 mbar, even during the epochs when the JESO was disrupted potentially due to energetic "Global Upheaval" events (Antuñano et al., 2021).A Pearson correlation analysis returns a low-to-moderate anticorrelation with a coefficient of −0.47 at the equator.The low correlation found in the Pearson correlation analysis is likely related to the lack of observations at some epochs.During these epochs the retrieved thermal profiles might not represent the real state of the atmosphere and therefore, one needs to be cautious when analyzing the results.The anticorrelation observed in Figure 18 is likely to be real as it also matches deviations from the usual state of JESO during the disturbed epochs, suggesting a potential tropospheric-stratospheric coupling.This is in agreement with the anticorrelation observed between brightness temperature changes at 7.9 μm and temperatures at 330 mbar between 1981 and 2011 (Orton et al., 2023).This anticorrelation between temperature temporal variability at 10 and 330 mbar does not extend deeper into the troposphere, as retrieved temperature and aerosol opacity at the 500 mbar level do not show any correlation between their variability and JESO.
A study of Jupiter's stratospheric temperatures between 2012 and late 2019 using high-vertical resolution spectra from TEXES (Giles et al., 2020) shows the descending pattern of warm and cold temperature anomalies and their evolution with time (their Figure 3).Using that figure, and assuming a mean scale-height of 23 km at around 10 mbar (Leovy et al., 1991), we obtain a descending rate of ∼13.6 km/year.A ∼13.6 km/year descending-rate of the JESO means that the temperatures at 10 mbar take ∼5.25 years to reach the 330 mbar level.A Pearson correlation analysis between stratospheric and upper tropospheric temperatures shows a correlation maxima of 0.4 for a lag of 2.5 years.This is roughly half of the descending time of a particular temperature anomaly, and if genuine, it could indicate that more than one temperature anomaly (positive or negative) are present between 10 and 330 mbar.Therefore, we believe that our results reveal for the first time that Jupiter's upper-tropospheric temperature variations are influenced by the equatorial stratospheric temperature oscillation, with the vertically stacked chain of warm/cool airmasses propagating downward into the upper troposphere.By changing the equatorial temperature contrasts, this could also influence the shear on the equatorial jets near the tropopause via the thermal wind equation (although this has not yet been observed due to the challenges of measuring winds at these pressure levels) and the aerosol condensation.This, however, has previously been observed in Earth's Quasi-Biennial Oscillation, while Saturn's Equatorial Stratospheric Oscillation has been observed to almost reach the tropopause (Schinder et al., 2011).

NEB-SEB Anticorrelation
The aerosol and thermal variability of the SEB and NEB have been characterized by numerous studies, focusing on observations during SEB fade/revivals and NEB expansion events.However, the lack of a consistent long-term analysis approach made it difficult to investigate the long-term behavior of the clouds, aerosols and temperatures.2019) reported a potential anti-correlation between the 5-μm brightness temperature changes at these two belts, suggesting a potential coupling between belts that are separated by ∼15° latitude.
Although difficult to detect in Figure 11, our retrieved aerosol opacity and tropospheric temperatures at the NEB and SEB, also seem to vary in antiphase with one another.This is more clearly observed in Figure 13c, where the tropospheric temperatures at 500 mbar (dashed line) and 330 mbar (solid line) are shown for 16° latitude (north and south), and Figure 19, where the tropospheric temperature and aerosol opacity profiles are represented between 10° and 16° latitude.
A Pearson correlation analysis of associate conjugate latitudes at the SEB and NEB (shown in Figure S5 in Supporting Information S1) returns a moderate negative correlation of the aerosol opacity and tropospheric temperatures, with maximum correlation coefficients of −0.53 at 14° latitude in the aerosol opacity, −0.52 at 12° latitude in the tropospheric temperatures at 500 mbar and between −0.6 and −0.67 at 14° and 16° latitude in the upper-tropospheric temperatures at 330 mbar, hinting in a potential connection between these two belts as found at 5 μm.
The continued long-term connection between belts separated by ∼15° latitude is still not understood.Rogers (1995) describes epochs where a number of belts and zones in Jupiter (mainly the SEB, EZ, NEB, and NTB) seem to display a disturbed appearance at visible wavelengths almost contemporaneously.These are called "Global Upheavals."During these events, Jupiter's equatorial and temperate latitudes display major disturbances as if one of the perturbations would trigger the rest of the disturbances, hinting in some kind of deep connections between these latitudes.Although one could think that the NEB-SEB anticorrelation are related to the Global Upheavals, these seem to occur only at specific epochs and are thought to be more frequent during SEB revival events Rogers (2007), unlike the long-term coupling of the NEB and SEB.
A potential cause for an anticorrelation in the tropospheric temperatures and aerosols opacity between belts situated in the same hemisphere, could be the meridional propagation of tropospheric waves, following a similar mechanism to that observed in Saturn during 2011-2014 (Fletcher, Guerlet, et al., 2017).At that epoch, stratospheric waves emanating from a hot and large stratospheric vortex near 40°N disrupted Saturn's Equatorial Stratospheric Oscillation (Fletcher, Guerlet, et al., 2017;Guerlet et al., 2018).However, meridional propagation of waves between hemispheres is unlikely due to the change in the Coriolis force.
Here we speculate that the NEB-SEB anticorrelation observed between 2 and 5 bar and 330 mbar originates in the deep troposphere, via a deep connection of these belts via cylinders parallel to the rotation axis.These would be most effective equatorward of ±16° where the cylinders do not intersect with the dynamics of a region of metallic hydrogen (Cao & Stevenson, 2017;Liu et al., 2008), in agreement with the highest negative correlation coefficients found at ∼14° latitude.Future studies analyzing the temporal variability of Jupiter's deep atmosphere using Juno Microwave Radiometer could shed some light on the causes of the NEB and SEB anticorrelation.

Conclusions
The continued monitoring of Jupiter in the mid-infrared from ground-based observatories during the last four decades provides essential long-term context of Jupiter's atmospheric variability.In this study, we made use of a large dataset of images captured between 7 and 25 μm between early 1984 and late 2019 to explore Jupiter's climate over three jovian years.In particular, we (a) characterize the long-term variability of Jupiter's atmosphere in each of the wavelengths analyzed; (b) retrieve stratospheric and tropospheric temperatures, as well as tropospheric aerosol opacity, to explore hemispherical asymmetries, both seasonal and non-seasonal changes, and equatorial oscillations; and (c) investigate thermal and aerosol opacity changes during cyclic and non-cyclic disturbances of the belt/zone structure.Although some of the mid-infrared images used in this study have been previously published (e.g., Antuñano et al., 2021;Fletcher, Orton, et al., 2018;Fletcher, Orton, Sinclair, et al., 2017, Orton et al., 1994, 2023), this is the first study that analyses simultaneously long-term ground-based images captured in 5 and 8 filters in a systematic fashion.The conclusions of this study are summarized below: • Hemispherical asymmetry: Stratospheric temperatures at 10 mbar retrieved from the 5-point and 8-point spectral cubes show persistent warmer northern latitudes poleward of 30°, compared to their conjugate latitudes in the south.The stratospheric temperature asymmetry between hemispheres increases with latitude, reaching differences of 3.1 ± 2.6 K at 40° and 2.7 ± 2.5 K at 30°.The nature of this asymmetry is still unknown, although it could be related to (a) the higher number density of aerosols found poleward of 30°N at 50 mbar (Zhang et al., 2013) leading to stronger radiative heating, and (b) the differences in hemispheric solar insolation (Guerlet et al., 2020).In contrast, the southern upper troposphere at 330 mbar appears to be overall warmer than the northern latitudes poleward of 30°.However, thermal differences between hemispheres at this pressure are only around 1 K.At deeper levels this asymmetry is no longer detected.• The temporal variability of stratospheric and upper tropospheric temperatures also display a remarkable asymmetry between hemispheres at latitudes poleward of 30°.At these pressures and latitudes, temperatures are observed to oscillate with a 10 to 14-year period (i.e., very close to the jovian year) in anti-phase between hemispheres, as reported by Orton et al. (2023).Comparisons of these thermal oscillations to seasonal variations reveal that they are unlikely to be solely related to radiative heating as temperature peaks do not coincide with solstices nor with the expected seasonal lag (Orton et al., 2023).Numerical radiative-climate models (Guerlet et al., 2020) predict smaller peak-to-peak thermal oscillations in the southern hemisphere, suggesting that the thermal oscillation that we observe is not strictly seasonal.• Belt/zone structure: Our results confirm the anticorrelation between the aerosol opacity and tropospheric temperatures at 500 mbar -typically white and cold zones present high aerosol opacity, while usually warm and low-albedo belts display low aerosol opacity.The NEB displays the lowest aerosol opacity and the largest tropospheric temperatures, while the contrary is true for the EZ.
• In this study we also confirm thermal and/or aerosol opacity changes during NEB expansion and SEB fading events.During some, but not all, NEB expansions, the NTrZ at 17°-22°N displays subtle increases of <1.5 K at 500 mbar accompanied by small decreases in the aerosol opacity, in agreement with Fletcher, Guerlet, et al. (2017).During SEB fading events, our retrievals show that, although the aerosol opacity seems to start increasing almost contemporaneously at all latitudes in the SEB, the most equatorward latitudes (8°-10°S) reach their maxima in aerosol opacity months before the poleward latitudes of the SEB (12°-18°S).Retrieved tropospheric temperatures at 330 and 500 mbar do not show a clear variability during SEB fading and revival events.• At the EZ, this study confirms the variability of the aerosol opacity observed during the EZ disturbances in Antuñano et al. (2018Antuñano et al. ( , 2020)).However, unlike in previous studies, tropospheric temperatures at 500 mbar appear to increase by less than 1.5 K immediately prior to disturbances at 5 μm, although this is hard to confirm due to the lack of data during some of the EZ disturbances.Changes in aerosol opacity appear to lag 180 ± 60 days behind tropospheric temperature changes.Finally, the aerosol opacity at the EZ shows a secondary decrease halfway between disturbances, with a periodicity of 3.5-4 years.This unexpected secondary decrease does not follow the coloration events nor changes in the NEB and SEB.• NEB-SEB anticorrelation: Antuñano et al. ( 2021) reported a continued anticorrelation in the variability of the 5-μm radiance at the NEB and SEB.In this study, we confirm that this long-term anticorrelation is also present higher up in the atmosphere, with maximum Pearson coefficients of −0.53 in the aerosol opacity between conjugate latitudes at 14° at 500 mbar, −0.52 in the tropospheric temperatures at 12° at 500 mbar and −0.6 in the upper-tropospheric temperatures at 14°.This long-term coupling between the NEB and SEB spans all the dates analyzed in this study  but, perplexingly, not all of the temperature/aerosol variations are directly related to the significant meteorological events (fades, revivals, expansions, etc.) characterized in visible light.Put another way, the deep-seated anticorrelated connection between the NEB and SEB may be a triggering mechanism for global upheavals, but they do not always lead to upheavals of the same magnitude.• Stratosphere-troposphere coupling: Retrieved stratospheric and upper tropospheric temperatures reveal, for the first time, that JESO extends down to the 330-mbar pressure level.In particular, retrieved temperatures at 330 mbar are overall observed to be moderately anticorrelated to those at 10 mbar, with a Pearson correlation coefficient of −0.47.The Pearson correlation analysis also returns a correlation maxima of 0.4 for a lag of 2.5 years between stratospheric and upper-tropospheric temperature maxima.This lag is half of the expected descending rate of ∼13.6 km/year derived from Giles et al. (2020), which would indicate a lag between 10 and 330 mbar of ∼5.25 years.
These long-term records of thermal variability help to place snapshots from visiting spacecraft (Galileo, Cassini, Juno, etc.) into a broader context, revealing the uniqueness of the jovian climate at a particular moment.Further ground-based spectroscopic data sets will allow us to better disentangle temperature, ammonia and aerosol changes.By extending this time series into the era of the James Webb Space Telescope and the forthcoming JUICE and Europa Clipper missions, new atmospheric discoveries in the 2020s and 2030s can be understood in terms of the meteorological variability observed since the 1980s.

Figure 1 .
Figure 1.Temporal coverage of observations as a function of wavelength.Colors represent different instruments.Note that the wavelength of the observations might not exactly match the wavelength annotated here.See Section 2.2.1 for information on the wavelength shifting performed in this study.A summary of these data is given in Table1and a full description is available in TableS1of Supporting Information S1.

Figure 2 .
Figure 2. Four examples of Jupiter images captured 8.6 μm, showing the evolution of ground-based mid-infrared imaging over the past 4 decades (a).From left to right, images in (a) were captured by IRTF/BOLO-1 (1989), IRTF/MIRAC (1999), and Very Large Telescope (VLT)/VISIR (2007, 2018).Cylindrical maps of Jupiter at 7.9-24.5 μm captured by VISIR and Cooled Mid Infrared Camera and Spectrometer instrument on 25-27 May 2018 (b-i).At 20.5 μm (i) images from April 1 and 21 August 2018 were also used.VISIR images in this Figure are reconstructed to remove the obscured regions resulting from the small field of view of VISIR and the chopping amplitude limitation of 25″ of the VLT.This correction leads to unrecoverable lost signal in thin regions seen as black arcs or strips in the maps (b-i).

Figure 3 .
Figure3.Zonal-mean brightness temperature maps of Jupiter as a function of time between ±48° latitude at eight different wavelengths, showing the atmospheric variability at each of the wavelengths analyzed in this study.Note that the spectral coverage increases in the mid-1990s with the introduction of the MIRAC instrument.Gray shadowed regions represent epochs where no data are available for 2 years or longer.In these cases, we use the value of the brightness temperature of the last available epoch and assume it remains constant during the period without any available data.We note that this is only for the representation and in our analysis we use the averaged and smoothed profiles shown in Figure4.

Figure 4 .
Figure 4. Examples of the average smoothed 8.6-μm radiance profiles smoothed over 3.95 years (red solid line) at the equator (top) and 16° south (bottom), compared to the 8.6-μm zonal-mean radiance (black dots).Black error bars represent measurement errors computed by adding calibration uncertainties and the standard deviation of the zonal mean profiles in quadrature.Pink shadowed regions represent the 1σ uncertainty of the averaged smoothed radiance profiles, as described in the main text.

Figure 5 .
Figure 5. Brightness temperature anomaly maps of Jupiter in the N-band as a function of time between ±48° (a-d, left).Middle and right panels represent the temporal variance and the Lomb-Scargle periodogram of the brightness temperatures, respectively.Note that only periods with 98% of significance or higher are shown.Brighter periods in the right panel correspond to higher spectral powers.Different boxes represent the periods analyzed in this study.Residual brightness temperatures at each wavelength are computed by subtracting the average brightness temperature over all latitudes and dates from the smoothed zonal-mean brightness temperature profile (see Section 2.2.3).Gray shadowed regions represent epochs where no data are available for 2 years or longer.

Figure 6 .
Figure 6.Same as Figure 5, but for the Q-band wavelengths.

Figure 7 .
Figure 7. Smoothed brightness temperature profiles at 8.6 μm (pink), 10.7 μm (gray), and 13.0 μm (yellow) for the equator at 0° latitude, averaged over 30° longitude of the minimum emission angle.Pink, yellow and gray shadowed regions represent the 1σ uncertainty.Blue shadowed regions represent gaps in our dataset.Dashed black lines represent the 13.0-μm brightness temperature minima.10.7 and 13.0 μm brightness temperature profiles are offset by 10 and 12 K, respectively, for clarity.Note the increase in the 8.6-μm brightness temperature during the Equatorial Zone disturbances in 2000, 2007 and late 2018, indicated by black arrows.The increase in 2012-2013 also indicated by a black arrow did not lead to a disturbance at 5 μm.

Figure 8 .
Figure 8.Comparison of the minimum emission angle at 40°N as a function of date (dashed line) and the 20.5-μm brightness temperature variations at the same latitude (solid line), showing a clear anticorrelation between them.The shadowed gray region represent the measurement errors described in Section 2.2.

Figure 11 .
Figure 11.Retrieved temperatures, their variance, and the Lomb-Scargle periodograms at (a) 10 mbar, (b) 330 mbar, and (c) 500 mbar and aerosol opacity at (d) 400-600 mbar as a function of time and latitude.The middle panel in (a-d) shows the temporal variance of the temperatures and aerosol opacity, while right panel in (a-d) represents the Lomb-Scargle periodogram, showing Jupiter's atmospheric variability between June 1983 and November 2019 and potential cyclic changes.Note that only periods with 98% of significance or higher are shown.Brighter periods in the right panel correspond to higher power spectrum.

Figure 12 .
Figure 12.Retrieved stratospheric temperature profiles at 40° (top), 30° (middle), and 24° (bottom) for the north (dashed lines) and south (solid lines) as a function of time.Light and dark gray shadowed regions represent the formal uncertainties of the retrieved temperatures.Black dots and stars in the bottom panel indicate epochs of North Equatorial Belt expansions and North Temperate Belt outbreaks, respectively.

Figure 13 .
Figure 13.Retrieved tropospheric temperature time series at 330 mbar (solid lines) and 500 mbar (dashed lines) for latitudes (a) 40°, (b) 24°, (c) 16°, and (d) 2° in the northern (red) and southern (green) hemisphere as a function of time.Green and red shadowed regions represent the formal uncertainties of the retrieved temperatures.Black dots in (b) indicate the eruption of a convective storm (start of the North Temperate Belt disturbances).Black and red horizontal lines in (c) represent the epochs of North Equatorial Belt expansions and South Equatorial Belt fadings, respectively.The blue horizontal lines in (d) represent the epochs of Equatorial Zone (EZ) disturbances at 5 μm (Antuñano et al., 2018), while the red horizontal line in the same panel indicates a reddening of the EZ without a 5-μm brightening.

Figure 14 .
Figure 14.Retrieved upper-tropospheric temperature contours at (a) 330 mbar, (b) 500 mbar, and (c) aerosol opacity at 400-600 mbar showing the North Equatorial Belt (NEB) and the NTrZ.Retrieved tropospheric temperature at 500 mbar and aerosol profiles (d and e), respectively, are also shown for clarification.White dashed lines in (a-c) and the horizontal black lines in (d and e) represent the beginning and end of the NEB expansions described in Fletcher, Orton, Sinclair, et al. (2017).

Figure 15 .
Figure 15.Retrieved upper-tropospheric temperature contours at (a) 330 mbar, (b) 500 mbar, and (c) aerosol opacity at 400-600 mbar showing the South Equatorial Belt (SEB).Retrieved tropospheric temperature at 500 mbar and aerosol profiles (d and e), respectively, are also shown for clarification.White dashed lines in (a-c) and the horizontal black lines in (d and e) represent the beginning and end of the SEB fading described in Antuñano et al. (2019).

Figure 16 .
Figure 16.Retrieved aerosol opacity profiles at 400-600 mbar for the Equatorial Zone (EZ), showing the aerosol opacity decreases during the EZ disturbances (represented by the dashed red squares) and a secondary peak found around 3.5-4-years after the EZ disturbances.The aerosol opacity decreases in 2012-2013 did not lead to severe coloration events and 5-μm brightening.

Figure 17 .
Figure 17.Retrieved temperature profiles at 10 mbar (solid lines) for 40°N (top) and 40°S (bottom) compared to predicted temperatures in the "all aerosols" model (dashed lines) in Guerlet et al. (2020) and the seasonal insolation cycle (dotted lines).To match the retrieved temperatures, model temperatures in the northern and southern hemispheres have been decreased by 1.0 K and increased by 0.3 K, respectively.

Figure 18 .
Figure 18.Retrieved temperature profiles at 10 mbar (top) and 330 mbar (bottom) for the equator, showing a clear anticorrelation between changes at these pressures.Dark and light gray shadowed regions indicate the uncertainty in the retrieved temperatures.Red shadowed regions indicate epochs without data.

Figure 19 .
Figure 19.Retrieved aerosol opacity (dashed-dotted line), temperatures at 500 mbar (dashed line) and at 330 mbar (solid line) profiles for the North Equatorial Belt and South Equatorial Belt, showing a clear anticorrelation in the changes at these two belts, mainly at 12°-16° latitude.

Table 1
Summary of the Instruments, Date of Observation, and Wavelength Range and Plate Scale of Each Instrument

Table 2
Summary of the Periodicities Described in Section 5