Strong Resemblance between Surface and Deep Zonal Winds inside Jupiter Revealed by High-degree Gravity Moments

Jupiter's visual appearance is dominated by belts, zones, and vortices of varying sizes. Shear in the east–west zonal flows and the associated upwellings and downwellings are closely linked to the belts and zones, although the origin of the alternating bands of zonal winds is still a topic of debate. Moist convection in the shallow atmosphere (Vasavada & Showman 2005; Schneider & Liu 2009; Liu & Schneider 2010), deep convection in the molecular envelope (Busse 1976; Heimpel et al. 2005; Gastine & Wicht 2021), and tides from the moons of Jupiter (Lindzen 1991; Tyler 2022) have been proposed as possible underlying mechanisms. Numerical models for each of these proposals have been constructed, with varying degrees of success in reproducing the observed zonal wind patterns and amplitude.

Determination of the depth of the rapid atmospheric zonal flows provides a critical observational clue regarding their origin and dynamical balance, even though this depth would not directly translate to the driving mechanism (e.g., Showman et al. 2006; Christensen et al. 2020). Measurements of the gravitational (Hubbard 1999; Iess et al. 2018; Kaspi et al. 2018) and magnetic field (Cao & Stevenson 2017a; Moore et al. 2019; Cao et al. 2020; Bloxham et al. 2022) can be used to infer structure and dynamics inside giant planets. Gravitational harmonics are increasingly more sensitive to the outer layers as the spherical harmonic (SH) degree increases (Guillot et al. 2004). The magnetic field is sensitive to depth with sufficient electrical conductivity (Liu et al. 2008; Cao & Stevenson 2017a), e.g., greater than ∼2000 km inside Jupiter (French et al. 2012).

Previous gravity solutions from the Juno radio tracking experiments (Folkner et al. 2017; Iess et al. 2018; Durante et al. 2020) have provided the zonal gravitational harmonics of Jupiter up to SH degree ∼10. These gravity solutions were employed to constrain the background structure (Debras & Chabrier 2019; Miguel et al. 2022; Militzer et al. 2022) as well as dynamics, in particular the deep zonal wind structure, inside Jupiter (Kaspi et al. 2018; Kong et al. 2018; Kulowski et al. 2021). One continuing debate regarding the deep zonal winds inside Jupiter is to what extent they resemble the observed surface winds, as deep winds that are either quite similar to (Kaspi et al. 2018; Kong et al. 2018; Kulowski et al. 2021; Militzer et al. 2022) or very different from (Kong et al. 2018) the surface winds have been constructed to match the gravity field up to SH degree ∼10.

1.1. Jupiter's Surface Winds in Spectral Space

Here we point out that an important part of this debate stems from the limited spatial resolution of the previous gravity solutions. The typical latitudinal length scale of the observed surface winds corresponds to SH degrees much higher than 10 (Figure 1). As a result, slower and latitudinally smoother deep winds can be constructed if one only tries to match the large-scale gravity field with SH degree n ≤ 10.

Zoom In Zoom Out Reset image size

To illustrate this point, like the way we decompose the gravitational field onto a set of spherical harmonics, here we decompose the zonal wind angular velocity (in the rotating frame) onto the Legendre polynomials

$$\begin{eqnarray}&&\omega (\theta )=\displaystyle \sum _{l=0}{\omega }_{l}{P}_{l}(\cos \theta ),\end{eqnarray} \tag{ 1 }$$

where l is the SH degree, θ is the colatitude measured from the north pole, ${P}_{l}(\cos \theta )$ is the Legendre polynomial of degree l, and ω_l is the coefficient. Any residual solid-body rotation (SBR) component is described with the degree-0 coefficient ω₀. Even degree modes correspond to north–south symmetric winds, while odd degree modes correspond to north–south antisymmetric winds.

The zonal wind velocity U_ϕ in the rotating frame is connected to the angular velocity via

$$\begin{eqnarray}&&{U}_{\phi }(r,\theta )=\omega (r,\theta )s=\omega (r,\theta )r\sin \theta ,\end{eqnarray} \tag{ 2 }$$

where r is the spherical radial distance from the center of mass, and $s=r\sin \theta$ is the cylindrical radial distance to the spin axis.

Figure 1 shows Jupiter's average surface zonal wind angular velocity (Tollefson et al. 2017) in the System III rotating frame as a function of planetocentric latitude (panel (a)) as well as its decomposition onto the Legendre polynomials (panel (b)). It can be seen from Figure 1(b) that the spectral content of the surface zonal winds peak around SH degree 24. Gravitational harmonics up to a similar wavenumber would be needed to determine whether Jupiter's deep zonal wind pattern resembles that of the surface wind. Throughout this paper, we use l to refer to the SH degree of the zonal wind spectral decomposition and n to refer to the SH degree of the gravity harmonics.

The degree-0 coefficient ω₀ is nonzero for Jupiter's surface zonal wind when viewed in the System III rotating frame. This corresponds to a nonzero positive SBR component of Jupiter's surface zonal winds in this rotating frame defined by Jupiter's internal magnetic field. One can also see from Figure 1(a) that the surface zonal wind viewed in this frame has a net positive (eastward) component, e.g., the eastward winds are stronger and there are more eastward winds in particular at low latitudes.

1.2. An Updated Gravity Field for Jupiter

A new gravity field solution for Jupiter (Kaspi et al. 2023) was recently derived by analyzing the Juno Doppler tracking data up to orbit PJ 37 (Figure 2). In this new gravity solution, the small-scale gravitational accelerations beyond SH degree 12 were constrained to be zero near the poles (e.g., Konopliv et al. 2020; Park et al. 2020). More specifically, for latitudes between (90°S, 40°S) and (70°N, 90°N), we created a grid point for every two degrees in latitude. Then, for each grid point, we assumed the a priori surface acceleration values due to zonal harmonics J₁₃ to J₄₀ are zero. The a priori uncertainties for these surface acceleration points are empirically determined so that the mapped surface acceleration reaches the uncertainty 1 mGal (mili-Gal, where 1 Gal=1 cm s⁻²).

Zoom In Zoom Out Reset image size

Figure 2 shows the zonal harmonics of this new gravity solution, J_n , in black circles and the derived uncertainties in black dashed line up to SH degree 40. The filled (open) circles represent positive (negative) values. The derived uncertainties shown here are the square root of the diagonal terms of the full covariance matrix. The derived uncertainties are considered nonexplicit factors in the formal inversion process and are about 1.5 times the formal uncertainty (1σ). Thus, twice the derived uncertainty corresponds to about 3 times the formal uncertainty (3σ). It can be seen from Figure 2 that this new gravity solution resolves Jupiter's J_n up to SH degree 32 above the derived uncertainty, and up to SH degree 24 above twice the derived uncertainty. This new solution now resolves the length scales in Jupiter's gravity field that correspond to the typical (latitudinal) length scales in Jupiter's surface zonal flows (∼SH degree 24). Furthermore, this new high-resolution gravity solution affords an independent reconstruction/inversion of Jupiter's deep zonal winds without any input about the surface winds, which is the philosophy we adopt in this study. Our philosophy is in marked contrast to previous studies that project the observed surface winds into Jupiter's interior (e.g., Kaspi et al. 2018, 2023).

Contributions from the SBR dominate the even degree J_n of Jupiter up to SH degree ∼10 but become insignificant shortly after, as shown by the red circles in Figure 2, which correspond to the latest Jupiter interior model by Militzer et al. (2022). Most of the even degrees measured J_n with n ≥ 12 and all odd degrees J_n must have a dynamical origin (Hubbard 1999; Kaspi 2013). They correspond to axisymmetric density/shape perturbations that are most likely in balance with deep zonal flows. It shall be noted that the details of this density-flow balance remain somewhat uncertain, in particular for the part where the wind decays with depth (e.g., see the discussion in Kulowski et al. 2021). At present, the only dynamical model for deep wind truncation inside Jupiter (Christensen et al. 2020) asserts that the zonal winds are z-invariant from the near-surface until a truncation depth where the wind amplitude decays sharply within a vertical range of 150–300 km via a thermal wind shear mechanism. Here, the z-axis is aligned with the spin axis of Jupiter. The latitudinal density gradient in this model (Christensen et al. 2020) results from meridional circulation acting on the background nonadiabatic density gradient in a hypothesized stably stratified layer (SSL). Maxwell stress in the semiconducting region could provide the driving force of the meridional circulation (Christensen et al. 2020). Reynolds stress is likely another option.

Even with ∼30 gravitational harmonics, inferring the deep wind structure inside Jupiter is still a nonunique problem given the integral nature of the gravitational harmonics. A priori assumptions about the spatial profile of the wind (latitudinal and/or vertical) and the force/vorticity balance are necessary ingredients of the solution, regardless of whether they are made explicit or not. Here we adopt the Christensen et al. (2020) deep zonal wind decay mechanism but do not require the latitudinal profile of the deep zonal winds to be that of the observed surface wind. In fact, no information on the observed surface zonal winds was fed into our inversion analysis. In addition, we allow the truncation depth of the deep zonal winds to be a free parameter and scan the range between 1000 and 3000 km. By adopting this philosophy, we could utilize the latest Jupiter gravity solution to clarify the resemblance between the surface and deep zonal winds above the highly conducting dynamo region (Bloxham et al. 2022).

2.1. Construction of an Ensemble of Deep Zonal Wind Modes

We first construct an ensemble of deep zonal wind modes with Legendre polynomials in the latitudinal direction (Equation (3)) and hyperbolic tangent functions in the vertical direction (Equation (4)) and then forward compute the J_n associated with these wind modes with the thermal wind balance (see Section 2.2 for more details regarding the thermal wind balance adopted in this study). Guided by the Christensen et al. (2020) model, the thickness of the (vertical) transition layer is kept in the range of 150–300 km. Figures 3(a) and (b) showcase two examples of individual deep zonal wind modes with truncation depth of 2500 km, corresponding to SH degrees 20 and 21, respectively.

Zoom In Zoom Out Reset image size

The construction of an individual deep zonal wind mode starts with a surface wind angular velocity prescribed by a single Legendre polynomial

$$\begin{eqnarray}&&\omega ({\theta }_{S})={\omega }_{l}{P}_{l}(\cos {\theta }_{S}),\end{eqnarray} \tag{ 3 }$$

where θ_S is the colatitude at the surface of the planet and ω_l is the coefficient. For the north–south symmetric wind modes, which are even degree modes, we first extend the surface winds along the spin-axis direction (z) invariantly and then apply the following hyperbolic vertical truncation function

$$\begin{eqnarray}&&f(r)=\displaystyle \frac{1}{2}\displaystyle \frac{\tanh \left(\tfrac{r+H-{r}_{J}}{{\rm{\Delta }}H}+1\right)+1}{\tanh \left(\tfrac{H}{{{\rm{\Delta }}}_{H}}+2\right)},\end{eqnarray} \tag{ 4 }$$

where tanh is the hyperbolic tangent function, r_J is the radius of Jupiter, H is the truncation depth, and ΔH is the parameter setting the thickness of the transition layer. ΔH is set to values between 50 and 100 km for the cases presented in this study to achieve a transition thickness between 150 and 300 km as suggested by the theoretical model of Christensen et al. (2020). Figure 3(c) shows an example of the vertical truncation function with H = 2500 km and ΔH = 75 km.

For north–south antisymmetric wind modes, which are the odd degree modes, we extend the surface winds along the spin-axis direction (z) invariantly in their respective hemisphere, then apply the following equatorial smoothing function, and further apply the vertical truncation function (Equation (4)),

$$\begin{eqnarray}&&\begin{array}{c}\omega (s,z)={\omega }_{N}({s}_{\mathrm{surf}}=s)\frac{1+\mathrm{erf}\left(\frac{z}{{\rm{\Delta }}z}\right)}{2}\\ \qquad +\ {\omega }_{S}({s}_{\mathrm{surf}}=s)\frac{1-\mathrm{erf}\left(\frac{z}{{\rm{\Delta }}z}\right)}{2},\end{array}\end{eqnarray} \tag{ 5 }$$

where ω_N is the surface angular velocity in the northern hemisphere and ω_S is the surface angular velocity in the southern hemisphere. This ensures that the zonal wind in each hemisphere is nearly z-invariant except within a thin layer near the equator. The influence of the choice of Δz on the values of odd J_n up to SH degree 40 associated with each individual mode is minimal for all Δz < 0.08r_J . In all the cases presented in this study, Δz is set to 0.075 r_J . Moreover, the reconstructed deep winds are found to be mostly north–south symmetric within ±15° from the equator, which implies that solutions with a truncation depth ≤2500 km are essentially unaffected by the equatorial smoothing function. Figure 4 shows an example of the equatorial smoothing function corresponding to 30° surface latitude. Here ω_N is normalized to 1 while ω_S is normalized to −1.

Zoom In Zoom Out Reset image size

2.2. Formulate the Deep Wind–Gravity Field Connection as a Linear Inverse Problem

For each deep zonal wind mode, referred to as ω_l for simplicity, we can forward compute the associated gravitational harmonics ${J}_{n}^{{\omega }_{l}}$ for all n of interest, which are dominated by the ones with the same symmetry properties around SH degree l. n_max is set to 40 as constrained by the available gravity solution, and we tested l_max between 40 and 60 and obtained broadly similar results. In all results presented in this paper, l_max is set to 50. In our forward computation, we adopt the thermal wind balance (Kaspi 2013; Cao & Stevenson 2017b; Kaspi et al. 2018) with spherical background density ρ₀(r) and internal gravitational acceleration g ₀(r) provided by the latest Jupiter interior model (Militzer et al. 2022). More specifically, we adopt the following thermal wind balance for both the z-invariant and the z-varying part of the deep zonal wind

$$\begin{eqnarray}&&2| {{\rm{\Omega }}}_{0}| \displaystyle \frac{\partial [{\rho }_{0}(r){U}_{\phi }{{\boldsymbol{e}}}_{\phi }]}{\partial z}=-\displaystyle \frac{1}{r}\displaystyle \frac{\partial \rho ^{\prime} }{\partial \theta }{{\boldsymbol{e}}}_{\theta }\times [-| {g}_{0}(r)| {{\boldsymbol{e}}}_{r}],\end{eqnarray} \tag{ 6 }$$

in which $\rho ^{\prime}$ is the wind-induced density perturbation and e _r , e _θ , and e _ϕ are unit vectors in the spherical radial, colatitudinal, and azimuthal directions, respectively. This choice is based on the theoretical work of Christensen et al. (2020) and supported by the 3D numerical modeling of Gastine & Wicht (2021), in which the z-varying part of the zonal winds is balanced by a latitudinal gradient in the nonadiabatic density (thermal wind shear). This is different from the barotropic zonal wind model adopted by Kulowski et al. (2021), in which Equation (6) is only applied to the z-invariant part of the zonal winds.

With a given deep zonal wind profile, U_ϕ (r, θ), one can integrate Equation (6) in the θ direction to get the wind-induced density perturbation, $\rho ^{\prime} (r,\theta )$, and then compute the wind-induced gravity harmonics via the volume integral

$$\begin{eqnarray}&&{J}_{n}=-\displaystyle \frac{1}{{{Ma}}^{n}}{\int }_{V}\rho ^{\prime} ({\boldsymbol{r}}){r}^{n}{P}_{n}(\cos \theta ){dV},\end{eqnarray} \tag{ 7 }$$

where M is the total mass of Jupiter, a is the 1 bar equatorial radius of Jupiter, and V is the entire volume of Jupiter. Interested readers can refer to Cao & Stevenson (2017b) for the derivation of Equation (6) and related technical details.

For each truncation depth, we can then formulate the deep wind–gravity field connection as a linear inverse problem. Formally, one could write the forward problem as

$$\begin{eqnarray}&&\mathrm{obs}=G\,\mathrm{model};\end{eqnarray} \tag{ 8 }$$

here the observations, obs, are the ${\rm{\Delta }}{J}_{n}^{\mathrm{obs}}$, which is the difference between the measured J_n and those associated with SBR J_n ^SBR, the model parameters, model, are the coefficients ω_l associated with the deep zonal winds, while the forward matrix, G, consists of the (forward computed) wind-induced gravity harmonics, ${J}_{n}^{{\omega }_{l}}$, associated with the individual zonal wind mode ω_l . To be more specific, one can write the above forward problem as

$$\begin{eqnarray}&&\left[\begin{array}{l}{\rm{\Delta }}{J}_{2}\\ {\rm{\Delta }}{J}_{3}\\ {\rm{\Delta }}{J}_{4}\\ \vdots \\ {\rm{\Delta }}{J}_{n}\end{array}\right]=\left[\begin{array}{ccccc}{J}_{2}^{{\omega }_{0}} & {J}_{2}^{{\omega }_{1}} & {J}_{2}^{{\omega }_{2}} & ... & {J}_{2}^{{\omega }_{l}}\\ {J}_{3}^{{\omega }_{0}} & {J}_{3}^{{\omega }_{1}} & {J}_{3}^{{\omega }_{2}} & ... & {J}_{3}^{{\omega }_{l}}\\ {J}_{4}^{{\omega }_{0}} & {J}_{4}^{{\omega }_{1}} & {J}_{4}^{{\omega }_{2}} & ... & {J}_{4}^{{\omega }_{l}}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ {J}_{n}^{{\omega }_{0}} & {J}_{n}^{{\omega }_{1}} & {J}_{n}^{{\omega }_{2}} & ... & {J}_{n}^{{\omega }_{l}}\end{array}\right]\left[\begin{array}{c}{\omega }_{0}\\ {\omega }_{1}\\ {\omega }_{2}\\ \vdots \\ {\omega }_{l}\end{array}\right],\end{eqnarray} \tag{ 9 }$$

for which the inverse can be computed with standard linear inversion techniques such as the singular value decomposition (SVD).

Formulating the deep wind–gravity field connection as a linear inverse problem also allows us to quantify the uncertainties of the solution based on the uncertainties in the observations. The error propagation starts with the full covariance matrix of the derived J_n up to SH degree 40, which we denote as data_cov. The model covariance matrix, mod_cov, which describes the uncertainties and the correlations of the model parameters can be computed via

$$\begin{eqnarray}&&{\mathrm{mod}}_{\mathrm{cov}}\,=\,{G}^{-1}\,{\mathrm{data}}_{\mathrm{cov}}{\left({G}^{-1}\right)}^{T},\end{eqnarray} \tag{ 10 }$$

where G⁻¹ is the inverse matrix and the superscript T represents matrix transpose. Here, the model parameters are the coefficients ω_l associated with the reconstructed zonal winds.

One can then compute the uncertainties associated with the solution in real space, Δω(θ), via

$$\begin{eqnarray}&&{\left[{\rm{\Delta }}\omega ({\theta }_{S})\right]}^{2}=\mathrm{Fwd}({\theta }_{S})\,{\mathrm{mod}}_{\mathrm{cov}}{\left[\mathrm{Fwd}({\theta }_{S})\right]}^{T},\end{eqnarray} \tag{ 11 }$$

where Fwd is the forward matrix connecting ω(θ) and ω_l evaluated at θ_S

$$\begin{eqnarray}&&\mathrm{Fwd}({\theta }_{S})=\left[{P}_{0}(\cos {\theta }_{S}),{P}_{1}(\cos {\theta }_{S}),...,{P}_{l}(\cos {\theta }_{S})\right].\end{eqnarray} \tag{ 12 }$$

With our methodology, we can solve for Jupiter's deep zonal winds at any truncation depth from a given set of ${\rm{\Delta }}{J}_{n}^{\mathrm{obs}}$ by inverting Equation (9). We remove the SBR gravitational harmonics, J_n ^SBR, corresponding to the Militzer et al. (2022) Jupiter model up to SH degree 16 from the observed J_n of Jupiter (Kaspi et al. 2023) to get the ${\rm{\Delta }}{J}_{n}^{\mathrm{obs}}$. There are still some uncertainties in the SBR J_n of Jupiter, in particular for J₆ (see Miguel et al. 2022). The Militzer et al. (2022) Jupiter structural model with a dilute core requires a dynamical contribution to J₆ of −0.27 × 10⁻⁶. However, low-degree ΔJ_n , even or odd, correspond to long-wavelength features in the deep zonal flows and have a relatively weak connection to the observed narrow zonal jets off the equator (see Figure 1).

Figure 5 compares the angular velocity of the reconstructed deep zonal winds with a truncation depth of 2500 km (solid red line) to the observed surface zonal winds (solid blue line) as a function of surface planetocentric latitude. Also shown are the uncertainties associated with this reconstruction (magenta dashed lines), computed/propagated from the full covariance matrix of the gravity solution following Equations (10) and (11). The uncertainties here correspond to twice the derived uncertainties in J_n (∼3 times the formal uncertainties). It can be seen that there is a strong resemblance between the nominal reconstructed deep zonal winds and the observed surface winds, in particular the northern off-equatorial jet (NOEJ) around 20°N and the southern off-equatorial jet (SOEJ) at a similar latitude in the southern hemisphere (Kulowski et al. 2021). For this particular truncation depth, the amplitude of the reconstructed prograde NOEJ also matches that of the observed NOEJ, which could indicate a local barotropic atmosphere at this location. The strong resemblance between the reconstructed deep wind (the nominal solution with H = 2500 km) and the observed surface wind is mostly confined within ±35° latitude from the equator. Another feature of the reconstructed deep winds is the more pronounced north–south symmetry of the first retrograde jet adjacent to the prograde equatorial jet compared to the observed surface winds at Jupiter. This feature resembles the observed surface winds and inferred deep winds at Saturn (García-Melendo et al. 2011; Galanti et al. 2019; Militzer et al. 2019), as well as nearly all numerical deep convection models for Jupiter (and Saturn; e.g., Heimpel et al. 2005; Gastine & Wicht 2021). It is interesting to note that, while the pattern of the reconstructed mid-to-high latitude winds seems to deviate from that observed, the amplitude of the two remain very close to each other. Here we emphasize that no information about the observed surface winds was used in the inversion. Our reconstructed wind is an emergent solution from inverting Equation (9) with the gravity harmonics of Jupiter and our forward model.

Zoom In Zoom Out Reset image size

Two aspects of the uncertainties associated with the reconstructed deep winds are worthy of comment. First, the uncertainties in the southern hemisphere are appreciably larger than that in the northern hemisphere. This results from the north–south asymmetric nature of Juno's orbits with increasingly more low-altitude coverage in the northern hemisphere due to the northward precession of Juno's periapsis. Second, the uncertainties within ±12° of the equator are as large as those in the southern midlatitudes. This is due to geometric effects associated with the shape of the planet and z-projection of the zonal winds, as (1) the z-derivative of the background density decreases toward zero as one approaches the lower latitude and (2) less mass is involved in a z-aligned column (annulus) with the same horizontal width at lower latitude. The observed prograde equatorial jet is within the uncertainties of the reconstructed deep zonal winds. Furthermore, since we are working with a spherical harmonic reconstruction of the deep winds and our knowledge of Jupiter's gravity harmonics (and their uncertainties) are currently limited to SH degree 40 (instead of, e.g., SH degree 100), these limit the computed uncertainties very close to the equator (e.g., within a few degrees from the equator). The true uncertainties very close to the equator are expected to be higher, as the geometric effect monotonically strengthens toward the equator.

As our wind reconstruction procedure is not bounded by the observed surface winds, we could also investigate the effect of the truncation depth on the deep zonal flow amplitude profile. Figure 6 displays an ensemble of reconstructed deep zonal wind as a function of the truncation depth. It can be seen that the latitudinal profile of the reconstructed deep zonal winds with different truncation depths remain broadly similar, in particular at latitudes northward of ∼30°S and the prominent prograde NOEJ (keeping in mind that the reconstruction features higher uncertainties in the southern hemisphere). For truncation depth > 2500 km, the amplitude of the reconstructed wind remains on the order of 100 m s⁻¹. The deeper-most portion of these winds would violate the ohmic dissipation constraint that results from the rapidly increasing electrical conductivity as a function of depth (Liu et al. 2008; Cao & Stevenson 2017a). The amplitude of the reconstructed deep winds increases as the truncation depth decreases, with faster rates at shallower truncation depth. The peak amplitude of the prograde NOEJ reaches ∼280 m s⁻¹, which is about twice that of the surface NOEJ, when the truncation depth is ∼1500 km. The peak amplitude reaches ∼780 m s⁻¹ (480 m s⁻¹) for a truncation depth of 1000 km (1250 km), due to the rapidly decreasing background density at shallower depth. For comparison, the sound speed at 1 bar is ∼1 km s⁻¹ (Lorenz 1998), and increases to ∼4 km s⁻¹ at a depth ∼1000 km (French et al. 2012).

Zoom In Zoom Out Reset image size

The Galileo probe, which was dropped into Jupiter around 6.5°N and 4.4°W (measured in the Jupiter System III coordinate) on 1995 December 7 measured a factor of 2 increase in local zonal wind speed, from ∼85 m s⁻¹ at 1 bar to ∼170 m s⁻¹ at 5 bar, which then remains at the high-speed level until about 21 bar (Atkinson et al. 1997). If similar near-surface shear also prevails at other latitudes (in particular near the NOEJ), our analysis indicates that the depth of such stronger winds would be shallower, ∼1500 km, though still significantly deeper than the expected water cloud layers. As discussed earlier, deep winds very close to the equator have relatively weak contributions to the gravity field due to geometrical effects. As a result, the true uncertainties of the derived deep winds very close to the equator, e.g., at the Galileo probe latitude, are expected to be large. Moreover, the uncertainties in the SBR contribution to low-degree even J_ns can affect the value of the derived equatorial jet as well as the pattern of the large-scale deep zonal wind in the latitudinal direction.

Adopting vertical profiles guided by the Christensen model (Christensen et al. 2020; Gastine & Wicht 2021) but making no a priori assumptions about the latitudinal profile, we constructed a suite of deep zonal winds inside Jupiter based on the latest gravity solution up to SH degree 40 (Kaspi et al. 2023). The latitudinal profile of our reconstructed deep zonal flows bears a strong resemblance to the observed surface zonal winds within ± 35° latitude from the equator, in particular the prograde NOEJ around 20° N. The NOEJ has been shown to dominate the low-degree odd gravity harmonics (J₅, J₇, J₉) even when assuming the vertical wind shear is balanced by Reynolds/Maxwell stress instead of latitudinal entropy/composition shear (thermal wind balance) by Kulowski et al. (2021). Moreover, our reconstruction demonstrated that latitudinally narrow zonal flows on the order of 100 m s⁻¹ or stronger reside in the top ∼2000 km of Jupiter. Our results thus strongly support the picture that the amplitude and pattern of the deep zonal winds in the nonelectrically conducting part of the molecular envelope, which extends much deeper than the water cloud layer, are comparable to those observed at the surface. It remains a critical open question whether this picture could be reconciled with dynamical models without an SSL at the truncation depth. At present, no structure models require/invoke an SSL near this depth (e.g., Debras & Chabrier 2019; Miguel et al. 2022; Militzer et al. 2022). Answering this question has important implications not only for rotating fluid dynamics but also for the structure and evolution of Jupiter. A near-surface SSL can significantly impact the temperature profile in the deeper part of Jupiter as well as the cooling efficiency of the entire planet.

All authors acknowledge support from the NASA Juno project. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. H.C. would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the program Frontiers in dynamo theory: from the Earth to the stars where some work on this paper was undertaken. This work was supported by EPSRC grant No. EP/R014604/1. In addition, this work was partially supported by a grant from the Simons Foundation.

All the Juno radio tracking data used to derive the gravity moments and their uncertainties in this study are available from the NASA Planetary Data System (Buccino et al. 2020; https://pds.nasa.gov).

H.C. conceptualized the study, carried out the reconstruction of deep winds from gravity harmonics, and drafted the manuscript. J.B., R.K.Y., and L.K. supported the analysis and interpretation of the results. R.S.P. derived the gravity harmonics of Jupiter and their uncertainties from the Juno Radio Tracking experiments. B.M. provided the gravity harmonics associated with the solid-body rotating Jupiter. All authors contributed to the discussion, as well as editing and revising the manuscript. D.J.S. leads the Interior Working Group of Juno, and S.J.B. is the principal investigator of the Juno Mission.