Linking Zonal Winds and Gravity II: explaining the equatorially antisymmetric gravity moments of Jupiter
Wieland Dietrich, Paula Wulff, Johannes Wicht, Ulrich R. Christensen
MMNRAS , 1– ?? (202-) Preprint 26 February 2021 Compiled using MNRAS L A TEX style file v3.0
Linking Zonal Winds and Gravity II: explaining the equatoriallyantisymmetric gravity moments of Jupiter
Wieland Dietrich ★ , Paula Wulff , , Johannes Wicht , and Ulrich R. Christensen Max Planck Institute for Solar System Research, Justus-von-Liebig-Weg 3, 37077, Göttingen, Germany Georg August University, Institute for Geophysics, Friedrich-Hund-Platz 1, 37077 Göttingen
Last updated 26 February 2021; in original form 26 February 2021
ABSTRACT
The recent gravity field measurements of Jupiter (Juno) and Saturn (Cassini) confirm theexistence of deep zonal flows reaching to a depth of 5% and 15% of the respective radius.Relating the zonal wind induced density perturbations to the gravity moments has becomea major tool to characterise the interior dynamics of gas giants. Previous studies differ withrespect to the assumptions made on how the wind velocity relates to density anomalies, onthe functional form of its decay with depth, and on the continuity of antisymmetric windsacross the equatorial plane. Most of the suggested vertical structures exhibit a rather smoothradial decay of the zonal wind, which seems at odds with the observed secular variation of themagnetic field and the prevailing geostrophy of the zonal winds. Moreover, the results reliedon an artificial equatorial regularisation or ignored the equatorial discontinuity altogether.We favour an alternative structure, where the equatorially antisymmetric zonal wind in anequatorial latitude belt between ± ◦ remains so shallow that it does not contribute to thegravity signal. The winds at higher latitudes suffice to convincingly explain the measuredgravity moments. Our results indicate that the winds are geostrophic, i.e. constant alongcylinders, in the outer 3000 km and decay rapidly below. The preferred wind structure is 50%deeper than previously thought, agrees with the measured gravity moment, is compliant withthe magnetic constraints and the requirement of an adiabatic atmosphere and unbiased by thetreatment of the equatorial discontinuity. Key words: planets and satellites: gaseous planets – gravitation – hydrodynamics
The spacecrafts
Juno and
Cassini delivered high-accuracy mea-surements of the gravity potentials for Jupiter and Saturn (Kaspiet al. 2018; Iess et al. 2018; Guillot et al. 2018; Galanti et al. 2019;Iess et al. 2019), which provide valuable constraints on the interiorstructure and dynamics of their atmospheres. For the first time, ithas been possible to resolve the tiny undulations in the gravity po-tential induced by zonal winds. This has ended the long-standingdebate on whether the zonal winds on Jupiter and Saturn are shallowweather phenomena or reach deeper into the planets’ convective en-velopes. The gravity data suggest that the winds extend down to 5%and 15% of Jupiter’s and Saturn’s radii, respectively (Kaspi et al.2018; Galanti et al. 2019; Kaspi et al. 2020). However, modellingthe deep-reaching zonal mass fluxes on a gas planet with all theircomplexity and relating them to anomalies in the gravity field is adifficult, non-unique problem.The recent attempts of constraining the deep-reaching windswith gravity data are based on simple parametrizations of the ver-tical wind structure. The observed surface winds are thus con- ★ Contact e-mail: [email protected] tinued downward along cylinders that are aligned with the rota-tion axis. The resulting geostrophic flow is then multiplied with aparametrized radial profile to model the decay with depth. The pa-rameters of this profile are constrained by the gravity measurements(Kaspi et al. 2010, 2013, 2018; Kong et al. 2018; Galanti et al. 2019).More recently, constraints from Jupiter’s magnetic field (Duer et al.2019) or it’s temporal evolution (Galanti & Kaspi 2020) have alsobeen used in addition to this. However, some studies found it neces-sary to partially modify the zonal flow compared to what is observedat cloud level (Kong et al. 2018; Galanti & Kaspi 2020) in order tomatch the modelled gravity anomalies to the observed ones.Gravity moments of odd degree exclusively contain the im-pact of equatorially antisymmetric zonal flows, since they are notobscured by the rotational deformation of the planet as opposed totheir even counterparts. The geostrophic, downward continuationof the antisymmetric zonal flows from each hemisphere yields op-posite signs when they each reach the equatorial plane and henceintroduces a problematic discontinuity or equatorial ‘step‘. This hasbeen reported to potentially bias the results (Kong et al. 2016). Tocircumvent this issue, several published studies, such as Kaspi et al.(2018); Galanti & Kaspi (2020) calculate the antisymmetric grav-ity signal in each hemisphere, but ignore the resulting equatorial © a r X i v : . [ a s t r o - ph . E P ] F e b W. Dietrich et al.
Figure 1.
Various proposed radial decay functions alongside the expectedamplitude of zonal flows that are compliant with the secular variation ofJupiter’s magnetic field (Moore et al. 2019). shear region. Alternatively Kong et al. (2017); Kong et al. (2018)suggested to smooth the equatorial region artificially.Zonal winds induce perturbations in the gravity potential viapressure perturbations that change the density distribution. The gov-erning leading order Navier-Stokes (or Euler) equation, a balancebetween pressure gradient; Coriolis force; and gravity forces; es-tablishes the link between zonal flows and the induced densityperturbation. The gravity term has two contributions, one due tothe zonal-flow induced density perturbation and a second due tothe dynamic gravity perturbation. If the latter, termed dynamic selfgravity (DSG) by Wicht et al. (2020), is ignored, the balance reducesto the classic thermal wind equation (TWE) which can directly besolved for the density perturbation and hence the gravity pertur-bation. This approach has been commonly used for modelling thegravity anomalies of Jupiter and Saturn and characterising theirdeep interior structure (Guillot et al. 2018; Kaspi et al. 2018; Iesset al. 2019; Galanti et al. 2019). Several authors claim that the DSG-term is indeed negligible (Galanti et al. 2017; Kaspi et al. 2018), butthis has been disputed by Zhang et al. (2015), Kong et al. (2017),and Wicht et al. (2020).Retaining the DSG-related term yields the so called thermo-gravitational wind equation (TGWE), which is harder to handlemathematically and numerically. Zhang et al. (2015) reformulatethe TGWE into an integro-differential equation for the density per-turbation. Their approach is restricted to polytropic interiors withindex unity. More recently, Wicht et al. (2020) show that the TGWEcan be formulated as an inhomogeneous Helmholtz equation for thegravity potential. The results show that the relative impact of theDSG decreases with the spherical harmonic degree ℓ of the gravityperturbation. Being of order one for ℓ =
2, it decreases to 10% atabout ℓ = ℓ >
25. However, thepolytrope provides only a crude approximation of Jupiter’s interior.Here we introduce a solution method for the TGWE that is basedon formalism of Wicht et al. (2020) but can also handle more real-istic interior models that are reflected by a radius-dependent DSGcoefficient.The discrepancy between existing studies of Jupiter’s odd grav-ity moments is related to differences in the modelling assumptions.Kaspi et al. (2018) explained the first four odd gravity momentsof Jupiter by extending the observed zonal wind profile obtained during the perijove (Tollefson et al. 2017) along cylinders, invert-ing the TWE equation, assuming a realistic Jupiter interior model(Guillot et al. 2003) and a relatively smooth radial decay while ig-noring the discontinuity at the equatorial plane (see fig. 1, green).This has been challenged by Kong et al. (2017), arguing that thethermo-gravitational wind equation (TGWE) should be used. Thesubsequent application to Jupiter by Kong et al. (2018) is based ona polytropic interior and the surface flow measurements from theCassini mission (Porco et al. 2003). In order to carry this out thesurface zonal flow was also modified by putting a cap on the ampli-tude of the spike in the antisymmetric flow at 𝛽 = ◦ latitude, andapplied a second decay function to smooth out the discontinuity atthe equatorial plane (orange profile).Both studies favour a rather smooth decay with depth in order tomatch the observed gravity signal, which implies quite significantflow amplitudes below 0 . 𝑅 . The radial decay profile of Kaspiet al. (2018) implies a remaining 10% flow amplitude at 0 . 𝑅 ,whereas Kong et al. (2018) found more than half of the surfaceamplitude remaining at this depth. The weak secular variation ofJupiter’s observed magnetic field over the past decades, however,can only be explained if the amplitude of the zonal wind at thedepth of the horizontal advection of the radial field, i.e. 4% to 6%of Jupiter’s radius, is as weak as cm/s compared to tens or hundredsof meters per second at the surface (Moore et al. 2019) (blackdots in the figure). More recently, Galanti & Kaspi (2020) showedthat a combined gravito-magnetic analysis favoured a sharper radialprofile with a 2000 km deep barotropic part and 600 km decay region(red profile in fig. 1). However, also this analysis neglected the flowdiscontinuity at the equatorial plane.In this study we revisit the key model assumptions and proper-ties on which the analyses of odd-degree gravity moments are basedand quantify their impact on the solutions. This includes justifyingthe fundamental equation, the problematic treatment of the equator,various surface zonal flows profiles and different models of Jupiter’sinterior. The gravity potential 𝛹 is directly related to the density distribution 𝜌 via the Poisson equation ∇ 𝛹 = 𝜋𝐺 𝜌 , (1)where 𝐺 is the gravitational constant. The general solution is 𝛹 ( 𝑟, 𝜃 ) = 𝜋𝐺 ∫ 𝑟 ∫ 𝜃 𝜌 ( ˜ 𝑟, ˜ 𝜃 )| 𝒓 − ˜ 𝒓 | ˜ 𝑟 sin ˜ 𝜃 𝑑 ˜ 𝑟 𝑑 ˜ 𝜃 (2)and the associated gravity force is 𝒈 = − ∇ 𝛹 . It is useful to separatethe external gravity potential into a spherically symmetric part ofa non-rotating, solid body and a series of higher order terms origi-nating from density perturbations and the rotational deformation ofthe planet: 𝛹 = − 𝐺 𝑀𝑟 (cid:34) − ∞ ∑︁ ℓ = 𝐽 ℓ (cid:18) 𝑅𝑟 (cid:19) ℓ 𝑃 ℓ ( cos 𝜃 ) (cid:35) , (3)where 𝑅 is the equatorial planetary radius, 𝑀 the total mass and 𝑃 ℓ are the Legendre polynomials of degree ℓ . The ℓ = 𝐽 ℓ MNRAS , 1– ????
25. However, thepolytrope provides only a crude approximation of Jupiter’s interior.Here we introduce a solution method for the TGWE that is basedon formalism of Wicht et al. (2020) but can also handle more real-istic interior models that are reflected by a radius-dependent DSGcoefficient.The discrepancy between existing studies of Jupiter’s odd grav-ity moments is related to differences in the modelling assumptions.Kaspi et al. (2018) explained the first four odd gravity momentsof Jupiter by extending the observed zonal wind profile obtained during the perijove (Tollefson et al. 2017) along cylinders, invert-ing the TWE equation, assuming a realistic Jupiter interior model(Guillot et al. 2003) and a relatively smooth radial decay while ig-noring the discontinuity at the equatorial plane (see fig. 1, green).This has been challenged by Kong et al. (2017), arguing that thethermo-gravitational wind equation (TGWE) should be used. Thesubsequent application to Jupiter by Kong et al. (2018) is based ona polytropic interior and the surface flow measurements from theCassini mission (Porco et al. 2003). In order to carry this out thesurface zonal flow was also modified by putting a cap on the ampli-tude of the spike in the antisymmetric flow at 𝛽 = ◦ latitude, andapplied a second decay function to smooth out the discontinuity atthe equatorial plane (orange profile).Both studies favour a rather smooth decay with depth in order tomatch the observed gravity signal, which implies quite significantflow amplitudes below 0 . 𝑅 . The radial decay profile of Kaspiet al. (2018) implies a remaining 10% flow amplitude at 0 . 𝑅 ,whereas Kong et al. (2018) found more than half of the surfaceamplitude remaining at this depth. The weak secular variation ofJupiter’s observed magnetic field over the past decades, however,can only be explained if the amplitude of the zonal wind at thedepth of the horizontal advection of the radial field, i.e. 4% to 6%of Jupiter’s radius, is as weak as cm/s compared to tens or hundredsof meters per second at the surface (Moore et al. 2019) (blackdots in the figure). More recently, Galanti & Kaspi (2020) showedthat a combined gravito-magnetic analysis favoured a sharper radialprofile with a 2000 km deep barotropic part and 600 km decay region(red profile in fig. 1). However, also this analysis neglected the flowdiscontinuity at the equatorial plane.In this study we revisit the key model assumptions and proper-ties on which the analyses of odd-degree gravity moments are basedand quantify their impact on the solutions. This includes justifyingthe fundamental equation, the problematic treatment of the equator,various surface zonal flows profiles and different models of Jupiter’sinterior. The gravity potential 𝛹 is directly related to the density distribution 𝜌 via the Poisson equation ∇ 𝛹 = 𝜋𝐺 𝜌 , (1)where 𝐺 is the gravitational constant. The general solution is 𝛹 ( 𝑟, 𝜃 ) = 𝜋𝐺 ∫ 𝑟 ∫ 𝜃 𝜌 ( ˜ 𝑟, ˜ 𝜃 )| 𝒓 − ˜ 𝒓 | ˜ 𝑟 sin ˜ 𝜃 𝑑 ˜ 𝑟 𝑑 ˜ 𝜃 (2)and the associated gravity force is 𝒈 = − ∇ 𝛹 . It is useful to separatethe external gravity potential into a spherically symmetric part ofa non-rotating, solid body and a series of higher order terms origi-nating from density perturbations and the rotational deformation ofthe planet: 𝛹 = − 𝐺 𝑀𝑟 (cid:34) − ∞ ∑︁ ℓ = 𝐽 ℓ (cid:18) 𝑅𝑟 (cid:19) ℓ 𝑃 ℓ ( cos 𝜃 ) (cid:35) , (3)where 𝑅 is the equatorial planetary radius, 𝑀 the total mass and 𝑃 ℓ are the Legendre polynomials of degree ℓ . The ℓ = 𝐽 ℓ MNRAS , 1– ???? (202-) onal Winds and Gravity are given by: 𝐽 ℓ = − 𝜋𝑀 𝑅 ℓ ∫ 𝑅 ∫ 𝜋 𝜌 ( 𝑟, 𝜃 ) 𝑃 ℓ ( cos 𝜃 ) 𝑟 ℓ + sin 𝜃 𝑑𝜃 𝑑𝑟 . (4)While the odd moments contain the signal of the equatorially anti-symmetric component of zonal flows, the even moments are dom-inated by the effects of the rotational deformation of the planet.We therefore focus on the equatorially antisymmetric winds and therespective observed odd gravity moments 𝐽 , 𝐽 , 𝐽 and 𝐽 . Thegeneral relation between the zonal mass flux and a density anomalyis expressed by the reduced Navier-Stokes equation, that describesthe conservation of momentum for a steady, inviscid, non-magnetic,inertia-less flow rotating around the ˆ 𝒛 -axis with a rotation rate Ω .In the co-rotating frame of reference this reads:2 𝛀 × ( 𝜌 𝒖 ) = − ∇ 𝑝 + 𝜌 ∇ 𝛹 + 𝜌 𝛀 × 𝛀 × 𝒓 , (5)where the terms (from left to right) are the Coriolis force, the pres-sure gradient, the gravity and the centrifugal force. This leadingorder force balance applies to the quasi stationary zonal flows in theouter envelope where the electrical conductivity is so low that theLorentz force can be neglected. The centrifugal force in the leadingorder force balance represents the rotational deformation and ren-ders the steady background state two-dimensional and non-spherical(e.g. Cao & Stevenson 2017). However, the rotational deformationitself is rather insignificant for the antisymmetric problem (Konget al. 2016) and thus we ignore the centrifugal forces for now.Pressure, density and gravity are separated into a hydrostaticbackground that depends only on radius and a small perturbation,e.g. 𝜌 = 𝜌 ( 𝑟 ) + 𝜌 (cid:48) ( 𝑟, 𝜃 ) . The first order perturbation equation is thengiven by:2 𝛀 × ( 𝜌 𝒖 ) = − ∇ 𝑝 (cid:48) + 𝜌 (cid:48) ∇ 𝛹 + 𝜌 ∇ 𝛹 (cid:48) , (6)where the last two terms on the right hand side are gravity forcecontributions due to a dynamic (i.e. flow-induced) density anomaly( 𝜌 (cid:48) ∇ 𝛹 ) and the dynamic self-gravity (DSG, 𝜌 ∇ 𝛹 (cid:48) ) term. In theclassic thermal wind approach (e.g. Kaspi et al. 2018), the DSG isneglected. The density anomaly can then simply be found from thethermal wind equation (TWE), which is the azimuthal componentof the curl of eq. 6:2 Ω 𝜕 𝑧 (cid:0) 𝜌𝑈 𝜙 (cid:1) = − 𝑟 ∇ 𝛹 𝜕 𝜃 𝜌 (cid:48) . (7)An integration along latitude and division by the background gravityyields the anomalous density field, which is subsequently used tocalculate 𝐽 ℓ by eq. 4.Zhang et al. (2015) and Wicht et al. (2020) show that the DSGterm represents a first order effect and, for example, changes the 𝐽 -values by up to 30%. Keeping the DSG term in the azimuthalvorticity equation yields the so called thermo-gravitational windequation (TGWE):2 Ω 𝜕 𝑧 (cid:0) 𝜌𝑈 𝜙 (cid:1) = ˆ 𝝓 · (cid:16) ∇ × (cid:16) 𝜌 (cid:48) ∇ 𝛹 + 𝜌 ∇ 𝛹 (cid:48) (cid:17)(cid:17) (8) = − 𝑟 𝑑𝛹𝑑𝑟 𝜕 𝜃 𝜌 (cid:48) − 𝑟 𝑑𝜌𝑑𝑟 𝜕 𝜃 𝛹 (cid:48) . (9)Replacing 𝛹 (cid:48) by eq. 2 and integrating over latitude leads to theintegro-differential equation solved by Zhang et al. (2015)2 Ω ∫ 𝜃 𝜕 𝑧 (cid:0) 𝜌𝑈 𝜙 ( 𝑟, ˜ 𝜃 ) (cid:1) 𝑑 ˜ 𝜃 = − 𝑟 𝑑𝛹𝑑𝑟 𝜌 (cid:48) − 𝑟 𝑑𝜌𝑑𝑟 ∫ 𝑟 ∫ 𝜃 ˜ 𝑟 𝜌 (cid:48) | 𝒓 − ˜ 𝒓 | sin ˜ 𝜃 𝑑 ˜ 𝜃𝑑 ˜ 𝑟 + 𝐶 ( 𝑟 ) . (10) While the integration function 𝐶 ( 𝑟 ) renders 𝜌 (cid:48) mathematically non-unique, the gravity moments nevertheless remain unique (Kaspiet al. 2010; Zhang et al. 2015). The integration function 𝐶 ( 𝑟 ) wouldonly contribute to the spherical symmetric gravity harmonic whichis determined by the total mass and therefore we can set 𝐶 ( 𝑟 ) tozero. Treating eq. 7 is mathematically demanding and has so faronly been solved for a simple interior model, i.e. a polytrope ofindex unity (Zhang et al. 2015; Kong et al. 2018).A potential work-around was devised by Braginsky & Roberts(1995) in the framework of classic geodynamo theory. It was shownthat the DSG term can be absorbed into the so-called effective variables (density and pressure): 𝑝 𝑒 = 𝑝 (cid:48) + 𝜌𝛹 (cid:48) (11) 𝜌 𝑒 = 𝜌 (cid:48) + 𝜇 𝜋𝐺𝛹 (cid:48) . (12)The radial function 𝜇 is thereby characterised by the compressibilityof the considered medium: 𝜇 𝜋𝐺 = 𝜌 𝜕𝜌𝜕 𝑝 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝑠 = 𝜌𝑐 𝑠 = 𝑔 𝑑𝜌𝑑𝑟 , (13)where 𝑐 𝑠 is the sound speed. For a polytropic perfect gas this canbe further simplified, 𝜇 𝜋𝐺 = 𝑚𝑚 + 𝜌 𝑚 + 𝑚 𝑐 𝑝 𝑐 𝜌 𝑚 − 𝑚 , (14)where 𝑚 is the polytropic index, 𝜌 𝑐 and 𝑝 𝑐 are density and pressureat the centre of the planet. They depend on the specific solution ofthe Lane-Emden equation. Moreover, for a polytropic index unity, 𝜇 is a constant and amounts to 𝜋 / 𝑅 .The effective variables, 𝑝 𝑒 and 𝜌 𝑒 , are equal to 𝑝 (cid:48) and 𝜌 (cid:48) reduced by the contribution of the local elevation or depression ofthe associated equipotential surface thus capturing the effect of theDSG. Using the effective variables and the definition of 𝜇 (eq. 13),the Navier-Stokes equation (eq. 6) simplifies to:2 𝛀 × ( 𝜌 𝒖 ) = − ∇ 𝑝 𝑒 + 𝜌 𝑒 ∇ 𝛹 . (15)Taking the azimuthal component of the curl, integrating along co-latitude 𝜃 and replacing ∇ 𝛹 = 𝑔 , then yields the effective densityperturbation: 𝜌 𝑒 ( 𝑟, 𝜃 ) = Ω 𝑟𝑔 ∫ 𝜃 𝜕 𝑧 (cid:0) 𝜌𝑈 𝜙 (cid:1) 𝑑 ˜ 𝜃 (cid:66) 𝜌 𝑈 . (16)Since this effective density disturbance formulates the zonal flowimpact, it has been called 𝜌 𝑈 by Wicht et al. (2020) and we adoptthis name here. Note that this is not the true density disturbance 𝜌 (cid:48) that could serve to calculate 𝐽 ℓ via eq. 4. In particular, 𝜌 𝑈 and 𝜌 (cid:48) are related as defined by eq. 12.Now we can find an equation for the gravity potential by re-placing 𝜌 (cid:48) with 𝛹 (cid:48) using eq. 1 in eq. 9, integrating along latitudeand making use of eq. 16: ∇ 𝛹 (cid:48) + 𝜇𝛹 (cid:48) = 𝜋𝐺 𝜌 𝑈 . (17)This is a two-dimensional, inhomogeneous PDE of secondorder, which describes the wind-induced anomalies in the gravitypotential of a gas planet (Wicht et al. 2020). The effective densityperturbation 𝜌 𝑈 derived from eq. 16 acts as the source term for thisHelmholtz-like equation. Only when 𝜇 = 𝑐𝑜𝑛𝑠𝑡 does this equationbecome an inhomogeneous Helmholtz equation and can be solved ina semi-analytical way (Wicht et al. 2020). In the more general casewhere 𝜇 depends on the radius we refer to the numerical methodsdiscussed in sec. 4. MNRAS , 1– ?? (202-) W. Dietrich et al.
It is important to note that the 2nd order differential equation(eq. 17) and the integro-differential form of the TGWE (eq. 10)describe the same physical problem. The main difference is thateq. 17 solves for 𝛹 (cid:48) , while eq. 10 solves for the density anomaly 𝜌 (cid:48) . Eq. 17 not only directly provides the gravity potential we areinterested in, but is also much easier to solve. The effective variables can be further exploited to show that theflow ( 𝑈 𝜙 ) and not the mass flux ( 𝜌𝑈 𝜙 ) is geostrophic and shouldbe initially extended along cylinders. To emphasise under whichconditions geostrophic winds can be modulated along the rotationaxis, we express the density as a function of pressure and entropy: 𝑠 (cid:48) 𝑐 𝑝 = 𝑐 𝑣 𝑐 𝑝 𝑝 𝑒 𝑝 − 𝜌 𝑒 𝜌 = 𝑐 𝑣 𝑐 𝑝 𝑝 (cid:48) 𝑝 − 𝜌 (cid:48) 𝜌 . (18)Dividing the Navier-Stokes equation (eq. 15) by the backgrounddensity and defining a reduced pressure 𝑝 ★ = 𝑝 𝑒 / 𝜌 , eq. 15 yields:2 𝛀 × 𝒖 = − 𝜌 ∇ 𝑝 𝑒 + 𝜌 𝑒 𝜌 ∇ 𝛹 = − ∇ (cid:18) 𝑝 𝑒 𝜌 (cid:19) − 𝑝 𝑒 𝜌 ∇ 𝜌 + (cid:18) 𝑐 𝑉 𝑐 𝑝 𝑝 𝑒 𝑝 − 𝑠 (cid:48) 𝑐 𝑝 (cid:19) ∇ 𝛹 = − ∇ 𝑝 ★ − 𝑠 (cid:48) 𝑐 𝑝 ∇ 𝛹 + (cid:32) 𝑐 𝑉 𝑐 𝑝 𝜌 𝑝 ∇ 𝛹 − ∇ 𝜌 (cid:33) 𝑝 𝑒 𝜌 (cid:39) − ∇ 𝑝 ★ − 𝑠 (cid:48) 𝑐 𝑝 ∇ 𝛹 . (19)Note, the term in the brackets scales with the non-adiabaticity of thebackground state and hence can safely be neglected for a vigorouslyconvecting atmosphere like Jupiter’s. This equation highlights thatthe winds 𝑈 𝜙 and not the mass fluxes are cylindrically invariant.Only if there are sizeable latitudinal gradients in the zonally aver-aged entropy, are deviations from geostrophy possible, if we restrictthe consideration to regions of negligible electrical conductivity.The associated temperature fluctuations required to drive the zonalwind out of its cylindrical invariance can be estimated by azimuthalcomponent of the curl of eq. 19:2 Ω 𝜕 𝑧 𝑈 𝜙 = − 𝑔𝑐 𝑝 𝑟 𝜕 𝜃 𝑠 (cid:48) , (20)where 𝑑𝛹 / 𝑑𝑟 = 𝑔 has been used. Assuming that the vertical andlatitudinal derivative can be approximated with the same lengthscale (e.g. close to the equator), thus 𝜕 𝑧 𝑈 𝜙 ≈ Δ 𝑈 / 𝛿 and 1 / 𝑟 𝜕 𝜃 𝑠 (cid:48) ≈ 𝛿𝑠 / 𝛿 . Furthermore, the entropy fluctuations can be (to first order)approximated by temperature fluctuations, thus 𝛿𝑠 ≈ 𝑐 𝑝 / 𝑇 𝛿𝑇 . Thisthen yields 𝛿𝑇 ≈ Ω 𝑔 𝑇 Δ 𝑈 . (21)For Jupiter Ω ≈ . · − s − and 𝑔 ≈
25 m / s . Then, in orderto induce a vertical variation of the wind Δ 𝑈 =
10 m / s at a tem-perature of 𝑇 ( 𝑟 = . 𝑅 ) = · K (Nettelmann et al. 2012),we find 𝛿𝑇 ≈ . − K (e.g. Jones 2007). These estimatesare applicable to the convective part of the atmosphere, i.e. below 𝑝 = 𝑈 𝜙 (and not the mass flux) arecylindrically invariant. The dynamic density perturbation (eq. 16) is governed by the as-sumed interior structure of the planet via the background densityand gravity, 𝜌 and 𝑔 , respectively, and the 𝑧 -gradient of the zonalmass flux. Theoretical considerations and numerical simulationssuggest a geostrophic zonal flow structure for the fast rotation andlow viscosity gas planets (e.g. Taylor 1917; Dietrich & Jones 2018;Gastine & Wicht 2021). This means that the flow depends only onthe distance 𝑠 = 𝑟 sin 𝜃 to the rotation axis, where 𝜃 is the colatitude.We could then simply downward continue the surface zonal flow 𝑈 along cylinders.However, the gravity measurements and the secular variationof the Jupiter’s magnetic field show that the wind speed must signif-icantly decrease with depth and should be confined to outer 5% ofthe planetary radius (Kaspi et al. 2018; Moore et al. 2019; Galanti& Kaspi 2020). This decrease is commonly parametrized with anadditional radial decay function 𝑄 𝑟 , 𝑈 𝜙 ( 𝑟, 𝜃 ) = 𝑄 𝑟 ( 𝑟 ) 𝑈 𝑜 ( 𝑠 ) (22)with 𝑄 𝑟 ( 𝑅 ) =
1. The cause for the deviation from geostrophy re-mains unspecified. Electromagnetic effects have been alluded to.Buoyancy forces arising in a stably stratified region with the assis-tance of Lorentz forces are a promising mechanism (Christensenet al. 2020; Gastine & Wicht 2021). We start with discussing theinterior state, then the different models of the surface zonal profileand finally the calculation of the dynamic density source 𝜌 𝑈 . The source term (eq. 26) and the DSG coefficient, 𝜇 , depend onthe background density and gravity profile, 𝜌 and 𝑔 , respectively.We explore their influence via the chosen interior state model bycomparing three commonly used models.A simple model for Jupiter’s interior is a polytropic ideal gasof index unity that provides a decent description of the pressuredependence on the density (Hubbard 1999), but not of the ( 𝑝 − 𝑇 )-curve. The background density is then: 𝜌 ( 𝑟 ) = 𝜌 𝑐 sin 𝜒𝜒 , with 𝜒 = 𝜋𝑟𝑅 , (23)where 𝜌 𝑐 is the central density. The associated background gravityis given by 𝑔 ( 𝑟 ) = − 𝑅𝐺 𝜌 𝑐 𝜒 cos 𝜒 − sin 𝜒𝜒 (24)and is hence directly proportional to the radial density gradient asdiscussed in Zhang et al. (2015). This proportionality ( 𝑔 ∝ 𝑑𝜌 / 𝑑𝑟 )is an exclusive property of a polytrope with index unity and implies aconstant 𝜇 = 𝜋 / 𝑅 . This has been exploited by Zhang et al. (2015)and Wicht et al. (2020) to solve the TGWE. The grey profiles infig. 2 illustrate the density, its radial gradient and 𝜇 for the polytropeof index 𝑚 = 𝜇 for a setup with interpolated hydrogen EOS, a 1-bar temperature of165 K, a core of 4.2 earth masses and heavy element abundance of 33 MNRAS , 1– ????
1. The cause for the deviation from geostrophy re-mains unspecified. Electromagnetic effects have been alluded to.Buoyancy forces arising in a stably stratified region with the assis-tance of Lorentz forces are a promising mechanism (Christensenet al. 2020; Gastine & Wicht 2021). We start with discussing theinterior state, then the different models of the surface zonal profileand finally the calculation of the dynamic density source 𝜌 𝑈 . The source term (eq. 26) and the DSG coefficient, 𝜇 , depend onthe background density and gravity profile, 𝜌 and 𝑔 , respectively.We explore their influence via the chosen interior state model bycomparing three commonly used models.A simple model for Jupiter’s interior is a polytropic ideal gasof index unity that provides a decent description of the pressuredependence on the density (Hubbard 1999), but not of the ( 𝑝 − 𝑇 )-curve. The background density is then: 𝜌 ( 𝑟 ) = 𝜌 𝑐 sin 𝜒𝜒 , with 𝜒 = 𝜋𝑟𝑅 , (23)where 𝜌 𝑐 is the central density. The associated background gravityis given by 𝑔 ( 𝑟 ) = − 𝑅𝐺 𝜌 𝑐 𝜒 cos 𝜒 − sin 𝜒𝜒 (24)and is hence directly proportional to the radial density gradient asdiscussed in Zhang et al. (2015). This proportionality ( 𝑔 ∝ 𝑑𝜌 / 𝑑𝑟 )is an exclusive property of a polytrope with index unity and implies aconstant 𝜇 = 𝜋 / 𝑅 . This has been exploited by Zhang et al. (2015)and Wicht et al. (2020) to solve the TGWE. The grey profiles infig. 2 illustrate the density, its radial gradient and 𝜇 for the polytropeof index 𝑚 = 𝜇 for a setup with interpolated hydrogen EOS, a 1-bar temperature of165 K, a core of 4.2 earth masses and heavy element abundance of 33 MNRAS , 1– ???? (202-) onal Winds and Gravity r [ k g / m ] polytropeN12G03 0.85 0.90 0.95 1.00 r [ k g / m ] polytropeN12G03 0.85 0.90 0.95 1.00 r R polytropeN12G03 Figure 2.
Jupiter’s interior density 𝜌 (left), radial gradient ∇ 𝜌 and 𝜇 = ∇ 𝜌 / 𝑔 in the outer 15% of radius from two Jupiter interior models in comparison to apolytropic perfect gas with index unity. The interior models based on 3-layer internal structure are adopted from Guillot et al. (2003) (blue) and Nettelmannet al. (2012) (red). The dots on the curve are based in ab-initio calculations of the sound speed ( 𝑐 𝑠 ∝ 𝜇 ) by French et al. (2012). earth masses (Guillot et al. 2003). In comparison to the rather simplepolytropic model (grey), the densities are substantially higher for 𝑟 < . 𝑅 and slightly lower above this radius. Thus the densitygradient shows a pronounced maximum around 𝑟 = . 𝑅 (fig. 2,middle panel). The DSG coefficient ( 𝜇 ) increases with radius andreaches a 2.5 times larger value than for the polytrope.Alternatively, we use the more recent calculations from Nettel-mann et al. (2012) based on the updated H-REOS2 model (‘N12‘).The DSG coefficient is related to the sound speed (eq. 13), whichwas calculated for the same Jupiter model by French et al. (2012).For this model (termed J11-8a) the depth of the molecular-metallicphase transition is at 𝑝 = 𝑟 < . 𝑅 and falls betweenthe G03 model and the polytrope at larger radii. The density gradi-ent and hence 𝜇 show two distinct maxima with smaller amplitudethan the G03 model. The surface flow of Jupiter is deduced from tracking cloud features,either with the Hubble Space Telescope (HST) or from space crafts.Fig. 3 displays several zonal flow measurements that have beenobtained over the last decades. The oldest illustrated flows (blueand yellow) are based on observations by Voyager I and II (Limayeet al. 1982; Limaye 1986) and represent the flow in 1979. HSTmonitored the wind structure several times between 1995 and 1998(green profile, ? ). Cassini measurements of the wind profile stemfrom late 2000 (red, Porco et al. 2003). Later, various differentHST-based profiles were obtained between 2009 and 206, duringthe perijove nine of the Juno space craft (Tollefson et al. 2017).Kaspi et al. (2018) and Galanti & Kaspi (2020) use this most recentflow profile for their gravity data analysis, while Kong et al. (2018)prefer the model based on Cassini images from late 2000 to early2001 presented in Porco et al. (2003).All flow models show the same principle structure but alsodiffer in some details. The equatorially antisymmetric flow contri-bution shown in fig. 3 is clearly dominated by the prograde jetwhich starts at a latitude around 𝛽 = ◦ and extends to about 𝛽 = ◦ . The amplitude of this jet varies between 65 m/s and 75 m/sfor the different models.Tollefson et al. (2017) conclude that the variations exceed themodel errors, at least for some epochs and at some latitudes. TheirLomb-Scargle periodogram analysis reveals dominant variation pe-riods of about 7 yr and 14 yr. Since the latter is close enough to Jupiter’s orbital period of 11 . An obvious problem arises with the equatorially antisymmetric con-tributions to the dynamic density source. Geostrophically downwardcontinuing the surface flow in each hemisphere separately yields adiscontinuity, or step, at the equatorial plane. It can be large in casethe geostrophic continuation of the surface flow hits the equato-rial plane at a radius where 𝑄 𝑟 is significantly larger than zero.This equatorial step seems unphysical but can easily be dealt withmathematically.Using the product ansatz (eq. 22) in eq. 16, we can separatethe dynamic density perturbation into two contributions: 𝜌 𝑈 = − Ω 𝑟𝑔 (cid:20) 𝜕𝜕𝑟 ( 𝜌𝑄 𝑟 ) ∫ 𝜃 𝑈 𝑜 cos ˆ 𝜃𝑑 ˆ 𝜃 + 𝜌𝑄 𝑟 ∫ 𝜃 𝜕𝑈 𝑜 𝜕𝑧 𝑑 ˆ 𝜃 (cid:21) . (25)The second integral contributes only at the equator where the 𝑧 -derivative yields a delta-peak. This can be integrated analytically: 𝜌 𝑈 = − Ω 𝑟𝑔 (cid:20) 𝜕𝜕𝑟 ( 𝜌𝑄 𝑟 ) ∫ 𝜃 𝑈 𝑜 cos ˆ 𝜃𝑑 ˆ 𝜃 + 𝜌𝑄 𝑟 𝑟 𝐻 ( 𝜃 − 𝜋 / ) Δ 𝑈 (cid:21) , (26)where 𝐻 is the Heaviside step function and Δ 𝑈 = 𝑈 𝑜 ( 𝜃 → 𝜋 / ) .As we will show below, the equatorial step can yield an impor-tant contribution to the gravity signal which is problematic. Kaspiet al. (2018) and Galanti & Kaspi (2020) therefore ignore the respec-tive term in an approach they call ‘hemispheric‘. They calculate 𝜌 𝑈 in each hemisphere but dismiss the contribution at (or very closeto) the equator, which is equivalent to evaluating 𝜌 𝑈 = − Ω 𝑟𝑔 𝜕𝜕𝑟 ( 𝜌𝑄 𝑟 ) ∫ 𝜃 𝑈 𝑜 cos ˆ 𝜃𝑑 ˆ 𝜃 . (27) MNRAS , 1– ?? (202-) W. Dietrich et al.
Figure 3.
Comparison of different zonal flow models for Jupiter observed either by HST or in-situ by space crafts. Shown is the equatorial antisymmetric flowpart relevant for the odd gravity moments. The models are from the Voyager mission Limaye et al. (1982); Limaye (1986), Cassini Porco et al. (2003), andvarious Hubble Space Telescope (HST) campaigns ( ? Tollefson et al. 2017).
Figure 4.
Zonal flow from the full 3D numerical model (a) separated into theequatorial symmetric part (b) and the antisymmetric (c) part. The contourlevels are identical amongst the panels.
This approach seems to offer a simple fix but is also inconsistentsince the step is an integral part of the chosen flow model.Alternatively, Kong et al. (2018) use an additional linear 𝑧 -dependence of 𝑈 𝑜 , such that the vertical gradient at the equator isregular and smooth. Their zonal flow model is thus given by 𝑈 ( 𝑟, 𝜃 ) = 𝑄 𝑟 | 𝑧 || 𝑧 𝑜 | 𝑈 𝑜 ( 𝑠 ) , (28)where | 𝑧 𝑜 | = ( 𝑅 − 𝑠 ) / is the local distance from the equatorialplane to the surface. The equatorial step and the linear 𝑧 -dependenceare the steepest and smoothest end member, respectively, of thepossible functions that reconcile the northern and southern zonalflows.Here we introduce a novel approach guided by the physi-cal principles of atmospheric dynamics. The generation of deep-reaching zonal flows powered by the statistical correlations of theconvective flows (Reynolds stresses) are best captured in numeri-cal simulations (Heimpel et al. 2005; Christensen 2002; Dietrich& Jones 2018; Gastine & Wicht 2021). Fig. 4 illustrates the zonalflow in a typical simulation of convection in a fast rotating spheri-cal shell with an aspect ratio of 𝑟 𝑖 / 𝑅 = . 𝐸 = · − , 𝑃𝑟 = . , 𝑅𝑎 = · , 𝑁 𝜌 = . 𝑄 𝑠 that reduces the amplitudeoutside of the latitude 𝛽 TC , i.e. where the TC touches the outerboundary: 𝑈 ★𝑜 = 𝑈 𝑜 ( 𝑠 ) 𝑄 𝑠 , (29)with 𝑄 𝑠 = (cid:20) tanh (cid:18) 𝛽 − 𝛽 TC 𝛿 𝛽 (cid:19) − tanh (cid:18) 𝛽 + 𝛽 TC 𝛿 𝛽 (cid:19)(cid:21) + . (30)A second parameter which is introduced here is the width of thelatitudinal cut-off, 𝛿 𝛽 . Both functions, 𝑄 𝑠 ( 𝛽 TC , 𝛿 𝛽 ) and 𝑄 𝑟 ( ℎ, 𝛿 ℎ ) ,define an individual tangent cylinder and hence should be consistent.More details on selecting the best parameter combination for 𝛽 TC , 𝛿 𝛽 , ℎ and 𝛿 ℎ are discussed in sec. 5. When 𝛽 TC is sufficientlylarge and the radial decay function drops rapidly below a certaindepth, the resulting flow splits into an independent northern andsouthern part. Then the equatorial step contribution is identical tozero. The underlying assumption is, that part of the observed zonalflow at cloud level, in particular its antisymmetric parts at lowlatitude, are shallow. We term cases that apply such an attenuationof antisymmetric surface flows at low latitudes ’TC-models’. MNRAS , 1– ????
This approach seems to offer a simple fix but is also inconsistentsince the step is an integral part of the chosen flow model.Alternatively, Kong et al. (2018) use an additional linear 𝑧 -dependence of 𝑈 𝑜 , such that the vertical gradient at the equator isregular and smooth. Their zonal flow model is thus given by 𝑈 ( 𝑟, 𝜃 ) = 𝑄 𝑟 | 𝑧 || 𝑧 𝑜 | 𝑈 𝑜 ( 𝑠 ) , (28)where | 𝑧 𝑜 | = ( 𝑅 − 𝑠 ) / is the local distance from the equatorialplane to the surface. The equatorial step and the linear 𝑧 -dependenceare the steepest and smoothest end member, respectively, of thepossible functions that reconcile the northern and southern zonalflows.Here we introduce a novel approach guided by the physi-cal principles of atmospheric dynamics. The generation of deep-reaching zonal flows powered by the statistical correlations of theconvective flows (Reynolds stresses) are best captured in numeri-cal simulations (Heimpel et al. 2005; Christensen 2002; Dietrich& Jones 2018; Gastine & Wicht 2021). Fig. 4 illustrates the zonalflow in a typical simulation of convection in a fast rotating spheri-cal shell with an aspect ratio of 𝑟 𝑖 / 𝑅 = . 𝐸 = · − , 𝑃𝑟 = . , 𝑅𝑎 = · , 𝑁 𝜌 = . 𝑄 𝑠 that reduces the amplitudeoutside of the latitude 𝛽 TC , i.e. where the TC touches the outerboundary: 𝑈 ★𝑜 = 𝑈 𝑜 ( 𝑠 ) 𝑄 𝑠 , (29)with 𝑄 𝑠 = (cid:20) tanh (cid:18) 𝛽 − 𝛽 TC 𝛿 𝛽 (cid:19) − tanh (cid:18) 𝛽 + 𝛽 TC 𝛿 𝛽 (cid:19)(cid:21) + . (30)A second parameter which is introduced here is the width of thelatitudinal cut-off, 𝛿 𝛽 . Both functions, 𝑄 𝑠 ( 𝛽 TC , 𝛿 𝛽 ) and 𝑄 𝑟 ( ℎ, 𝛿 ℎ ) ,define an individual tangent cylinder and hence should be consistent.More details on selecting the best parameter combination for 𝛽 TC , 𝛿 𝛽 , ℎ and 𝛿 ℎ are discussed in sec. 5. When 𝛽 TC is sufficientlylarge and the radial decay function drops rapidly below a certaindepth, the resulting flow splits into an independent northern andsouthern part. Then the equatorial step contribution is identical tozero. The underlying assumption is, that part of the observed zonalflow at cloud level, in particular its antisymmetric parts at lowlatitude, are shallow. We term cases that apply such an attenuationof antisymmetric surface flows at low latitudes ’TC-models’. MNRAS , 1– ???? (202-) onal Winds and Gravity The different forms of the radial decay function 𝑄 ( 𝑟 ) that have beensuggested to explain the gravity observations are illustrated in Fig. 1.Kaspi et al. (2018) use a combination of an exponential decay anda hyperbolic tangent: 𝑄 𝑟 ( 𝛼, ℎ, 𝛿 ℎ ) = 𝛼 (cid:20) tanh (cid:18) 𝑟 − ( 𝑅 − ℎ ) 𝛿 ℎ (cid:19) + (cid:21) / (cid:20) tanh (cid:18) ℎ𝛿 ℎ (cid:19) + (cid:21) + ( − 𝛼 ) exp (cid:18) 𝑟 − 𝑅ℎ (cid:19) . (31)Here 𝛼 is the weight of the hyperbolic tangent contribution, ℎ thedecay depth and 𝛿 ℎ the decay width. For explaining the gravity ob-servations, Kaspi et al. (2018) propose a large 𝛼 = .
92, a relativelylarge 𝛿 ℎ = ℎ = 𝑄 𝑟 ( ℎ, 𝐻 ) = exp (cid:20) ℎ (cid:18) − 𝐻 𝐻 − ( − 𝑟 ) (cid:19)(cid:21) , (32)for 𝑟 > 𝐻 , and 𝑄 𝑟 = . 𝑟 (cid:54) 𝐻 . Kong et al. (2018) report that theobservations are best matched with a combination of 𝐻 =
10 484 kmand 𝛿 ℎ =
15 377 km. Fig. 1 shows that the respective profile (orange)decays also rather smoothly.Both the solution suggested by Kaspi et al. (2018) and by Konget al. (2018) are not compatible with the magnetic observations(Moore et al. 2019). Furthermore, vigorous convection provides analmost adiabatic environment and hence the decay functions shouldbe rather flat with a sharp drop at greater depth (see also sec. 2.2).Realizing this, Galanti & Kaspi (2020) proposed a steeper decayingalternative with a hyperbolic tangent outer, and an exponential innerbranch: 𝑄 𝑟 = 𝛼 tanh (cid:16) 𝑟 − 𝑟 𝑇 𝛿 ℎ (cid:17) / tanh (cid:16) ℎ𝛿 ℎ (cid:17) + − 𝛼 for 𝑟 (cid:62) 𝑟 𝑇 ( − 𝛼 ) exp (cid:16) 𝑟 − 𝑟 𝑇 𝛿 (cid:48) ℎ (cid:17) for 𝑟 < 𝑟 𝑇 . (33)They suggest 𝛼 = . ℎ = 𝑟 𝑇 / 𝑅 = .
972 and 𝛿 ℎ = 𝛿 (cid:48) ℎ =
204 km in conjunction with an optimised surface flowprofile to explain the gravity observations. The radial profile dropsalmost faster than the magnetic constraints require (see fig.1).We adopt the pure hyperbolic tangent profile: 𝑄 𝑟 ( ℎ, 𝛿 ℎ ) = (cid:20) tanh (cid:18) 𝑟 − ( 𝑅 − ℎ ) 𝛿 ℎ (cid:19) + (cid:21) / (cid:20) tanh (cid:18) ℎ𝛿 ℎ (cid:19) + (cid:21) . (34)For 𝛿 ℎ =
200 km and ℎ = ℎ = 𝛿 ℎ =
500 km. Fig. 1 shows the respective decay function inblue.
The semi-analytical method for solving the TGWE (eq. 17) devel-oped by Wicht et al. (2020) is restricted to the case of a constant 𝜇 .However, as shown in fig. 2, 𝜇 ( 𝑟 ) varies strongly and reaches muchhigher values than 𝜋 / 𝑅 , particularly in the outer 10% of Jupiter’sradius where the flow-induced gravity moments originate from. Wehave therefore developed a numerical method that can also handleradial variations in 𝜇 .We first formulate the dynamic density source 𝜌 𝑈 ( 𝑟, 𝜃 ) viaeq. 26 in spatial space by choosing a radial decay function 𝑄 𝑟 , an interior state model setting 𝜌 , 𝑔 and 𝜇 , and a surface zonal flowprofile 𝑈 𝑜 , that is cylindrically downward continued.We then expand the latitudinal dependence of the gravity po-tential perturbation 𝛹 (cid:48) and the dynamic density 𝜌 𝑈 in Legendrepolynomials, e.g. for the former this is given by 𝛹 (cid:48) ( 𝑟, 𝜃 ) = 𝐿 ∑︁ ℓ = 𝛹 (cid:48) ℓ ( 𝑟 ) 𝑃 ℓ ( cos 𝜃 ) , (35)with 𝛹 (cid:48) ℓ ( 𝑟 ) = ∫ 𝜋 𝛹 (cid:48) ( 𝑟, 𝜃 ) 𝑃 ℓ ( cos 𝜃 ) sin 𝜃𝑑𝜃 . (36)Since the background state, and hence 𝜇 , are only a function ofradius, the TGWE (eq. 17) decouples for each degree ℓ , and we arethen left with a set of radial ODEs, one for each spherical harmonicdegree ℓ : (cid:18) 𝜕 𝑟 + 𝑟 𝜕 𝑟 − ℓ ( ℓ + ) 𝑟 + 𝜇 (cid:19) 𝛹 (cid:48) ℓ ( 𝑟 ) = 𝜋𝐺 𝜌 𝑈ℓ ( 𝑟 ) . (37)At the outer boundary the potential must match to the solutionof a source free region, i.e. ∇ 𝛹 𝑝𝑜𝑡 =
0, with the characteristicradial dependence 𝑟 −( ℓ + ) . This leads to the mixed outer boundarycondition: 𝛹 (cid:48) ℓ ( 𝑅 ) = − 𝑅ℓ + 𝜕 𝑟 𝛹 (cid:48) ℓ ( 𝑅 ) , (38)and 𝛹 (cid:48) ℓ ( ) = ℓ we construct a set of N nor-malised orthogonal functions 𝑗 ★ℓ𝑛 ( 𝑘 ℓ𝑛 𝑟 ) . The radial scales, 𝑘 ℓ𝑛 ,are determined by finding the first N radii (starting at the origin), atwhich 𝑗 ★ℓ𝑛 ( 𝑘 ℓ𝑛 𝑟 ) fulfils the boundary condition (eq. 38). The radialfunctions 𝛹 (cid:48) ℓ ( 𝑟 ) are expanded in 𝑗 ★ℓ𝑛 : 𝛹 (cid:48) ℓ ( 𝑟 ) = 𝑁 ∑︁ 𝑛 = 𝛹 (cid:48) ℓ𝑛 𝑗 ★ℓ𝑛 ( 𝑘 ℓ𝑛 𝑟 ) , (39)where the expansion coefficients 𝛹 (cid:48) ℓ𝑛 are given by 𝛹 (cid:48) ℓ𝑛 = ∫ 𝑅 𝛹 (cid:48) ℓ ( 𝑟 ) 𝑗 ★ℓ𝑛 ( 𝑘 ℓ𝑛 𝑟 ) 𝑟 𝑑𝑟 . (40)Since the modified spherical Bessel functions are eigenfunctions ofthe Laplace operator, eq. 37 transforms into a set of linear algebraicequations for the expansion coefficients 𝛹 (cid:48) ℓ𝑛 : 𝑁 ∑︁ 𝑛 = (cid:16) 𝜇 ( 𝑟 ) − 𝑘 ℓ𝑛 (cid:17) 𝛹 (cid:48) ℓ𝑛 𝑗 ★ℓ𝑛 ( 𝑘𝑟 ) = 𝜋𝐺 𝜌 𝑈ℓ ( 𝑟 ) . (41)Eq. 41 defines a coupled set of algebraic equations and solved byusing a matrix formalism. The source term, 𝜌 𝑈ℓ , on the RHS isdiscretized along 𝑁 radial grid points, 𝑟 𝑖 . The solution of the matrixequation, 𝛹 (cid:48) ℓ𝑛 , contains the Bessel function expansion coefficientsof the gravity field perturbation. We thus introduce a square matrix H ℓ defined by H ℓ𝑛 𝑖 = (cid:16) 𝜇 ( 𝑟 𝑖 ) − 𝑘 ℓ𝑛 (cid:17) 𝑗 ★ℓ𝑛 ( 𝑘 ℓ𝑛 𝑟 𝑖 ) , (42)where 𝑛, 𝑖 ∈ [ , 𝑁 ] . Then the TGWE with radially varying DSGcoefficient can be written in the symbolic matrix form H ℓ 𝜳 (cid:48) ℓ𝑛 = 𝜋𝐺 𝝆 𝑈ℓ , (43) MNRAS , 1– ?? (202-) W. Dietrich et al. d eca y w i d t h δ h [ k m ] decay depth h [km] r a d i u s / d r op [ R ] Figure 5.
Trade-off between decay depth ℎ and decay width 𝛿 ℎ . Thecoloured contours indicate combinations, where individual gravity momentsare compatible within the uncertainties with the measurements ( 𝐽 blue, 𝐽 orange, 𝐽 green and 𝐽 red), while the grey shaded areas indicate a flowdrop of three orders of magnitude at various depths as required by Jupiter’ssecular variation (Moore et al. 2019). where 𝝆 𝑈ℓ represents the source vector per degree along all radialgrid points and 𝜳 (cid:48) ℓ𝑛 is the solution vector containing the expansioncoefficients. The matrix equation is solved via LU factorization. Fordealiasing, we thus ignore the highest 10% of the Bessel coefficients.We use an equidistant radial mesh with up to 𝑁 = 𝑁 𝜃 = ℓ = , , , 𝛹 (cid:48) ℓ𝑛 in Bessel space, eq. 39 yieldsthe radial representation 𝛹 (cid:48) ℓ . The gravity moments are then simplygiven by 𝐽 ℓ = 𝑅𝐺 𝑀 ℓ + 𝛹 (cid:48) ℓ ( 𝑅 ) , (44)If required, the density perturbation can be calculated by solving ∇ 𝛹 (cid:48) = 𝜋𝐺 𝜌 (cid:48) .To test our method, we compare the solution for constant 𝜇 withresults using the semi-analytical method introduced by Wicht et al.(2020). To verify our method with radius-dependent 𝜇 , we solveeq. 37 with an independent solver based on the shooting method.This is a standard tool for solving initial value problems and canreadily be applied here, making use of a variable transformation toaccount for the asymptotic behaviour at the origin when applyingthis method. We start by exploring the challenges of modelling the odd gravityperturbations for a reference model that combines the interior modelby Guillot et al. (2003) with the average flow introduced in sec. 3.2and employs the TGWE method (eq. 37 or 41). Furthermore thedynamic density source is calculated via eq. 26 and thus we assume ahemispherically geostrophic flow, keep the equatorial step and applythe hyperbolic tangent radial decay profile for 𝑄 𝑟 in accordance witheq. 34. Fig. 5 shows an attempt to model the observed gravity harmonics byvarying the two parameters in 𝑄 𝑟 , the decay depth ℎ and the decaywidth 𝛿 ℎ . A latitude-dependent attenuation in the form of eq. 28 oreq. 29 is not applied at this point. To assess the quality of the mod-elled gravity moments, we compare them individually for each de-gree ℓ to the observations rather than minimising an ℓ -independent,global cost function (Kaspi et al. 2018). Semi-transparent, coloured,regions indicate the parameter combinations where the respectivegravity harmonic stays within 3 𝜎 of the observations (Durante et al.2020). There is a trade-off between both parameters since increas-ing either boosts the impact of deeper regions. Unfortunately, thereis no parameter combination where all stripes overlap, i.e. no com-bination of ℎ and 𝛿 ℎ that would explain all four odd harmonicswith one single decay function. Gravity harmonic 𝐽 seems to beparticularly problematic as it always maintains some distance fromthe other three, since to match it requires particularly deep sources.The grey background contours show at which radius the radialdecay function 𝑄 𝑟 would have decayed by three orders of magni-tude. The SV constraint by Moore et al. (2019) suggests that thisshould happen somewhere between 𝑟 = . 𝑅 and 𝑟 = . 𝑅 .Modelling the observed 𝐽 always requires deeper sources that areincompatible with this constraint. For the smooth decay functionssuggested by Kaspi et al. (2018) or Kong et al. (2018) the threeorders of magnitude drop lies far below 𝑟 = . 𝑅 (see also fig. 1).Numerical simulations and theoretical consideration suggestthat the flow remains geostrophic until a stably stratified layer,Lorentz forces or a combination of both, drastically quench theamplitude over a few hundred kilometres (Christensen et al. 2020;Wicht & Gastine 2020). We therefore restrict our analysis to 𝛿 ℎ =
500 km in the following, a value that represents the scaleheight of the electrical conductivity in the outer atmosphere ofJupiter quite well (French et al. 2012; Wicht et al. 2019) and willbe further validated from the results of the TC model (see sec. 5.2).Fig. 6 illustrates how the gravity harmonics change when vary-ing ℎ while keeping 𝛿 ℎ =
500 km fixed. For the reference model(Fig. ?? , a), the harmonics 𝐽 , 𝐽 , and 𝐽 all agree with the re-spective observations for ℎ values between 1500 km and 2300 km.However, 𝐽 requires much deeper flows with ℎ ≈ ℎ and the respective uncertainties for all consid-ered models. We use the spherical harmonic degree ℓ as an indexto denote the different values of ℎ ℓ required to explain the differentobservations. For a successful model, all four ℎ -values must agreewithin the errors. This is clearly not the case for our reference model.Only for degree ℓ = ℓ = ℎ ≈ ℎ between2000 km and 4000 km the modelled gravity moment 𝐽 is of thewrong sign, whereas 𝐽 is much larger than the observed value inthat range of ℎ .Panel b) in fig. 6 illustrates that the situation improves whenwe follow the approach by Kaspi et al. (2018) and ignore the contri-bution due to the equatorial step in eq. 26. While the values of ℎ , ℎ and ℎ remain nearly unchanged, ℎ decreases by half to about3200 km (see also tab. 1). This shows the large impact of the equa-torial discontinuity, in particular on degree ℓ =
3. If in addition tothis we assume the smoothly decaying radial function 𝑄 𝑟 suggestedby Kaspi et al. (2018), ignore the DSG term and utilise the perijovesurface flow model, we can largely reproduce their results and find areasonably good agreement of the decay depths across the degrees.Panel c) of fig. 6) shows the results when using the additionallinear z-dependence suggested by Kong et al. (2018). This leads to MNRAS , 1– ????
3. If in addition tothis we assume the smoothly decaying radial function 𝑄 𝑟 suggestedby Kaspi et al. (2018), ignore the DSG term and utilise the perijovesurface flow model, we can largely reproduce their results and find areasonably good agreement of the decay depths across the degrees.Panel c) of fig. 6) shows the results when using the additionallinear z-dependence suggested by Kong et al. (2018). This leads to MNRAS , 1– ???? (202-) onal Winds and Gravity Figure 6.
Odd gravity moments (degree ℓ = , , , 𝜎 uncertainties (Durante et al. 2020). The panel a) represents the reference model using the full surface flow, a geostrophic downward continuationand inlcuding the equatorial step, whereas in b) the equatorial step is ignored. In c) the additional linear 𝑧 -decay is applied, and the d) utilises a zero flowoutside TC. More details in the text or in tab. 1. shallower winds closer to the equator, a smooth vertical gradient andthus no equatorial step. Again, ℎ decreases significantly, but notas much as for the model that ignores the equatorial step. The otherharmonics are somewhat more affected, such that all ℎ ℓ increase byroughly 10% with respect to the reference model (tab. 1). Thoughthis additional modification of the vertical flow structure avoidsthe equatorial discontinuity, it is hard to justify physically. Againwe would have to also adopt the other model ingredients (TGWEmodel with polytropic interior, Cassini flow with capped spike,inverse Gauss decay profile) to reproduce the results by Kong et al.(2018) and match all harmonics. Finally we explore our TC model where the flow amplitude outsideof an assumed tangent cylinder is reduced by applying an additionalattenuation function (eq. 30). The main parameters to examine arethe latitude of the TC, 𝛽 𝑇 𝐶 , and the decay depth ℎ (fig. 7). Likein fig. 5 the transparent stripes show where individual gravity har-monics are in agreement with the observed values (Durante et al.2020). The TC decay width and the radial decay width are fixed to 𝛿 𝛽 = ◦ and 𝛿 =
500 km respectively. For small 𝛽 𝑇 𝐶 the solutionis not affected, implying that the antisymmetric zonal flow betweena latitude of ± ◦ is irrelevant for the gravity signal. This changesonce the TC angle is increased beyond 𝛽 𝑇 𝐶 = ◦ since we start re-ducing the amplitude of the dominant prograde jet at about 𝛽 = ◦ latitude. The required decay depth ℎ ℓ decreases for ℓ = ℓ =
5, 7, and 9. At ℎ = 𝛽 𝑇 𝐶 = . ◦ themodel can finally explain all observed gravity harmonics (see also tab. 1). The corresponding surface flow of the TC model is shownin the inset of fig. 7 (blue profile) indicating that the dominant jetaround 20 ◦ is thinner than in the reference model (black) and re-duced to a peak amplitude of 35 m/s. Also Kong et al. (2018) foundit necessary to reduce the amplitude of the prograde jet.Reducing the flow in an equatorial latitude band reduces theimpact of the equatorial step on the density anomaly. It even vanisheswhen the TC is positioned at sufficiently large latitudes such thatthe flow is split into two independent parts. This must be consistentwith the assumed radial decay function 𝑄 𝑟 : a deeper-reaching flow(larger ℎ or 𝛿 ℎ ) requires a larger 𝛽 TC to separate the northern andsouthern hemisphere. Thus, a geometric relation between radialdecay function 𝑄 𝑟 and the width of the equatorial cut-off 𝛽 TC canbe formulated. The dark curves in fig. 7 indicate the position of a TCattached to various depths; i.e. where the associated drop-off equals10 − (light grey), 10 − (dark grey) or 10 − (black). Interestingly,the black line exactly hits the intersection point of the colouredisocontours. This means that the TC defined by reducing the flowamplitude in an equatorial band of ± . ◦ coincides with the TCdefined by the 10 − -drop-off of the radial decay function 𝑄 𝑟 . Thusthis marks the point in the ( 𝛽 TC - ℎ )-parameter space, where the flowis just split into two separate hemispheric flows at the minimumpossible 𝛽 TC .Moreover, the other two parameters, 𝛿 ℎ and 𝛿 𝛽 , can be con-strained with this result. Higher values of 𝛿 𝛽 do not lead to anintersection of the individual gravity solutions (colours in fig. 7),whereas larger or smaller values of 𝛿 ℎ shift the black line out of theintersection point.Consequently, the TC gravity model not only explains all grav- MNRAS , 1– ?? (202-) W. Dietrich et al.
Figure 7.
Two-dimensional parameter scan with respect to the decay depth ℎ and the assumed TC angle 𝛽 𝑇 𝐶 . The coloured isocontours indicate wherethe modelled gravity moments agree with the measurements and its uncer-tainties. The dark lines show the surface latitude of a TC attached to tenfolddrop-offs of the associated decay function. The optimal parameter choice is 𝛽 TC = . ◦ and ℎ = ity measurements; it is also independent of the handling of theequatorial step since the geostrophic extension of the antisymmet-ric flows in either hemisphere do not reach the equatorial plane.Fig. 6 d) illustrates how the gravity harmonics change for 𝛽 𝑇 𝐶 = . ◦ when increasing ℎ . The TC scenario most stronglyaffects 𝐽 which now remains negative for all ℎ . However, the otherharmonics are more significantly changed than by the different ap-proaches to treat the equatorial discontinuity illustrated in fig. 6.The subsequent ℎ is much smaller than in other models, whereasthe other ℎ increase drastically. All ℎ ℓ are in agreement within theuncertainties at approximately ℎ = 𝐽 Evidently the 𝐽 gravity moment is highly sensitive to different mod-els for the equatorial treatment. To understand this, we analyse therespective source - the dynamic density perturbation 𝜌 𝑈 - in moredetail. Separating the first term in eq. 26 into two and transformingin Legendre space yields three source contributions for ℓ = 𝜌 𝑈 = − Ω 𝑟𝑔 𝑄 ( 𝜕 𝑟 𝜌 ) F (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) 𝜌 𝑈 ,𝑐 − Ω 𝑟𝑔 𝜌 ( 𝜕 𝑟 𝑄 ) F (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) 𝜌 𝑈 ,𝑑 − Ω 𝑔 𝜌𝑄 H (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) 𝜌 𝑈 ,𝑒 , (45)where F and H exhibit the Legendre-transformations of the inte-grated flow and the step contribution: F = ∫ 𝜋 𝑃 ∫ 𝜃 𝑈 𝑜 cos ˆ 𝜃𝑑 ˆ 𝜃 sin 𝜃𝑑𝜃 (46) H = Δ 𝑈 ∫ 𝜋 𝑃 𝐻 ( 𝜃 − 𝜋 / ) sin 𝜃𝑑𝜃 . (47)The first part, 𝜌 𝑈 ,𝑐 mainly represents the geostrophic part of theflow (where 𝑄 𝑟 is constant), the second term, 𝜌 𝑈 ,𝑑 is dominatedby the decay region and the last, 𝜌 𝑈 ,𝑒 shows the influence of theequatorial step. In fig. 8 the left panel shows the resulting profiles forthe reference model with ℎ = 𝐽 . The right panel displays the TC modelwith ℎ = F and H are plotted inthe insets of the respective model.For the reference model, the first contribution 𝜌 𝑈𝑐 (blue) isproportional to the radial derivative of 𝜌 and is therefore negativesince F is also negative. It dominates at shallow depths 𝑟 > . 𝑅 .The second source contribution 𝜌 𝑈𝑑 (orange) is proportional to theradial derivative of the depth profile 𝑄 𝑟 and therefore positive. It isonly significant between 0 . 𝑅 and 0 . 𝑅 where the large ℎ -valuepositions the steep decrease in 𝑄 𝑟 . The flow amplitude representedby F is rather small at this depth hence the small positive bump in 𝜌 𝑈𝑑 . The third contribution, 𝜌 𝑈𝑒 (green) is mostly positive (due to H ) and thus has the wrong sign for explaining the negative 𝐽 . Fora successful model, the negative 𝜌 𝑈𝑐 has to overcompensate 𝜌 𝑈𝑑 and 𝜌 𝑈𝑒 which only succeeds when 𝜌 𝑈𝑐 reaches deep enough. Both 𝜌 𝑈𝑐 and 𝜌 𝑈𝑒 individually would yield much larger 𝐽 amplitudes thanrequired. Only the delicate balance between the three contributionyields the correct sign and amplitude. This explains the sensitivityof 𝐽 to the reference model setup.The source analysis also readily explains why ignoring 𝜌 𝑈𝑒 inthe hemispherical model (panel b) of fig. 6) allows one to explain 𝐽 with a much smaller value of ℎ . The linear 𝑧 -dependence assumedfor panel c) of fig.6 only decreases the amplitude of 𝜌 𝑈𝑒 and istherefore less effective.For the TC model (right panel of fig 8) the equatorial stepcontribution, 𝜌 𝑈 ,𝑒 is identical to zero validating the results from theprevious chapter (see also fig. 6, d). The density anomaly entirelycomprises of a shallow, negative 𝜌 𝑈𝑐 and deep-rooted positive decaypart 𝜌 𝑈𝑑 . Evidently, 𝜌 𝑈𝑑 is larger because F is larger at the depth ofthe decay (close to 0 . 𝑅 ). We have mentioned above that some authors ignore the dynamicself-gravity or DSG-term (Kaspi et al. 2018; Galanti & Kaspi 2020),which allows them to solve the TWE rather than the TGWE equa-tion. We model TWE-based gravity moments by setting 𝜇 = 𝜇 -profile we model the constant coefficient TGWE equation usedby Zhang et al. (2015); Kong et al. (2018); Wicht et al. (2020) bycalculating gravity moments with a constant 𝜇 = 𝜋 / 𝑅 .Tab. 1 lists the decay depth ℎ ℓ that is required to to match theobserved gravity signal. Ignoring the DSG term (TWE model, or 𝜇 =
0) shows a reduction of ℎ by almost 1000 km (or 15%) and anincrease of ℎ by 170 km (10%). In general, the absolute differenceshrinks with the degree ℓ . The model with 𝜇 = 𝜋 / 𝑅 shows anintermediate behaviour indicating that a more realistic radial profileof 𝜇 is an important part of the model.Alternatively, we calculate the deviation of the gravity mo-ments 𝐽 ℓ between different models for fixed degree ℓ and decaydepth ℎ . We quantify the error in 𝐽 ℓ -error of neglecting the DSGterm as a function of decay depth ℎ in terms of 𝑒 TWE ℓ ( ℎ ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − 𝐽 TWE ℓ ( ℎ ) 𝐽 ref ℓ ( ℎ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (48)The measure 𝑒 cTGWE ℓ , defined in analogy to eq. 48, quantifies theerror for the constant coefficient TGWE. Fig. 9 shows how botherrors change for our reference model when varying ℎ between1000 km and 5000 km (we exclude larger ℎ values because here 𝐽 changes sign and the definition in eq. 48 becomes useless). Except MNRAS , 1– ????
0) shows a reduction of ℎ by almost 1000 km (or 15%) and anincrease of ℎ by 170 km (10%). In general, the absolute differenceshrinks with the degree ℓ . The model with 𝜇 = 𝜋 / 𝑅 shows anintermediate behaviour indicating that a more realistic radial profileof 𝜇 is an important part of the model.Alternatively, we calculate the deviation of the gravity mo-ments 𝐽 ℓ between different models for fixed degree ℓ and decaydepth ℎ . We quantify the error in 𝐽 ℓ -error of neglecting the DSGterm as a function of decay depth ℎ in terms of 𝑒 TWE ℓ ( ℎ ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − 𝐽 TWE ℓ ( ℎ ) 𝐽 ref ℓ ( ℎ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (48)The measure 𝑒 cTGWE ℓ , defined in analogy to eq. 48, quantifies theerror for the constant coefficient TGWE. Fig. 9 shows how botherrors change for our reference model when varying ℎ between1000 km and 5000 km (we exclude larger ℎ values because here 𝐽 changes sign and the definition in eq. 48 becomes useless). Except MNRAS , 1– ???? (202-) onal Winds and Gravity reference model TC model Figure 8.
Radial profiles of the dynamic density source for degree ℓ =
3, separated into three parts (constant, decay and equatorial step part, see eq. 45). Theleft panel shows the reference model with ℎ = ℎ = F and H . Figure 9.
Error of the modelled gravity moments for the TWE (solid)and constant coefficient TGWE solver (dot-dashed) relative to the referencemodel as function of the decay depth. Colors indicate different odd degrees. for ℓ = ℎ in this range. For ℓ = 𝑒 TWE3 always exceeds 60%. For ℓ = ℓ = ℓ =
9, the error reachesabout 10% for larger values of ℎ . These errors are even larger thanquantified by Zhang et al. (2015) or Wicht et al. (2020) and showthat the DSG term with a radius dependent DSG coefficient shouldnot be neglected. The cTGWE approach, i.e. where 𝜇 is a constant,captures roughly half of the effect of the full TGWE approach. Forexample, for ℓ = 𝑒 cTGWE3 amounts to 40%, whereas 𝑒 TWE3 exceeds 70% at ℎ = We have discussed above that the different zonal wind models pub-lished over the years likely are likely to reflect seasonal variations,which may be restricted to the very shallow atmosphere and there-fore do not affect the gravity signal. This motivates us to use a meanflow averaged over the available profiles from different epochs (seesec. 3.2 and fig. 4).Kaspi et al. (2018) relied on the HST observations obtainedduring Juno’s perijove 9 pass (Tollefson et al. 2017), while Kong et al. (2018) based their analysis on the Cassini measurements(Porco et al. 2003). We implement both, the Perijove and the Cassinisurface flow model, in our standard setup and additionally apply theamplitude cap by Kong et al. (2018) in a model we call ‘cap‘. Tab. 1lists the change in ℎ ℓ required when using either flow and shows thatthe largest influence is again found for degree ℓ =
3, where in partic-ular the Cassini measurements require deeper sources. The Perijoveprofile, on the other hand yields larger offsets in the higher degrees.The results of Kong et al. (2018) are based in a modified version ofthe Cassini profile (Porco et al. 2003), i.e. capping the large spikeat 20 ◦ -latitude to a maximum amplitude of 𝑈 spike = . / s. Thisflow modification is a necessary part of their model in order to matchthe gravity moments. Tab. 1 show that the major effect is once morea reduction of ℎ and a smaller increase for the higher degrees. Noteonce again that the successful reproduction of Jupiter’s odd-degreegravity moments by Kong et al. (2018) is due to a combinationof applying this particular flow modification, the additional linear 𝑧 -dependence, the smooth radial decay function, assuming a poly-tropic interior and solving the TGWE.Generally, the impact of different surface flow profiles reaches10% and is thus smaller than, for example, the effect of the equatorialstep or neglecting the DSG. Finally we also test the impact of using different interior structuremodels. We replace the model by Guillot et al. (2003) in our ref-erence setup by the simple polytrope of index unity or the modelby Nettelmann et al. (2012). Details of the models are outlined insec. 3.1. The chosen interior model factors in the dynamic densitysource (eq. 26) as the background density and gravity, 𝜌 and 𝑔 , andin the DSG coefficient 𝜇 (eq. 13). Fig. 2 shows the relevant profilesof density, density gradient and 𝜇 . Tab. 1 shows that the impact onthe ℎ ℓ is even smaller than the impact of the surface flow model,always remaining below 10%. This indicates that the odd-degreegravity moments are merely insensitive to details of the Jupiter’s in-terior structure, such as properties of the core or the heavy elementabundance. MNRAS , 1– ?? (202-) W. Dietrich et al.
Table 1.
List of the models with the different treatment of the equator, different solvers, interior models and surface flows. Further given is the decay depth ℎ ℓ at which the modelled gravity moment agrees with the observations.name equator solver interior surface flow ℎ [km] ℎ [km] ℎ [km] ℎ [km] reference full TGWE G03 mean 6541 + − + − + − + − hemispheric hemispheric TGWE G03 mean 3179 + − + − + − + − linear z linear z TGWE G03 mean 5127 + − + − + − + − TC TC TGWE G03 TC cut 2950 + − + − + − + − TWE full TWE G03 mean 5670 + − + − + − + − cTGWE full cTGWE G03 ∗ mean 6190 + − + − + − + − perijove full TGWE G03 HST, perijove 6480 + − + − + − + − Cassini full TGWE G03 Cassini 7005 + − + − + − + − cap full TGWE G03 capped Cassini 6263 + − + − + − + − poly full TGWE polytropic mean 6213 + − + − + − + − N12 full TGWE N12 mean 6215 + − + − + − + − Table 2.
Jupiter’s odd-degree gravity moments from our TC model ( 𝐽 ℓ )and the observations ( 𝐽 obs ℓ ) alongside the 3 𝜎 uncertainties (Durante et al.2020). degree ℓ 𝐽 obs ℓ [10 − ] 𝐽 ℓ [10 − ]3 -4.5 ± ± ± ± ± For the first time, the precise measurements of Jupiter’s gravity fieldby the Juno mission allowed quantifying the tiny undulations causedby the zonal winds. First attempts to invert the odd-degree gravitymoments suggested a rather smooth decay with depth, reachingabout 2 m/s or even 20 m/s at a radius of 𝑟 = . 𝑅 (Kaspi et al.2018; Kong et al. 2018). This seems at odds with the magnetic ob-servations, suggesting that the speed at this depth should not exceedthe centimetre per second level (Moore et al. 2019). The gradualdecay is also difficult to reconcile with the results from numericalsimulations and theoretical considerations. They suggest the windsshould remain geostrophic in the outer part of the atmosphere untilthey are abruptly quenched at some depth (Christensen et al. 2020;Gastine & Wicht 2021).In addition to the too smooth radial decay, the handling of thediscontinuity at the equatorial plane that can arise when asymmet-ric surface winds are downward continued in z-direction (parallelto the rotation axis). Such a discontinuity is unphysical, but if itexisted it would contribute substantially to the modelled gravitysignal. To avoid this problem, Kong et al. (2018) therefore some-what arbitrarily used an additional linear variation of the zonal with 𝑧 to eliminate the discontinuity, which is obviously at odds witha largely geostropic flow. Kaspi et al. (2018) chose an alternativeapproach and simply ignore the contribution originating from the discontinuity. Our analysis suggest that this was essential to success-fully model the observed odd gravity harmonics but at the expenseof physical consistency. More recently, Galanti & Kaspi (2020) sug-gest a different radial profile that seems consistent with the magneticobservations consisting of a geostrophic outer envelope and a steeprapid decay at 𝑟 = . 𝑅 . They also adopt the arguable approachof ignoring the equatorial discontinuity from Kaspi et al. (2018).The differences in the findings might be caused by various differ-ent modelling assumptions regarding how the zonal flow relates toanomalies on the gravity potential, the treatment of the equatorialstep, the applied surface zonal flow model and the interior model.Here we aim to quantify the impact of each of those.We model the odd gravity moments of Jupiter by using thethermo-gravitational wind equation (TGWE) which correctly de-scribes the flow induced perturbations of the gravity potential. Themajor difference to the more commonly used thermal wind equationis the neglect of one gravity force contribution, the dynamic self-gravity (DSG). So far solving the TGWE has been possible onlyfor polytropic interior models with index unity Zhang et al. (2015);Wicht et al. (2020), thus we introduce a method to solve a moregeneral form that allows to use realistic interior state models andcaptures the important radial variation of the DSG coefficient.This enables us to measure the impact of the DSG term thatseveral authors neglect in their study (Kaspi et al. 2018; Galantiet al. 2017; Galanti & Kaspi 2020). Our results confirm that theDSG is indeed a first order effect as already predicted by Zhanget al. (2015), Kong et al. (2018) or Wicht et al. (2020). The gravityharmonic 𝐽 is most severely affected and changes by up to 90% inthe models explored here. The impact then decreases with sphericalharmonic order and remains below 10% for 𝐽 . Moreover, we showthat impact of the radius-dependent form of the DSG coefficient isroughly twice as large as in the polytropic (constant DSG coefficient)form of the TGWE (Wicht et al. 2020).Presumably the wind structure responsible for generating thegravity anomalies is of larger scale and rather stationary, whereasall of the published jovian surface zonal wind profiles are snapshotswith respect to the Jovian orbital period of 12 yrs. We thus generatean average zonal flow profile that represent the deep flow structure
MNRAS , 1– ????
MNRAS , 1– ???? (202-) onal Winds and Gravity more faithfully. Our results suggest that using a specific individualprofile rather than the mean significantly biases the analysis.We also quantify the impact of different interior models andfind that its influence is smaller than the equatorial treatment, theDSG term or a specific flow profile. Hence constraining the proper-ties of the deep interior structure (Guillot et al. 2018), such as coremass or composition, on the basis of the odd gravity moments seemsrather challenging as they are far more sensitive to other model as-sumptions, such as the equatorial treatment, considering the DSGterm or assumptions about the structure of the zonal winds.Here we suggest an alternative scenario that is motivated bythe observation that the equatorially antisymmetric flow contribu-tions are much weaker outside than inside the tangent cylinder (TC)in numerical simulations. The TC is the imaginary cylinder with aradial distance from the rotation axis that is defined by the depth atthe equator where deviations from a barotropic state and/or Lorentzforces enforce a sharp drop of the zonal flow. This seems consistentwith the simple fact that equatorially antisymmetric flows are nec-essarily ageostrophic outside the TC, while there is no such conflictinside the TC. In our model we therefore considerably damp theantisymmetric flow outside the TC, smoothing the transition with atangent hyperbolic function. This relieves us from the problem ofthe equatorial discontinuity without having to assume an unrealistic 𝑧 -dependence of the flow. We thus favour a model where the an-tisymmetric surface flow in an equatorial latitude band of ± ◦ isa shallow phenomena, not contributing to the gravity signal. This,in turn, requires strong baroclinicity in the first few hundred kilo-meters below the cloud level generated from latitudinal gradients intemperature or composition. Interestingly, such pronounced equato-rial variations in the nadir brightness temperature and the ammoniamixing ratio have been detected by the Juno mission (Bolton et al.2017).We use the four significant odd moments, 𝐽 to 𝐽 to indi-vidually constrain the vertical and latitudinal structure of the zonalwinds. We can match the modelled gravity moments to the ob-servations within the 3 𝜎 uncertainties (Durante et al. 2020) whenusing a radial decay profile that models a geostrophic flow insideTC for 𝑟 > . 𝑅 but then decays rapidly with depth (see alsotab. 2). Outside TC the flow must be shallow. Our preferred decayprofile suggests that the winds are 50% deeper than suggested byGalanti & Kaspi (2020) and hence the geostrophic part reaches tosubstantially larger pressures (830 kbar rather than 300 kbar). At 𝑟 = . 𝑅 the profile reaches finally 10 − what is in agreementwith the constraints from Moore et al. (2019). This depth equalsalso to the lower boundary of a TC touching the surface at ± ◦ implying that our preferred wind structure exhibits a splitting ofthe gravity perturbation into an independent northern and southernhemisphere, thus the contribution from the equatorial step naturallyvanishes.In order to explain the odd gravity moments with a physicallyconsistent handling of the flow near the equatorial plane, we find itnecessary to exclude part of the observed cloud level zonal windsfrom the analysis by assuming that it is restricted to shallow depth.Similarly, Kong et al. (2018) had to cap the amplitude of the strongjet centred at a latitude of 𝛽 = ◦ N to fit the gravity moments intheir model. While it is physically reasonable that deep antisymmet-ric winds do not exist outside the tangent cylinder (defined with thedepth at the equator where deep winds drop sharply in amplitude),there is obviously no guarantee that the antisymmetric cloud-levelwinds inside the tangent cylinder are fully representative of deepflow and do not also comprise a shallow component. This uncer-tainty is probably the severest limitation for utilising the gravity anomalies to infer the precise depth extent of the zonal flow. A bet-ter understanding of the atmospheric dynamics at different latitudesin the top few hundred kilometres below the clouds is necessary toovercome this limitation.
ACKNOWLEDGEMENTS
We thank Eli Galanti for very constructive discussions. This workwas supported by the Deutsche Forschungsgemeinschaft (DFG)in the framework of the priority program SPP 1992 ’Diversity ofExoplanets’. The MagIC-code is available at an online repository(https://github.com/magic-sph/magic).
REFERENCES
Bolton S. J., et al., 2017, Science, 356, 821Braginsky S. I., Roberts P. H., 1995, Geophysical and Astrophysical FluidDynamics, 79, 1Cao H., Stevenson D. J., 2017, Journal of Geophysical Research: Planets,122, 686Christensen U. R., 2002, Journal of Fluid Mechanics, 470, 115Christensen U. R., Wicht J., Dietrich W., 2020, ApJ, 890, 61Dietrich W., Jones C. A., 2018, Icarus, 305, 15Duer K., Galanti E., Kaspi Y., 2019, ApJ, 879, L22Durante D., et al., 2020, Geophys. Res. Lett., 47, e86572French M., Becker A., Lorenzen W., Nettelmann N., Bethkenhagen M.,Wicht J., Redmer R., 2012, The Astrophysical Journal Supplement Se-ries, 202, 5Galanti E., Kaspi Y., 2020, Monthly Notices of the Royal AstronomicalSociety, 501, 2352Galanti E., Kaspi Y., Tziperman E., 2017, Journal of Fluid Mechanics, 810,175–195Galanti E., Kaspi Y., Miguel Y., Guillot T., Durante D., Racioppa P., Iess L.,2019, Geophysical Research Letters, 46, 616Gastine T., Wicht J., 2012, Icarus, 219, 428Gastine T., Wicht J., 2021, MNRASGuillot T., Chabrier G., Morel P., Gautier D., 1994, Icarus, 112, 354Guillot T., Stevenson D. J., Hubbard W. B., Saumon D., 2003, The Interiorof Jupiter. Springer International Publishing (Cham), pp 35–58Guillot T., et al., 2018, Nature, 555, 227Heimpel M., Aurnou J., Wicht J., 2005, Nature, 438, 193Heimpel M., Gastine T., Wicht J., 2016, Nature Geoscience, 9, 19Hubbard W. B., 1999, Icarus, 137, 357Iess L., et al., 2018, Nature, 555, 220Iess L., et al., 2019, ScienceJones C. A., 2007, in Olson P., ed., Treatise on Geophysics, Treatise onGeophysics Volume 8: Core Dynamics. Elsevier, pp 131–185Kaspi Y., Hubbard W. B., Showman A. P., Flierl G. R., 2010, Geo-phys. Res. Lett., 37, L01204Kaspi Y., Showman A. P., Hubbard W. B., Aharonson O., Helled R.,2013, in AAS/Division for Planetary Sciences Meeting Abstracts , 1– ?? (202-) W. Dietrich et al.
Nettelmann N., Becker A., Holst B., Redmer R., 2012, The AstrophysicalJournal, 750, 52Porco C. C., et al., 2003, Science, 299, 1541Schaeffer N., 2013, Geochemistry, Geophysics, Geosystems, 14, 751Taylor G. I., 1917, Proc.~R.~Soc.~Lond.~A, 93, 99Tollefson J., et al., 2017, Icarus, 296, 163Wicht J., 2002, Physics of the Earth and Planetary Interiors, 132, 281Wicht J., Gastine T., 2020, Nature Communications, 11, 2886Wicht J., Gastine T., Duarte L. D. V., Dietrich W., 2019, A&A, 629, A125Wicht J., Dietrich W., Wulff P., Christensen U. R., 2020, MNRAS, 492,3364Zhang K., Kong D., Schubert G., 2015, ApJ, 806, 270 MNRAS , 1– ????