Depth Dependent Dynamics Explain the Equatorial Jet Difference Between Jupiter and Saturn

Jupiter's equatorial eastward zonal flows reach wind velocities of ∼100 m s−1, while on Saturn they are three times as strong and extend about twice as wide in latitude, despite the two planets being overall dynamically similar. Recent gravity measurements obtained by the Juno and Cassini spacecraft uncovered that the depth of zonal flows on Saturn is about three times greater than on Jupiter. Here we show, using 3D deep convection simulations, that the atmospheric depth is the determining factor controlling both the strength and latitudinal extent of the equatorial zonal flows, consistent with the measurements for both planets. We show that the atmospheric depth is proportional to the convectively driven eddy momentum flux, which controls the strength of the zonal flows. These results provide a mechanistic explanation for the observed differences in the equatorial regions of Jupiter and Saturn, and offer new understandings about the dynamics of gas giants beyond the Solar System.


Introduction
Measurements from recent decades have revealed several key similarities between the atmospheres of the two gas planets of the Solar System, Jupiter and Saturn.They are both dominated by east-west jet-streams (zonal flows) surrounding the planet (Garcia-Melendo et al., 2011;Tollefson et al., 2017), which have been fairly stable over the past 50 years, since the Voyager era (Fletcher et al., 2020).The jets on both planets penetrate thousands of kilometers into the interior, where the pressure exceeds 10 5 bar (Galanti et al., 2019;Kaspi et al., 2018), to a depth where the Lorentz force might become significant (Kaspi et al., 2020;Liu et al., 2008).Both planets have wide regions possessing strong turbulent activity, like the famous "storm alley" (Sromovsky et al., 2018) or Great White Spot (Sánchez-Lavega, Fischer, et al., 2019) on Saturn, or the Great Red Spot (Parisi et al., 2021;Wong et al., 2021) and ovals (Sánchez-Lavega et al., 2001) on Jupiter.The turbulence at the cloud level was shown to be related to the atmospheric mean flow, on both planets, hence the midlatitude jets are considered to be eddy-driven (Del Genio & Barbara, 2012;Duer et al., 2021Duer et al., , 2023;;Salyk et al., 2006).
The dynamical similarities between the planets are not surprising, as they are both of similar size (Hubbard & Smoluchowski, 1973), have an overall similar composition (Encrenaz et al., 1999), rotate with nearly the same rotation rate (Helled et al., 2015) and have comparable internal heat sources (Guillot et al., 2022), all of which affect the planetary atmospheric dynamics.However, alongside the similarities, the two planets possess notable differences.Some of these differences have comprehensive explanations and theories, like the cloud optical depth and condensation level (Atreya et al., 1999;West et al., 2004), the temperature vertical structure (Mueller-Wodarg et al., 2008), the penetration depth of the zonal winds, being inversely proportional to the planetary mass (Kaspi et al., 2020), and others.Nevertheless, there remain certain unexplained disparities, with one of the most prominent being the distinct differences observed in the equatorial zonal flows on both planets.planetary surface to give the cloud-level latitude of the tangent cylinder (Figures 1a and 1b, dashed white lines).The tangent cylinder is an imaginary cylinder, parallel to the planet's rotation axis, at a specific equatorial depth, which differentiates the equatorial region from the midlatitudes (Figures 1c and 1d, dashed black lines).At the equatorial region, corresponding to outside from the tangent cylinder, both planets have a strong prograde flow (in the direction of rotation), which is superrotating within ∼6°(∼10°) latitude of the equator (Imamura et al., 2020) and reaching ∼100 m s 1 (∼300 m s 1 ) on Jupiter (Saturn) (Garcia-Melendo et al., 2011;Tollefson et al., 2017) (Figures 1a and 1b).Note that the wind profile of Saturn is shown according to the latest estimate for the rotation rate (Mankovich et al., 2023).This estimate aligns with findings from the past two decades, indicating Saturn's rotation period to be approximately between 10 hr 32 min and 10 hr 34 min (Anderson & Schubert, 2007;Helled et al., 2015;Mankovich et al., 2019Mankovich et al., , 2023;;Read et al., 2009).Closer to the tangent cylinder there are strong retrograde flows reaching ∼30 m s 1 on Jupiter (Tollefson et al., 2017) and ∼90 m s 1 on Saturn (Garcia-Melendo et al., 2011) (Figures 1a and 1b).Notably, the equatorial winds on Saturn are stronger and latitudinally wider than on Jupiter (Figure 1).
The extent and strength of the equatorial zonal flows on giant planets are influenced by several factors, which have been studied using numerical simulations and theoretical arguments (Sánchez-Lavega, Sromovsky, et al., 2019;Vasavada & Showman, 2005).The planetary rotation rate directly impacts the direction and strength of the equatorial zonal flows, with rapidly rotating planets having stronger prograde zonal winds (Camisassa & Featherstone, 2022;Gastine et al., 2013;Kaspi et al., 2009).Heating sources and the thermal diffusivity control the convection and the onset of convection, which in turn affects the energy transferred to the mean zonal wind  (Garcia-Melendo et al., 2011;Tollefson et al., 2017), along with the latitude of the tangent cylinder (dashed white lines).(c, d) The zonal winds of Jupiter and Saturn projected onto a sphere with the radial decay profiles from the Juno (Kaspi et al., 2018) and Cassini (Galanti et al., 2019) results, respectively.Jupiter's inner surface is at 0.955 R J (3,000 km) and Saturn's at 0.84 R S (9,000 km), representing the depth of the two atmospheres according to the gravity analysis.The latitudinal slice is from 90°N to 20°S on Jupiter and 40°S on Saturn.The colorbar is shared between panels (c, d) emphasizing that Saturn's zonal jets are stronger than Jupiter's.The two planets are to scale.(Busse, 1976;Christensen, 2001Christensen, , 2002;;Kaspi et al., 2009).The depth of the convective layer has been shown to influence the position of the tangent cylinder and hence the latitudinal extent of equatorial dynamics (Busse, 1983;Gastine et al., 2014;Heimpel & Aurnou, 2007;Heimpel & Gómez Pérez, 2011).Numerical simulations have also shown that the numerical viscosity parameter (Kaspi et al., 2009;Showman et al., 2011), boundary conditions (Jones et al., 2003), and factors such as background reference entropy state and compressibility (and the possible presence of a stable layer) (e.g., Gastine & Wicht, 2012, 2021;Heimpel et al., 2022;O'Neill et al., 2017;Wicht & Gastine, 2020;Wulff et al., 2022) affect the strength of the jets in simulations, including the equatorial region, while magnetic field effects may also play a role in the structure of zonal winds (Gastine et al., 2014;Yadav et al., 2020Yadav et al., , 2022)).Nonetheless, these mechanisms cannot fully explain the differences in the equatorial zonal flows between Jupiter and Saturn.Notably, Jupiter exhibits a faster rotation rate and a stronger internal heat source compared to Saturn.While these characteristics are anticipated to enhance cylindrically oriented convection, driving momentum through upgradient fluxes and resulting in the generation of a stronger equatorial zonal mean flow (e.g., Kaspi et al., 2009), the observed flow strength appears to be weaker on Jupiter.Furthermore, it is likely that thermal diffusivity, viscosity, and boundary conditions are fairly consistent between the two planets, and as such, these factors alone cannot account for the observed differences.
In this study, using the recent estimations of the different atmospheric depth of Jupiter and Saturn (Galanti et al., 2019;Kaspi et al., 2018), we investigate the cause of the difference in equatorial zonal flow between Jupiter and Saturn.Utilizing high-resolution numerical simulations, we explore the relationship between the depth of the atmosphere and the formation of these flows.By manipulating the domain size, while keeping other simulation parameters constant, we are able to provide an explanation for this phenomenon.Additionally, we delve into the underlying mechanism responsible for the formation of equatorial zonal flows on gas giants in general.
The hydrodynamic set of equations can be written with dimensionless control parameters: the modified Rayleigh number Ra * = g o ΔS c p Ω 2 L , the Ekman number Ek = ν ΩL 2 , the Prandtl number Pr = ν κ , and the dissipation number Di = g o L c p T o = η e N ρ /n 1) , in addition to the polytropic index n (see Supporting Information S1), the number of scale heights of density across the shell N ρ = ln( ρ i ρ o ) and the radius ratio η = r i r o , where the subscripts i and o denote the inner and outer boundaries, respectively (e.g., Jones et al., 2011).In the Rayleigh number, g o is the gravitational acceleration at the top of the domain, ΔS is the entropy difference across the domain, c p is the specific heat capacity at constant pressure, Ω is the planetary rotation rate, and L is a typical length unit.The modified Rayleigh number can also be formulated by incorporating the classical Rayleigh number (Ra = g o ΔSL 3 c p κν ) , such that Ra * = Ek 2 Pr Ra.In the Ekman number, ν is the kinematic viscosity and in the Prandtl number κ is the thermal diffusivity constant, ρ is the density and r is radius.For the non-dimensional equations we adopt the following units: length → L (shell depth), time → 1/Ω, temperature → T o , density → ρ o , and entropy → ΔS.
Using the parameters defined above we can write the non-dimensional momentum and thermodynamic equations (e.g., Heimpel et al., 2016): Geophysical Research Letters (2) In the momentum equation (Equation 1), u is the 3D velocity vector (u in the zonal direction, v in the meridional direction, and w in the radial direction), t is time, ẑ is the vertical coordinate (parallel to the axis of rotation), r is the radial coordinate, p is pressure, ρ is the background density, S in entropy, and D is the viscous stress tensor, , and e ij = 1 2 ( ∂u i ∂x j + ∂u j ∂x i ) .In the thermodynamic equation (Equation 2), T is the background temperature profile, Q i is a radially dependent internal heating (Jones et al., 2011), and Π is the viscous heating term, where Π = 2 ρ e ij e ji 1 3 (∇ ⋅ u) 2 ) .With the addition of the continuity equation and an equation of state (see Supporting Information S1) the system of equations is complete.
Within the scope of this study, we refrain from delving into the determination of the mechanism governing atmospheric depth, factors such as Ohmic dissipation, stable layers, and significant density variations have been proposed (e.g., Cao & Stevenson, 2017;Christensen et al., 2020;Kaspi et al., 2009;Liu et al., 2008;O'Neill et al., 2017;Wulff et al., 2022).Instead, our focus is directed toward examining the impact of this depth on the resulting dynamics.Hence, we do not include magnetic field equations and their impact on the flow fields, a unique background entropy profile or a complicated equation of state.Our experiments consists of a series of simulations, where we vary the domain depth (representing the "atmospheric" non-conductive region).In the more shallow setups we also vary the viscosity term (Ekman number) to allow numerical stability.As all numerical simulations solving for gas giants, our simulations are overforced and include numerical viscosity that is orders of magnitude larger then the molecular viscosity of Jupiter and Saturn to allow high Rayleigh numbers (Showman et al., 2011).Our goal is to keep the main control parameters of the model fixed, hence we adjust the energy source, the thermal diffusivity constant and the kinematic viscosity coefficient according to the chosen domain depth (L).The chosen set of parameters is well-within the customary values used in the benchmarks (Jones et al., 2011).For values of the parameterization see Tabs.S1 and S2 in Supporting Information S1.The simulations are calculated until a statistical steady state is reached (about 1,000 rotations), and all results are shown for either a snapshot or time-averaging over 300 rotations in steady state.The main three control parameters that are kept constant between the different simulations are the modified Rayleigh number, the Ekman number, and the fluid Prandtl number.As we are interested in the statistical steady state solution, which is assumed to be the current dynamical state for both Jupiter and Saturn, the results are presented long after the onset problem, which is influenced by the domain depth (Barik et al., 2023;Dormy et al., 2004).For consistency, all of the simulations run with Jupiter's radius (R) and rotation rate.

Eddies Outside the Tangent Cylinder
Convection in fast-rotating spherical bodies tends to align in cylinders parallel to the rotation axis.These cylinders are tilted in the direction of rotation and, hence, transfer positive momentum outward, driving the equatorial zonal flows (e.g., Busse, 1976Busse, , 1982Busse, , 2002;;Zhang & Busse, 1987).To investigate this mechanism, we focus on the convergence of eddy fluxes oriented perpendicular to the axis of rotation within the entire region outside the tangent cylinder (Figure 2).Examining the zonally averaged zonal wind (Figure 2b) along with the convergence of the zonally averaged perpendicular eddy momentum flux (Equation S7 in Supporting Information S1, Figure 2c) in an idealized simulation (r min = 0.84 R, Ek = 5 × 10 −4 , Pr = 2, and Ra* = 0.0132), reveals that the zonal wind is being driven by these fluxes, which are aligned in cylinders throughout the entire domain (except at the outer boundary near the equator and near the tangent cylinder).Thus, the leading momentum source in the zonal momentum equation can be described by where v′ ⊥ = w ′ cos θ + v ′ sin θ and θ is the latitude (see Supporting Information S1 for the full equation in an anelastic form in spherical coordinates, Equations S6 and S7 in Supporting Information S1).The bar represents a zonal average and the prime denotes deviations from that average.The zonal momentum equation describing such dynamics includes additional forces, like the Coriolis force, the mean eddy fluxes and the viscosity term, none of which dictate the direction of the zonal wind (the Coriolis force cancels itself in the equatorial region, the mean eddy fluxes are small due to the flow being geostrophic (Galanti et al., 2017) and the viscosity terms acts against the reminder, (Duer et al., 2023;Kaspi et al., 2009), see Supporting Information S1).
Examining the instantaneous zonal wind reveals the structure and number of the convection columns throughout the simulation (Figure 2d).As this simulation is calculated with a relative high Ekman number (Ek = 5 × 10 4 ), the resulting columns are near ideal and spread from the domain's inner boundary to the domain's outer boundary (this can be seen in the eddy velocity field, Figure 2e).A more turbulent simulation for a domain with the same depth, in which the eddies are stronger but less organized, and they do not span the entire depth of the domain, is shown in Figure S1 in Supporting Information S1.Additionally, the zonally averaged convergence term changes sign at multiple depths.Yet, averaging the simulation over time shows that the momentum source keeps its sign and strength near the tangent cylinder (inner boundary) and near the equator (outer boundary), as these regions allow only tilted columns that diverge and converge, respectively.The middle region averages to small values, as convection columns are being built and destroyed while the simulation continues (Figure S1c in Supporting Information S1).
The columnar structure of the flow is illustrated in Figure 2a, showing the transition between the divergence and convergence of the eddy momentum fluxes together with the zonal wind direction dictated by it.As the divergence of the eddies occurs outside the tangent cylinder, the retrograde flanks of the zonal flow must also be positioned outside the tangent cylinder.This implies that the atmospheric depth needs to match the depth of the retrograde flow adjacent to the prograde equatorial jet in superrotating giants.

The Extent of the Equatorial Zonal Flows
Next, we examine how the depth of the atmospheric flow affects the extent of the equatorial zonal winds.The steady-state, zonally averaged, zonal wind field is compared between 10 different simulations, calculated with different domain depths (Figure 3).All the simulations in the top row (Figures 3a-3e) are calculated with identical control parameters: Ek = 5 × 10 5 , Pr = 2, and Ra* = 0.0132 and a scaled reference state (see Supporting Information S1), while adjusting the other system constants (see Supporting Information S1 for details).For domains smaller than r min = 0.85 R the Ekman number must be larger to constrain the eddy-related phenomena inside the tangent cylinder (poleward from the tangent cylinder), allowing to reach a dynamical steady state.Hence, the bottom five simulations (Figures 3f-3j) are calculated with a larger Ekman number, and their control parameters are: Ek = 5 × 10 4 , Pr = 2, and Ra* = 0.0132.For comparison, the domain of r min = 0.84 R is presented both with both Ekman numbers.
The latitudinal extent of the equatorial zonal wind field is compared between the different simulations, and also to Jupiter and Saturn (Figure 3k).It is defined according to the projection of atmospheric region on the outer radius (the tangent cylinder), represented by the latitude α.In the simulations, it can simply be set according to the shell's depth (L = R r min ) and the chosen planetary radius (R), such that α = arccos 1 L R ) .For Jupiter and Saturn, we set L according to the decay profiles provided from the gravity analysis, 3,000 km for Jupiter and 9,000 km for Saturn.This gives values of α Jup = 17°and α Sat = 32°(dashed white lines, Figures 1a and 1b).The latitude α separates the prograde equatorial zonal wind and partly the retrograde flow engulfing it from the alternating jets at midlatitudes (Figures 1c and 1d).Specifically, it nearly coincides with the maximal retrograde flow, on both planets (Figures 1a and 1b).For comparison with the simulations, we define the latitude β, which is the latitude where the retrograde jet adjacent to the equatorial prograde flow peaks at the outer boundary.This determination represents well the extent of the equatorial dynamics, since tilted convection columns are the momentum source of the equatorial zonal jets as detailed in the previous section.For asymmetric simulations (or observations), an averaged value was taken between the hemispheres.This gives values of β Jup = 16.6°andβ Sat = 34.1°.For both planets the values of α and β are nearly equal.Presenting them together with the 10 simulations of Figure 3, it is apparent that all 10 simulations also exhibit near equal values of α and β, with a R 2 value of 0.99.The robustness of our findings is assessed by varying the main parameters to ensure the general applicability of the results.Particularly, we conducted experiments by scaling the values of Pr and Ra* across five different domain depths, and reexamined the correlation depicted in Figure 3k for the additional 35 simulations.The results prove to be robust under these variations, as illustrated in Figure S2 in Supporting Information S1 where the values of α and β maintain the 1:1 ratio.This emphasizes the consistent connection between domain depth and the ensuing equatorial dynamics.
Figure 3 reveals the dominance of the domain depth in controlling the latitudinal extent of the equatorial zonal flows and the adjacent retrograde jets.It also shows, that the two retrograde jets adjacent to the prograde jet have the same origin and are linked outside the tangent cylinder.This relationship has been presented in previous studies (e.g., Busse, 1983;Christensen, 2001;Gastine et al., 2014;Heimpel & Aurnou, 2007), however, the comparison with observations for Jupiter and Saturn, which only became recently available, constrain the position of the tangent cylinder and allow fixating on the relative zonal winds region that participates in the equatorial dynamics.Note that the latitude of the tangent cylinder (α) of the deepest simulation reaches nearly 70°latitude, leaving no place for midlatitude dynamics.This could explain the absence of midlatitude jets, for example, in deep, convection-driven, simulations (e.g., Kaspi et al., 2009).The columnar structure of the deep flow, characterized by Taylor columns, is uninterrupted in the equatorial region and the adjacent retrograde jet up to its maximal absolute value.This means that the equatorial jet and its adjacent jets are close to north-south symmetric and continue throughout the planet's interior, while the density varies at different depths.Electrical conductivity, of course, changes with depth as well (French et al., 2012;Liu et al., 2008), and might be directly linked to what sets the depth participating in the equatorial dynamics (Kaspi et al., 2020).However, even without the magnetic field constraints, the agreement between the simulations and recent estimations for Jupiter and Saturn implies that the latitudinal extent of the equatorial zonal flows is dictated directly by the atmospheric depth.

The Strength of the Zonal Flows and the Eddies Driving It
Lastly, we examine the strength of the equatorial zonal winds.It is evident that in shallower domains the zonal winds weaken (Figure 3).As the main three control parameters in the simulations are kept constant (besides the Ekman number, which is increased for simulations shallower than r min = 0.84 R), the zonal wind velocity must be dependent on other factors that are different between the simulations.We examine eight simulations together on the same panels (Figure 4), five deep simulations (Figures 3a-3e) and three shallow simulations with higher Ekman number (Figures 3f,3h,and 3j).As before, the simulations with r min = 0.84 R, are presented with both Ekman values.The zonal winds of the high Ekman simulations are significantly weaker than those of the deep simulations, reaching wind velocities of a few meters per second at the equator.This is an expected outcome as the viscosity was increased by an order of magnitude, the zonal winds should be reduced by approximately the same ratio (the viscosity acts against the zonal wind acceleration term, see Equation 1 and Supporting Information S1).To present all eight simulations together, the three simulations with the high Ekman value are scaled (the velocities are multiplied by a factor of 6; dashed lines, Figure 4).
The wind velocity appears to be proportional to the simulation depth at all latitudes outside the tangent cylinder (Figure 4a) and at all simulation depths (Figure 4b).This suggests that the strength of the zonal winds must be directly influenced by the size of the domain (more precisely, its volume, see Supporting Information S1), allowing more momentum to be transferred to the acceleration term of the zonal wind (Equation 1).In fastrotating, convection-driven simulations, such as the simulations presented here, the mechanism driving momentum toward the equatorial zonal wind is the tilted convection columns that originate outside the tangent cylinder, as shown in Figure 2. At the equatorial plane, this can be represented by the radial eddy momentum flux (u′w′, Figure 4c).
In idealized simulations, the columns are spread almost evenly through the entire plane between the boundaries (e.g., Figures 2d and 2e).The tilt of the columns is in the prograde direction, hence, the eddy momentum flux is always positive (see illustration, Figure 2a).In more turbulent simulations, the columns are noisy and continue to break and reappear (Figures S1d and S1e in Supporting Information S1).Examining the eddy momentum flux at the equatorial plane reveals that indeed this term is always positive, and while in the more turbulent simulations it is more noisy (even when looking at time-averaged values), the convection columns are positioned close to the tangent cylinder and are dominant between the boundaries (Figure 4c).The eddy momentum fluxes are also proportional to the depth of the simulation.Deeper simulations allow stronger eddies to evolve, thus, giving more momentum to the zonal wind, resulting in higher velocities.This explains why Saturn's equatorial winds (both prograde and retrograde) are stronger than Jupiter's, as the tangent cylinder is positioned three times deeper (Figure 1).However, due to lack of measurements of eddy fluxes beneath the cloud level, comparison to the real atmospheres is unavailable.
The transition in the zonal wind direction, from retrograde to prograde (Figure 4b), should be positioned where the eddy momentum flux term changes its tendency, as the full term that appears in the momentum equation is the convergence of the eddy fluxes (Equation 3).To examine this, we calculate the correlation between the location where the zonal wind changes its direction ( ū = 0) and the location where the eddy momentum flux is maximal (equivalent to zero convergence).These values in all the simulations are very close to the 1:1 ratio, with r 2 value of 0.99 (Figure 4d), suggesting that the zonal wind field strength is directly related to the eddy momentum fluxes at the equatorial plane.While a quantitative comparison with Jupiter and Saturn is unavailable, given that we do not have measurements of eddy momentum fluxes at depth, they are projected on the 1:1 ratio according to the abscissa (Figure 4d), which is found according to the flow decay profiles obtained by the gravity measurements by Juno (Kaspi et al., 2018) and Cassini (Galanti et al., 2019).This relationship implies on the deep structure of eddies in the Jovian and Saturnian atmospheres.where the vertical eddy momentum flux is maximal (ordinate) and the depth where the zonal wind crosses the zero velocity (abscissa).The dashed line represent the 1:1 ratio.The r 2 value of the eight simulations is 0.99.Also shown are the abscissa values for Jupiter (Kaspi et al., 2018) and Saturn (Galanti et al., 2019) (gray), according to the flow decay profiles obtained by the gravity measurements, projected on the 1:1 ratio.All simulations are at a statistical steady state and the terms shown are time-averaged.The legend is shared between the panels.

Conclusions
In this study, we show that the penetration depth of the atmospheric dynamics directly influences the extent and strength of the cloud-level equatorial zonal flows on Jupiter and Saturn.The latitudinal extent of the equatorial zonal flows, defined here as the latitude of the maximum retrograde flows engulfing the prograde equatorial zonal wind, is dictated by the depth of the atmospheric dynamics and is compatible with the cloud-level wind structure and the zonal jets' depth given by the gravity measurements for Jupiter and Saturn.The depth of the atmospheres on these planets is inversely proportional to the mass of the planets.Jupiter's mass is approximately three times heavier than Saturn's, leading to denser gas closer to the outer radius of the planet (Guillot, 1999).The planets' mass affects the depth at which the gas becomes ionized and the potential depth for a stable layer to exist, hence, influencing the extent of the planetary atmospheric depth regardless of what mechanism might dissipate the zonal jets (Christensen et al., 2020;Gastine & Wicht, 2021;Liu & Schneider, 2010;Liu et al., 2008;Wulff et al., 2022).
The strength of the equatorial zonal flows might depend on parameters such as rotation rate, viscosity, internal heating, and boundary conditions.Yet, with a representative set of control parameters we show that the depth of the tangent cylinder is proportional to the strength of the eddy momentum flux, and hence to the strength of the zonal flow.This suggests that the difference between the equatorial zonal flows of Jupiter and Saturn could be a direct result of their atmospheric depth, mainly due to the similarities in the other parameters that might affect the zonal flow strength.The relation between Jupiter's and Saturn's equatorial zonal wind velocity is 1:3, which is the same ratio as their atmospheric depth, and nearly the same ratio as the resulting volumes outside the tangent cylinder (see Supporting Information S1).Although the 1:3 ratio for Jupiter's zonal wind strength compared to Saturn's might subtly vary depending on rotation rate and specific measurements (Sánchez-Lavega et al., 2023), all wind profiles consistently support the notion that Saturn's winds are demonstrably stronger and deeper than Jupiter's.These results cannot apply to Uranus and Neptune, as they are not superrotating and possess much shallower atmospheres (Kaspi et al., 2013;Soyuer et al., 2020), which might not even hold such convection columns.
We also show that the zonal flows outside the tangent cylinder are aligned in cylinders, along with the eddy momentum fluxes that drive them.Deeper atmospheres allow stronger eddies in the direction perpendicular to the axis of rotation, transferring more momentum to the zonal flows.This mechanism is restricted to the region outside the tangent cylinder, and cannot explain the midlatitudes and polar dynamics.The compatibility of these modeling results with the winds on Jupiter and Saturn suggests that tilted Taylor convection columns control the dynamics in gas giant atmospheres outside the tangent cylinder, as suggested decades ago (Busse, 1976).Finally, based on the results presented here, which are general for convectively driven giant planets, future estimation of the cloud-level zonal winds of superrotating exoplanets may allow to constrain their radial atmospheric depth, and hence, better understand their interior density structure and evolution.

Figure 1 .
Figure 1.The zonal winds on Jupiter and Saturn.(a) Jupiter's and (b) Saturn's zonally averaged cloud level zonal wind profiles(Garcia-Melendo et al., 2011;Tollefson et al., 2017), along with the latitude of the tangent cylinder (dashed white lines).(c, d) The zonal winds of Jupiter and Saturn projected onto a sphere with the radial decay profiles from the Juno(Kaspi et al., 2018) and Cassini(Galanti et al., 2019) results, respectively.Jupiter's inner surface is at 0.955 R J (3,000 km) and Saturn's at 0.84 R S (9,000 km), representing the depth of the two atmospheres according to the gravity analysis.The latitudinal slice is from 90°N to 20°S on Jupiter and 40°S on Saturn.The colorbar is shared between panels (c, d) emphasizing that Saturn's zonal jets are stronger than Jupiter's.The two planets are to scale.

Figure 2 .
Figure 2. Snapshot of the dynamics in an ideal simulation with r min = 0.84 R. (a) An illustration of a tilted convection column (the tilt is in the positive zonal direction, i. e., into the page, see inset), pushing positive momentum in a direction perpendicular to the axis of rotation, where Ω is the rotation axis direction and r⊥ is the perpendicular to the rotation axis direction.The dashed black line is aligned with the zero values in panels (b, c).Due to the tilt, the perpendicular momentum flux is always positive (u ′ v′ ⊥ > 0) inside the columns (see inset, where φ being the zonal direction).(b) The zonally averaged zonal wind (m s 1 ) and (c) the eddy momentum flux convergence (Equation S7 in Supporting Information S1, m s 2 ) that is perpendicular to the axis of rotation, outside the tangent cylinder (the abscissa range is [0.84-1] r/R).It is evident that the zonal wind is proportional to the eddy momentum flux convergence.(d, e) A snapshot demonstrating the zonal wind and the eddy velocity at the equatorial plane (m s 1 ), respectively.An equivalent figure for a more turbulent simulation with the same domain depth is provided in Figure S1 in Supporting Information S1.

Figure 3 .
Figure 3. Zonally averaged zonal wind (colors, m s 1 ) of 5 simulations dominated by a prograde equatorial jet and an adjacent retrograde jet, each extend to a different depth (r min ).All the simulations are calculated with identical control parameters: Ek = 5 × 10 5 , Pr = 2, and Ra* = 0.0132.(f-j) Shallower 5 simulations with equal Pr and Ra*, but with Ek = 5 × 10 4 .(k) β (ordinate) and α (abscissa), representing the latitude where the retrograde jet peaks in each simulation and the projection of r min (the tangent cylinder) on the outer boundary, respectively.The 10 simulations are presented (top-blue, bottom-cyan) along with the values for β and α of Jupiter (green) and Saturn (red).The r 2 value of panel (k) is 0.99.The dashed line represents the 1:1 ratio.

Figure 4 .
Figure 4. (a) The zonally averaged zonal wind as a function of latitude at the outer radius (r = R) of the five top simulations (Figures 3a-3e, solid lines) and three additional simulations with higher Ekman number (Figures 3f, 3h, and 3j, dashed lines).The three simulations with high Ekman number are scaled (the velocities are multiplied by 6) such that the eight simulations are presented in the same plot.All eight simulations have identical Rayleigh and Prandtl numbers.(b) The zonally averaged zonal wind as a function of normalized radius at the equator (θ = 0°) of the same eight simulations.(c) The vertical eddy momentum flux u′w′) with normalized radius at the equator of the eight simulations.(d) Comparison between the depthwhere the vertical eddy momentum flux is maximal (ordinate) and the depth where the zonal wind crosses the zero velocity (abscissa).The dashed line represent the 1:1 ratio.The r 2 value of the eight simulations is 0.99.Also shown are the abscissa values for Jupiter(Kaspi et al., 2018) and Saturn(Galanti et al., 2019) (gray), according to the flow decay profiles obtained by the gravity measurements, projected on the 1:1 ratio.All simulations are at a statistical steady state and the terms shown are time-averaged.The legend is shared between the panels.