Strong resemblance between surface and deep zonal winds inside Jupiter
revealed by high-degree gravity moments

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 (Fig. 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\leq 10$.

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

where $l$ is the SH degree, $\theta$ is the co-latitude measured from the north-pole, $P_{l}(\cos\theta)$ is the Legendre polynomial of degree $l$, and $\omega_{l}$ is the coefficient. Any residual solid-body rotation component is described with the degree-0 coefficient $\omega_{0}$. Even degree modes correspond to north-south symmetric winds, while odd degree modes correspond to north-south antisymmetric winds.

The zonal wind velocity $U_{\phi}$ in the rotating frame is connected to the angular velocity via

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.

Fig. 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 Fig. 1b that the spectral content of the surface zonal winds peak around SH degree 24. Gravitational harmonics up to a similar wave-number 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 $\omega_{0}$ is non-zero for Jupiter’s surface zonal wind when viewed in the System III rotating frame. This corresponds to a non-zero positive solid-body rotation component of Jupiter’s surface zonal winds in this rotating frame defined by Jupiter’s internal magnetic field. One can also see from Fig. 1a 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.

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 (Fig. 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_{13}$ to $J_{40}$ 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^{2}$).

Fig. 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 considered non-explicit factors in the formal inversion process, and are about 1.5 times the formal uncertainty ($1\sigma$). Thus, twice the derived uncertainty corresponds to about three times the formal uncertainty ($3\sigma$). It can be seen from Fig. 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 ($\sim$ 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 which project the observed surface winds into Jupiter’s interior (e.g., Kaspi et al., 2018, 2023).

Contributions from the solid-body-rotation (SBR) dominate the even degree $J_{n}$ of Jupiter up to SH degree $\sim$10 but become insignificant shortly after, as shown by the red circles in Fig. 2 which correspond to the latest Jupiter interior model by Militzer et al. (2022). Most of the even degree measured $J_{n}$ with $n\geq 12$ and all odd degrees $J_{n}$ must have a dynamical origin (Hubbard, 1999; Kaspi, 2013). They correspond to axisymmetric density/shape perturbations which 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 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 non-adiabatic density gradient in a hypothesized stably stratified layer (SSL). Maxwell stress in the semi-conducting region could provide the driving force of the meridional circulation (Christensen et al., 2020). Reynolds stress is likely another option.

We first construct an ensemble of deep zonal wind modes with Legendre polynomials in the latitudinal direction (Eq. 4) and hyperbolic tangent functions in the vertical direction (Eq. 5) 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. Fig. 3ab showcase two examples of individual deep zonal wind mode with truncation depth of 2500 km, corresponding to SH degree 20 and 21, respectively.

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

where $\theta_{S}$ is the co-latitude at the surface of the planet and $\omega_{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

where $\tanh$ is the hyperbolic tangent function, $r_{J}$ is the radius of Jupiter, $H$ is the truncation depth, and $\Delta H$ is the parameter setting the thickness of the transition layer. $\Delta H$ is set to values between 50 km 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). Fig. 3c shows an example of the vertical truncation function with $H=2500$ km and $\Delta 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 (Eq. 5),

where $\omega_{N}$ is the surface angular velocity in the northern hemisphere and $\omega_{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 $\Delta z$ on the values of odd $J_{n}$ up to SH degree 40 associated with each invidual mode is minimal for all $\Delta z<0.08r_{J}$. In all the cases presented in this study, $\Delta z$ is set to 0.075 $r_{J}$. Moreover, the reconstructed deep winds are found to be mostly north-south symmetric within $\pm$15^∘ from the equator, which implies that solutions with a truncation depth $\leq$ 2500 km are essentially unaffected by the equatorial smoothing function. Fig. 4 shows an example of the equatorial smoothing function corresponding to 30^∘ surface latitude. Here $\omega_{N}$ is normalized to 1 while $\omega_{S}$ is normalized to -1.

For each deep zonal wind mode, referred to as $\omega_{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 $\rho_{0}(r)$ and internal gravitational acceleration $\mathbf{g_{0}}(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

in which $\rho^{\prime}$ is the wind-induced density perturbation and $\mathbf{e}_{r}$, $\mathbf{e}_{\theta}$, and $\mathbf{e}_{\phi}$ are unit vectors in the spherical radial, co-latitudinal, 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 non-adiabatic density (thermal wind shear). This is different from the barotropic zonal wind model adopted by Kulowski et al. (2021), in which Eq. (7) is only applied to the $z$-invariant part of the zonal winds.

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

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 eqn. (7) 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

here the observations, $\mathbf{obs}$, are the $\Delta J_{n}^{obs}$, which is the difference between the measured $J_{n}$ and those associated with solid-body-rotation $J_{n}^{SBR}$, the model parameters, $\mathbf{model}$, are the the coefficients $\omega_{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 $\omega_{l}$. To be more specific, one can write the above forward problem as

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

where $G^{-1}$ is the inverse matrix and the superscript $T$ represents matrix transpose. Here, the model parameters are the coefficients $\omega_{l}$ associated with the reconstructed zonal winds.

One can then compute the uncertainties associated with the solution in real space, $\Delta\omega(\theta)$, via

where Fwd is the forward matrix connecting $\omega(\theta)$ and $\omega_{l}$ evaluated at $\theta_{S}$