Lower‐order zonal gravitational coefficients caused by zonal circulations inside gaseous planets: Convective flows and numerical comparison between modeling approaches

To infer the internal equilibrium structure of a gaseous planet, especially the equation of state (EOS) and size of its inner core, requires accurate determination of lower‐order zonal gravitational coefficients. Modeling of the gravitational signature associated with deep zonal circulation depends critically upon reliable subtraction of the dynamical components from totally derived gravitational coefficients. In the era of the Juno mission and the Grand Finale phase of the Cassini mission, it is timely and necessary to revisit and examine the so‐called ‘Thermal Wind Equation (TWE)’, which has been extensively utilized to diagnose the dynamical parts of the gravitational fields measured by the two spacecrafts. TWE treats as negligible a few terms in the full equation of balance. However, the self‐gravitational anomaly of the distorted fluid, unlike oblateness effects of solid‐body rotation, is not a priori minor and thus should not be neglected in the name of approximation. Another equation, the ‘Thermal Gravitational Wind Equation (TGWE)’, includes this important additional term; we compare it with the TWE and show that physically the TGWE models a fundamentally different balance from the TWE and delivers numerical results considerably different from models based on the TWE. We conclude that the TWE balance cannot be relied upon to produce realistic convection models. Only after the TGWE balance is obtained can the relative importance of terms be assessed. The calculations we report here are based on two types of zonal circulations that are produced by realistically possible convections inside planets, instead of being constructed or assumed.


Introduction
The Juno mission and the Grand Finale stage of the Cassini mission have achieved unprecedentedly accurate measurements of the gravitational fields of Jupiter and Saturn (Bolton et al., 2017;Edgington and Spilker, 2016). It has long been expected that highly accurate gravitation data would greatly improve our knowledge of internal states of gas giants (Stevenson, 1982;Guillot, 2005;Anderson and Schubert, 2007;Nettelmann et al., 2012;Kong DL et al., 2016a). After data became available (Durante, 2017;Iess et al., 2018Iess et al., , 2019, many efforts have been made to interpret the measurements in order to ribe the mysterious internal structures and equations of state (EOS) of these two giant planets (Miguel et al., 2016;Militzer et al., 2016;Wahl et al., 2017;Ni DD 2018;Iess et al., 2019;Debras and Chabrier, 2019;Kong DL et al., 2019). Such an objective demands a reliable knowledge of the hydrostatic component of the gravitational fields. However, it should be noted that a spacecraft measures only the total gravitational field, which includes dynamical signatures associated with internal fluid motions, including zonal circulations that can substantially perturb zonal gravitational coefficients. To extract the main hydrostatic component for studies of internal structure, accurately modeling the dynamical part becomes crucial; otherwise, uncertainties in flow-related gravitation can make it difficult to distinguish between different internal models (Wahl et al., 2017;Ni DD, 2019). The two most important sources of such uncertainties are (1) an unknown deep zonal circulation profile below the cloud-top, and (2) how to diagnose accurately the distortion caused by a zonal circulation.
Jupiter and Saturn are marked by fast zonal winds observed at their cloud-top (Porco et al., 2003). It has remained disputable (Kaspi et al., 2018;Kong DL et al., 2018b) whether such zonal wind extends deeply into the planetary interiors until being stopped by magnetic braking effects, or whether the wind represents just very shallow atmospheric dynamics. In the shallow scenario, there must be, beneath the zonal wind layer, a deeper circulation of a different profile but dynamically coupled with the upper jet stream (Dowling, 1995;Thomson and McIntyre, 2016). In either case, zonal flow exists in the conductivity-low envelope of a gaseous planet, produced by thermal convection under fast rotation (Zhang KK and Schubert, 2000;Zhang KK and Liao X, 2017). The fluid motion causes a density anomaly that further distorts the gravitational field. Externally, such a distorted planetary gravitational potential can be expanded by in which is the gravitational constant, is the planetary mass, is the equatorial radius of the planet. In the spherical polar coordinate system (r, θ, φ), r and θ are the radial and angular coordinates and the origin, r = 0, is set at the center of figure; θ = 0 marks the symmetry axis of rotation. P n (cos θ), n = 2, 3, ··· are Legendre polynomials. J n are zonal gravitational coefficients linked to the hydrostatic equilibrium structure of the planet, under a uniform rotation Ω; ∆J n are dynamical gravitational signatures associated with the internal zonal circulation, either equatorially symmetric (resulting in ∆J 2k , k = 1, 2, ···) or anti-symmetric (resulting in ∆J 2k+1 , k = 0, 1, 2, ···). g In modeling the zonal-circulation-related distortion, several authors have employed the so-called Thermal Wind Equation (TWE) for density ρ, gravity and pressure p (Kaspi et al., 2018;Iess et al., 2019;Debras and Chabrier, 2019;Galanti et al., 2019), where u is the zonal circulation velocity, Ω is the angular speed of planetary rotation, and the z-axis is parallel to the rotation axis. The subscript 0 denotes the leading-order hydrostatic equilibrium variables while the primed variables represent dynamical distortion. The TWE Equation (2) is an incomplete balance that has neglected a few terms in the full equation (to be discussed in the next section), some of which, such as the geometric non-spherical effects caused by rotation, can indeed be regarded, a priori, as negligible. But the self-gravitational anomaly is not only of a fundamentally different mathematical nature but also is not necessarily small. On one hand, the two terms, namely and , are generally of the same order of magnitude. It is mathematically valid to compare their relative importance only by solving the equation that includes both terms. However, to judge the size of the based on the partial balance obtained through the TWE is highly questionable, e. g. Galanti et al., (2017); Kaspi et al., (2019). Later in this article, it will be shown that when zonal circulations are not restricted to the shallow outer atmosphere but allowed in-g ′ ρ 0 g ′ side a significant bulk of interior, the two terms are indeed equally important numerically. On the other hand, from consideration of the physics and the actual measurement data, the gravitational signature is the quantity of direct interest. It is inappropriate to delete from the equation the very term that is supposed to be measured. Last but not least, Kong DL et al., (2016b) has demonstrated that even more profound algebraic differences result when the term is introduced.
In this article, modeling the dynamical gravitational signatures of zonal circulations is discussed. Adopted deep zonal circulations are produced by possible rotating thermal convections but not artificially constructed from cloud-top winds. The significant role of the self-gravitational anomaly is discussed and highlighted in a clear-cut manner. This important additional physics is included in the 'Thermal Gravitational Wind Equation (TGWE)' (Zhang KK et al., 2015), whose name reflects that inclusion; we compare it with the TWE and demonstrate a large numerical discrepancy between the two approaches when they are applied to the modeling of lower-order zonal gravitational coefficients of gaseous planets.
In what follows, the full Euler equation and relevant approximations are introduced and analyzed in Section 2. Then in Section 3, two examples are discussed of zonal circulations that could be produced by realistic rotating thermal convections. The TGWE and TWE models are directly compared in solving gravitational distortions caused by the flows. Conclusions and discussions appear in Section 4.

Models and Equations
g Since cloud-top zonal winds on Jupiter and Saturn have been fairly steady over decades of observations, and the flows have typically been characterized by small Rossby number, the full Euler equation for density ρ, self-gravity , and pressure p can be written as Precisely speaking, the domain is essentially non-spherical. The planetary outer level is an equipotential surface, described by Equation (8), whose departure from sphericity depends on the EOS Equation (4) and a small rotation parameter in which is the mean density of the planet. Because typical zonal flow speed is small relative to the solid-body rotation, namely , for both Jupiter and Saturn, Equations (3)-(8) can be analyzed via a perturbation approach 90 Earth and Planetary Physics doi: 10.26464/epp2020014 u u u where the subscript 0 denotes variables in the leading order of hydrostatic equilibrium under the solid-body rotation, and the primed variables are dynamical distortions caused by the zonal circulation . The leading order problem, though not of our primary concern in this article, needs to be explained, about a few important approximations that are carried to the next order. Substitution of Equations (9)-(11) into Equations (3)-(8) yields the equilibrium equations free from flow distortion, (12)-(16) and corresponding domain is numerically complicated (Hubbard 2013;Kong DL et al., 2015a;Nettlemann, 2017). However, in this study, our aim is merely to compare TWE and TGWE in the next order problem, under the assumption that the rotational effects resulting from the centrifugal force are neglected for two reasons: first, the rotational parameter is generally small for Jupiter and Saturn; second, departure from spherical symmetry is an a priori known effect that is common for both TWE and TGWE. Therefore, the domain is approximated by a sphere whose radius R is the mean radius of the planet; hence all leading-order variables are functions of radius r only. Another approximation involved is the effective polytropic EOS adopted as

Self-consistently solving Equations
for Equation (13), which has been used in several previous studies of gaseous planets (Hubbard, 1999;Kong DL et al., 2015b). In reality, the physical EOS of a giant planet can be rather complicated (Militzer et al., 2016). However, an effective polytropic EOS can always be found to best approximate the true EOS, following the polytropic equilibrium condition discussed in Chandrasekhar, (1939). The effective polytropic index for Jupiter has been found to be very close to unity (Horedt, 2006). Under the two approximations, Equations (12)-(16) adopt the solutions where ξ = r/L, The next order problem for perturbations and is defined in the same spherical domain.
Note that it can be clearly seen that The small deviation of the polytropic EOS from the physical EOS will only incur a second-order effect in Equation (19) and hence is neglected. Generally, there is no any a priori reason that the term should be small before the complete equations, Equations (19)- (20), are solved, unless under some very special circumstances when is not concerned or the equations are solved within a given very shallow layer of atmosphere. Applying to both sides of Equation (19), considering the zonal flow in the form of , we end up with the so-called 'Thermal-Gravitational Wind Equation', in which the arbitrary function C(r) does not enter the computation of zonal gravitational coefficients The TGWE is algebraically different from the TWE because of the kernel term (Kong DL et al., 2016b;Zhang KK et al., 2017a). A numerical scheme to solve the integral equation has been discussed at length in Zhang KK et al., (2015).

Numerical Comparison Between TGWE and TWE
To demonstrate that TGWE and TWE yield major numerical differences, we have examined two typical types of zonal circulation inside a planet, an equatorially symmetric one and an equatorially anti-symmetric one. The adopted flows are neither constructed artificially nor extended from cloud-top winds, but are linked physically to possible thermal convections operating inside a planet. In both cases, the common parameters are listed in the Table 1. For illustrative purpose, the convective motion in a sphere is considered under the regime of Boussinesq approximation, which is still consistent with Equation (6) because all further discussions will be restrained within the context of zonal circulations. Consider the case of thermal convection in a fluid sphere with constant thermal diffusivity , thermal expansion coefficient , and kinematic viscosity (Kong DL et al., 2018a). The fluid sphere rotates uniformly with constant angular velocity in the presence of its own gravitational field , where is the position vector. The whole sphere is heated by a uniform distribution of heat sources, producing the unstable conducting temperature gradient , β being a positive constant. When β is sufficiently large, convective instability takes place and drives fluid motion in the sphere. Upon employing the radius of the sphere R as the length scale, Ω -1 as the unit of time, and as the unit of temperature fluctuation, the problem of thermal convection is governed by the dimensionless equations (Pr/Ek) (

∂Θ ∂t
where the three controlling dimensionless parameters, the Rayleigh number, Ekman number, and Prandtl number, are defined as The convection can be subject to either a no-slip boundary or a stress-free boundary condition. When, as is true for fast-rotating Jupiter and Saturn, the Ekman number is asymptotically small , the Prandtl number determines the regime of dynamics.

Equatorially Symmetric Differential Rotation in the Regime
In the cases of , convective motions in the form of spiraling columnar rolls are marked by large azimuthal variations ( ) from the linear analysis. Cylindrical differential rotation , which is equatorially symmetric, can be generated from small-scale convection velocity by the nonlinear Reynolds stress at the onset of convection. Solution of Equations (23)-(25) can be expanded by where represents a small amplitude, with . Substitution of Equations (26)-(29) into Equations (23)-(25) yields the linear convective motion at onset . The nonlinearly induced differential rotation then can be derived by For the equatorially symmetric case of zonal circulation, the differential rotation of the convection with Ek = 5×10 -5 and Pr = 0.025 is adopted (see the Section 19.2.4 of Zhang KK and Liao X, (2017)). Figure 1 depicts the zonal flow speed as a function of distance from the rotation axis. Table 2 gives the polynomial coefficients that recover the curve in Figure 1. The equatorially symmetric flow is respectively inserted into the TWE Equation (2) and the TGWE Equation (21) to compute and associated gravitational distortions . The r−θ grid discretization for numerical computations (see Zhang KK et al., (2015)) is 200 × 401. Results are presented in Table 3. As one can immediately see from the comparison, lower order zonal gravitational coefficients can differ by O(100%), as predicted by theoretical analysis. The implication of the contrast is that the term missing in the TWE is important for this particular problem.

Pr/Ek
For a moderate value of , both the inertial effect and the viscous effect contribute to the rotating convection. Equatorially anti-symmetric zonal circulation, in the form of torsional oscillation, can be yielded in some very special convection regimes (Zhang KK et al., 2017b). We adopt a simple but physically possible flow (e. g. the bifurcation A presented and discussed in Kong DL et al., (2018a)) For this zonal flow, Table 4 presents the results of gravitational coefficients computed respectively by TWE and TGWE. It can be seen that again there are significant numerical differences in The equatorially symmetric zonal circulation profile. The mean zonal flow is generated via nonlinear process from columnar convective motion whose Ekman number is 5×10 -5 and Prandtl number is 0.025.

92
Earth and Planetary Physics doi: 10.26464/epp2020014 lower order zonal gravitational coefficients. It is worth mentioning that there is no contribution from the hydrostatic equilibrium structure ρ 0 in odd-order gravitational coefficients (Kong DL et al., 2016b). Geometric flattening effects from fast rotation pose negligible influence. The comparison between TGWE and TWE hence is more straightforward. Huge numerical difference between TWE and TGWE modeling, as shown in the Table 4, may significantly affect gravitational sounding of internal flow structure and strength (Zhang KK et al., 2017a;Kong DL et al., 2018b).

Discussions and Conclusions
In this article, we have directly and clearly demonstrated that TWE and TGWE models can yield large numerical differences in gravitational distortions of zonal flows. One of our illustrative examples of zonal circulation is equatorially symmetric; the other is equatorially anti-symmetric. The important point is that both of these flows are physically realistic in the sense they can indeed be produced and sustained by thermal convections inside a rapidly rotating planet. The theoretical analyses have predicted that the two terms and in TGWE are generally of similar size, if the zonal circulation is not artificially assumed to be highly concentrated in the outer shallow layer of a planet. Our numerical comparisons have verified this fact, as listed and discussed in Tables 3 and 4.
There are two circumstances under which the term might be less important. The first is the case of a terrestrial planet whose atmospheric self-gravity is known to be much smaller than the gravity from the rocky planet. The second is the case that, after the TG-WE is solved for small-scale dominant zonal flows, one finds that the contribution of the ρ 0 g′ term to higher order is small because of positive-negative quadrature cancellation in the integral 22. However, in studies of gaseous planets, unlike terrestrial planets, self-gravity of the bulk of gas is crucially important and of direct concern to the research. Also, there is no a priori knowledge that deep zonal flows in Jupiter and Saturn are extended from cloud-top winds. That zonal circulations might exist in a vast volume of low-conductivity interior cannot be dismissed. So, generally, there is no guarantee that the term would be much smaller than the term, under the assumption that non-spherical effects of solid-body rotation are neglected. ρ 0 g ′ Technically speaking, it is mathematically invalid to "prove" that a term is small through solving the equation without the term. More specifically, for zonal-flow-distorted gaseous planets, the relative importance of the term can be judged only after the more complete TGWE, rather than just the TWE, is computed. The balance resulting from TWE is simply a different irrelevant physics from what is described by TGWE.
A more intuitive example may help explain the idea. In fluid dynamics, the Euler equation neglects the viscosity term in the Navier-Stokes equation. Euler's equation therefore cannot be employed to study any viscous effects, such as the problems of boundary layer and thermal dissipations. It should be obviously wrong to use a solution of the Euler equation to claim that a boundary layer effect is small.
Deep structure of zonal flows in fast-rotating gaseous planets is still a puzzle. Gravitational sounding itself is subject to the nonuniqueness problem, and to uncertainties related to the EOS, conductivity, and so on. In the future, more independent measurements are needed -such as seismology, microwave sounding, data on magnetic field secular variation, or even in-situ entry probing. Future study of this problem poses both challenges and opportunities.  Table 2. The fitting of the symmetric zonal flow speed as a polynomial of s = r sin θ/R, up to the 12th order. The speed curve shown in Figure 1 can be recovered by the summation .  Table 3. The even-order zonal gravitational coefficients that are induced by the equatorially symmetric zonal flow. Contrast between results computed by TGWE and TWE is clearly seen. There is a large discrepancy in lower order coefficients and , which is expected from theoretical order-of-magnitude analysis.