Deep, Closely-Packed, Long-Lived Cyclones on Jupiter's Poles

Juno Mission to Jupiter has found closely-packed cyclones at the planet's two poles. The observation that these cyclones coexist in very confined space, with outer rims almost touching each other but without merging, poses a big puzzle. In this work, we present numerical calculations showing that convectively sustained, closely-packed cyclones can form and survive without merging for a very long time in polar region of a deep rotating convection zone (for thousands of planetary rotation periods). Through an idealized application of the inertial stability criterion for axisymmetric circulations, it is found that the large Coriolis parameter near the pole plays a crucial role in allowing the cyclones to be packed closely.

ranging from 5600 km to 7000 km, but recent observations have found that the number of surrounding circumpolar cyclones increased to six and then changed back to five again (Adriani et al. 2020). Aside from the presence of small precessions around the poles, the cyclones are quite stationary. Measured at 1500 km away from their respective centers, the cyclonic rotation periods range from 27.5 to 60 hours . Stationary cyclonic vortices have also been observed on Saturn's polar regions, but only one strong cyclone dominates at each pole (Godfrey 1988;Vasavada et al. 2006).
It is well known that beta effect associated with rotation of a sphere can cause cyclones to migrate poleward and anticyclones to migrate equatorward (Schecter & Dubin 1999). Accumulation of cyclonic vorticity in polar region has been demonstrated by potential vorticity calculation based on contour dynamics (Scott 2011). Idealized aqua-planet general circulation model experiments have shown that tropical cyclones tend to cluster toward the pole (Merlis et al. 2016). For giant planets, polar cyclonic vortices have been examined by a shallow-water model in the setting of γ-plane (O'Neill et al. 2015;Nof 1990;O'Neill et al. 2016) (the Coriolis parameter is a quadratic function of distance to the pole, also called polar β-plane). With Saturn's parameters, the model, driven by hypothesized moist convection in the surface weather layer, created a strong polar cyclonic vortex. For Jupiter, the model predicted random vortices near the pole (O'Neill et al. 2016). Using a different shallow γ-plane model, Brueshaber et al. (2019) investigated the effects of randomly injected storms (local mass perturbations) on polar flows. They found that a large cyclonic polar vortex is usually formed when the Burger number (the square of the ratio of Rossby deformation radius to planetary radius) is large. Multiple small vortices can be formed when this number is small. However, the simulated cyclones do not remain in the pole for significant length of time.
In this paper we consider a dynamically self-consistent approach based on solving the compressible flow equations for a deep rotating convection zone (CZ). Juno's gravitational measurements indicate that Jupiter's zonal winds extend 3000 kilometers below the visible surface at low latitudes, and 0 − 500 kilometers at high latitudes (Kaspi et al. 2018). Despite the high uncertainty in the polar region, it does not contradict with the presence of a deep CZ. Numerical simulations (Chan 2007;Käpylä et al. 2011;Chan & Mayr 2013;Guervilly et al. 2014;Stellmach et al. 2014;Favier et al. 2014;Aurnou et al. 2015;Cai 2016;Heimpel et al. 2016;Guervilly & Hughes 2017;Yadav & Bloxham 2020) have demonstrated that large-scale, long-lived vortices (cyclonic or anticyclonic) can be generated spontaneously in a fast rotating CZ. Our numerical model, to be discussed in Section 2, is along this pathway.
Simulation of vortices in a shearing zonal flow indicated that two stable vortices tend to merge if their cross-zonal distance is smaller than a critical value approximately equal to the vortices' radii (Marcus 1990). Two same-sign vortices generally merge when the separation is below a certain threshold value (also of the order of the vortex radius) (Cerretelli & Williamson 2003;Folz & Nomura 2017). A principal puzzle, therefore, is why Jupiter's closely-packed polar cyclonic vortices do not merge. In this study, we provide numerical evidence that convection generated long-lived cyclones in the polar γ-box (three-dimensional) can be packed without merging for a long time, and an idealistic application of the inertial stability criterion (Yanai 1964;Schubert & Hack 1982;Kloosterziel et al. 2007) is invoked to provide some insights (Section 3).
While considering vortex packing in the polar regions of Jupiter as deep convective phenomenon, one cannot avoid the question of why Saturn displays only a single cyclone in each pole. We address this question in Section 4.

THE MODEL
Numerical simulations are performed by solving the fully compressible hydrodynamic equations for an ideal gas in a rectangular box:

Cai, Chan, & Mayr
where ρ is the density, v is the velocity, p is the pressure, Σ is the viscous stress tensor, g is the gravitational acceleration, Ω is the angular velocity, E = e + ρv 2 /2 is the total energy density, e is the internal energy, F d is the diffusive energy flux, T is the temperature, c p is the heat capacity under constant pressure, τ is the time scale of Newton cooling (to mimic fast radiative loss above the planet's optical surface), the subscript top denotes the value at the top of the box. An ideal gas law p = ρR s T , where R s is the specific gas constant, is used for the equation of state. The ratio of specific heats of the gas is γ ad = 1.47 (Horedt 2004). We adopt a 'large eddy simulation' approach in which sub-grid-scale (SGS) turbulent processes in the momentum equation are modeled by a SGS kinematic viscosity ν according to the Smagorinsky scheme (Smagorinsky 1963) where ∆x, ∆z are local sizes of the horizontal and vertical grids, σ is the strain rate tensor, the double dot denotes tensor contraction, and the coefficient c ν is chosen to be 0.28.
The hydrodynamic equations (1-3) are solved by an explicit finite-difference method. The numerical code is written in a way that conserves total mass, total energy, and momentum to round-off limits.
The computational domain is discretized by a staggered grid. Pressure, density, and temperature are located on grid levels, while velocities are located on half-grid levels. In the original version of this code (Chan & Sofia 1986), the equations were integrated by a semi-implicit method that can accommodate time steps a few times larger than the Courant-Friedrichs-Lewy restriction. For the simulations performed in this paper, the Mach number is on the order of O(0.1), and the speed-up achieved by the semi-implicit treatment of time integration is not advantageous versus efficiency of parallelization. For this reason, we instead use an explicit predictor-corrector method.
The three-dimensional numerical model, apart from the location-dependent Coriolis terms and some differences in parameters, is essentially the same as the one used in Chan & Mayr (2013).
The rectangular 'γ-box' is considered to represent a piece of spherical shell with polar axis through its center, along the vertical direction. The local gravitational acceleration g defines the vertical direction (z) for all points in the box (it is a constant in the present model). The components of the Coriolis term vary with location as the rotation vector Ω can be tilted with respect to g. The inclined angle between Ω and −g is the colatitude θ, and the latitude is φ = 90 o − θ. The form of the vector equations remain unchanged, but Ω now depends on location. As the box is threedimensional, the effects of both the vertical Coriolis parameter f = 2Ω sin φ (= 2Ω cos θ in terms of the colatitude) and the horizontal Coriolis parameter f = 2Ω cos φ (= 2Ω sin θ) are included. Here 1/2 and θ = ξ/R are the distance and the colatitude of an arbitrary point (x, y, z) to the polar axis (center line of the box where x/x c = y/y c = 1), and R is the outer radius of the sphere. The horizonal Coriolis parameter f is small near the pole and turns out not to be influential, but it is included for completeness of the three-dimensional equations. In the spherical shell to rectangular box mapping, the ratio H/R, where H is the thickness (or height) of the box, is a relevant parameter. For the approximation to be valid, both H/R and θ max (the colatitude angle spanned by the half-width of the box) need to be much less than 1. We use the term 'γ-box' to emphasize the finite thickness of the domain. If the thickness of the box can be considered vanishingly small compared to the radius of the spherical shell, this approximation reduces to the ordinary γ-plane approximation (f = 2Ω − Ω(ξ/R) 2 and f = 0).
In our calculation, the side boundaries of the box are periodic and the top and bottom boundaries are stress-free. These conditions ensure that aside from the Coriolis force, no unspecified momentum exchange occurs between the fluid layer and the planet. Furthermore, the stress-free boundary condition is more favorable for the formation of large-scale vortical structures (Guervilly & Hughes 2017;Guzmán et al. 2020). Constant, uniform energy flux F bot is fed at the bottom, and the temperature at the top is held constant and uniform. All quantities are made dimensionless by setting the thickness of the box, the initial temperature, density, and pressure at the top to 1. As a result, velocity is scaled by the isothermal sound speed (square root of the ratio of pressure to density) at the top. The box contains two layers: a convectively unstable layer in the lower part and a convectively stable layer in the upper part. In the convectively unstable layer, where the mean vertical temperature gradient is almost adiabatic (but slightly superadiabatic), the thermal conductivity κ is set to deliver only 50% of the total energy flux so that convection needs to pick up the rest of the flux. Due to the significant role of convection in energy transfer, this layer is identified as a convection zone (CZ). In the stable layer, the thermal conductivity is raised to a level that all of the energy flux can be fully delivered by a sub-adiabatic temperature gradient. It is to be labeled as a radiation zone (RZ). In addition, to model the fast loss of energy in the optically thin tropospheric layer, a Newtonian cooling term is smoothly introduced to replace the conduction term in the RZ. After thermal relaxation, it effectively creates a constant temperature layer near the top of the box. The change of atmospheric stability (and thermal conductivity) occurs at the interface between the CZ and the RZ. The location of this interface is at 0.95 of the total height of the box for Jupiter-like simulations (Cases A-B, see later discussion), and at 0.75 for Saturn-like simulations (Cases C-E).
The initial structures of the gas layers are polytropic, i.e. ρ ∝ T n , where n is the polytropic index (see Hurlburt et al. (1984) for a detailed description of polytropic structure in numerical modeling).
n takes the value 2.128 (equals to the adiabatic value n ad = 1/(γ ad − 1)) in the CZ. The value 9 is adopted in the RZ where a large value of polytropic index is used to mimic the rapid change of density relative to temperature. Note that as long as this number is large, the exact value is not important. The final equilibrium state of the RZ is mainly determined by the Newton cooling which enforces an isothermal layer in most of the RZ.
For the Jovian cases, the CZ and the RZ contain about 4.7 and 0.5 pressure scale heights, respectively. Although our simulations consider a fairly deep density-stratified region, the stratification is still much less than that of Jupiter's outer convection zone (where magnetic effects on the hydrodynamics can be considered small). The Jovian CZ may contain more than 10 pressure scale heights, which is prohibitive for current computer resources (see later discussion on thermal relaxation). Here we intend to take the compressibility effects into account as much as computer resources allow.
The box has a square base; λ is the aspect ratio (lateral dimension over height). The aspect ratio has an important effect on the formation of multiple vortices. In earlier simulations (Käpylä et al. 2011;Chan & Mayr 2013;Guervilly et al. 2014) where aspect ratios were small, a single vortex often filled the box. To produce multiple vortices, a computational box with large aspect ratio is needed; λ is chosen to be 16 in our polygon simulations. Investigations on rapidly rotating convection reveal that the convective Rossby number Ro c (root-mean-square velocity in the CZ of the non-rotating case, computed separately, divided by the thickness of the CZ and by the Coriolis parameter f ) and Reynolds number Re (root-mean-square velocity of the computed flow multiplied by the thickness of the CZ, and divided by the averaged kinematic viscosity) are two key factors determining the formation of large-scale vortices (Guervilly et al. 2014). Ro c compares the strength of convective driving to that of the Coriolis effect. Numerical simulations show that Ro c needs to be smaller than a critical value ∼ 0.25 for long-lived cyclones to be generated (Chan & Mayr 2013). Another necessary condition is that the Reynolds number Re should be 100 (Guervilly et al. 2014). As a result, simulation of multiple cyclones needs large λ, high Re, and low Ro c . All these conditions require large number of grids and long integration time. In our numerical experiments, we find that polygon patterns could only be formed in high resolution simulations with long periods of time evolution (see Appendix B). Each computed case typically requires millions of CPU core hours (tens of millions of time steps) to complete. In the current stage, it is difficult to explore the parameter space for a full delineation of the pattern formation boundaries. Our main purpose here is to use a couple of cases to demonstrate that large-scale circumpolar vortices can be produced and survive in polygon patterns for a long time in a deep convection zone.
Parameters of our simulation cases are listed in Table 1. For all the cases, we choose Ω = 0.5.
We first found the polygonal clustering of vortices through Case A. The parameters were changed to bracket some of Jupiter's physical parameters with Case B. We can compare our model parameters with the physical parameters of Jupiter by substituting physical values of the scaling quantities  Jones 2014). In our simulations, the internal heat flux has been amplified by a factor of about 10000 and the stratification is reduced to about 5 pressure scale heights (through reduction of the physical value of gravitational acceleration by the factors 12.2 and 6.15 for Cases A and B, respectively) to speed up the relaxation. The fraction of radiative flux in the CZ is also raised to speed up the relaxation process. Effectively, these modifications increase the Rossby number relative to that of Jupiter. In dimensionless parameter space, our cases are near the upper boundary region of the convective Rossby number where the clustering phenomenon can occur.
To reduce computational cost, the flow fields were first generated in a smaller f-plane box (aspect ratio λ = 4) with high grid resolution. At the beginning of calculation, a small random perturbation in the velocity field was introduced to initiate the development of convection. When well-defined small vortices were formed, the flow field of the smaller box was then periodically extended to initialize the flow field for the λ = 16 γ-box. The thermal structure, turbulence, and organized flows evolved in a mutually consistent way.

Vortex Configurations
In the polar γ-box, the cyclonic vortices clustered toward the pole, merged, and grew bigger.
Merging was the principal process that facilitated the growth of vortex sizes. After a long period of thermal and dynamical relaxations, the growth process ceased. A stable polygon pattern was finally formed. Figs. 1A and 1B show the temperature structures near the radiative-convective boundaries.
To illustrate the dynamical processes, time evolution of temperature structures in the small and large boxes for case B are shown in MovieS1 and MovieS2, respectively. The former shows the run in the λ = 4 box (for about 98 rotation periods) where intermediate cyclones first appeared, and the latter shows the extended run in the λ = 16 box (for over 2600 rotation periods). The bright (dark) color represents higher (lower) temperature.
MovieS2 tracks the evolution of the temperature field from the merger generation of size-saturated vortices to the formation of the hexagonal pattern. At the initial stage, small vortices merge to form larger vortices. The merging process essentially stops after the large vortices grow to saturation size.
Finally the beta effect around the pole pushes the cyclones in to form a stable hexagonal pattern.
From the movie, we also observed that a nine-vortex pattern, similar to but not as well organized as the octagonal pattern on Jupiter's north polar region, lasted for a while.
The While the appearance of different polygons in the two cases illustrates that the cluster patterns need not be unique, the general characteristics bear similarity to those observed by the Juno spacecraft in the southern polar region. The octagon pattern of cyclones in the northern polar region  has not yet been replicated in this paper, but it would not be surprising that a suitable change in some parameters can make that realized. We postpone this for future research.

Motions of Vortex Centers
The cyclones are labeled as PC0-PC5 for Case A and HC0-HC6 for Case B. The positions of the centers of the cyclones (locations of local temperature minimum) in the two cases are shown in Fig. 4.
The center of PC0 (the central vortex of Case A) is not exactly at the pole (x c = y c = 8); there is a small offset in the direction away from the largest cyclone PC4 (Fig. 4A). The small offset of the pentagon pattern in Jupiter's south polar region was also observed by Juno ).
Despite the stable structure, small regular movements can be detected. One interesting finding is that PC0 drifts around its time averaged center in an anticyclonic fashion, within a distance of 0.75 • .
In Fig. 5    As the overall structures of the vortices vary little vertically due to dominance of rotation (convective Rossby number is much smaller than one), one may consider the motions in the medium scale (vortex scale) as quasi-two-dimensional, so that some simple analysis can be applied to explore the essential factors of the dynamics.    centers, their averaged effect is to produce a circulation in the clockwise direction (anticyclonic).
The same procedure can be performed for the off-center vortices and the qualitative features are the same. One can make the analysis more quantitative by decomposing the tangential velocity on a circle into modes with different azimuthal wave numbers. On a circle with given radius, we perform a Third, at the radial distance where the circulation velocity peaks (around r = 1), the vortex Rossby numbers V c peak /(2Ωr peak ) of PC0 and HC0 at this height are about 0.4. They are larger than the value 0.24 estimated by Li et al. (2020) according to observation (Grassi et al. 2018).
Given that the vortex flow is almost axisymmetric up to a certain known radial distance, we can proceed to perform a stability analysis within the distance, based on analogy with an idealized axisymmetric flow having the same circulation velocity distribution. An inertial stability (Yanai 1964;Schubert & Hack 1982) criterion in rather general settings (stable and unstable backgrounds) has been given by Kloosterziel et al. (2007). Here we focus on a convectively unstable background for which a necessary condition for inertial stability is (f + 2V c /r) (f + V c /r + dV c /dr) > 0 (see Appendix C), where r is the distance to the center of the vortex, V c is the axisymmetric flow velocity around the center, f is the Coriolis parameter (or planetary vorticity), V c /r is the curvature vorticity (or mean relative angular velocity), and dV c /dr is the shear vorticity. For a system of multiple vortices, this linear stability criterion can be applied to flows around each of the vortices. If any of them violates the criterion, the system cannot be maintained. The combined requirement becomes a necessary condition for stability of the system.
We take Case A as example here. The mean tangential velocity V c versus r is plotted for all the cyclones in Fig. 9A (on the z = 0.5 level). V c is always positive in the near field. At larger distances, V c turns negative as the circle for picking up the axisymmetric part of the circulation velocity expands into neighboring vortices (Fig. 7A). The first term enclosed in parentheses on the left-hand-side of the stability inequality is always positive as the magnitude of V c /r is much less than f (∼ 1 for the current case) for large r (Fig. 9B). The second term enclosed in parentheses, however, contains consequential negative terms created by close packing of the vortices, but inertial stability requires positivity of this factor. The shear vorticity dV c /dr makes the greatest negative contribution to the second term enclosed in parentheses. The radial drop of V c is expedited by the shortened separations among the neighboring vortices. The tighter the packing, the more negative dV c /dr may reach. Despite the positivity of the curvature vorticity, the sum V c /r + dV c /dr is negative beyond r = 1 for PC0 and r = 1.8 for PC1 to PC5 (Fig. 9C). These radial distances are within the ranges where the vortex can safely be regarded as axisymmetric (see Fig. 8). f is the only term to keep the stability factor f + V c /r + dV c /dr everywhere positive (Fig. 9D). Thus it plays a critical role in allowing the cyclones to stay closely packed. A corollary is that cyclones cannot be packed closely for long in low latitudes. The analysis is similar for the hexagon case (see Fig. 10).
Recently, Li et al. (2020) found that in shallow water formulation anticyclonic 'rings' need to be introduced around individual cyclones to shield the cyclones from merging. This conclusion is compatible with the need for a significant Coriolis parameter. The rings cannot be stable without a large enough f; that was implicitly satisfied by the selected vortex profiles. Around the computed vortices of our cases, anticyclonic 'rings' do appear in the averaged sense (see Figs. 7, 9C, and 10C), and they arise naturally as a result of packing the vortices closely.
In the above discussion, the stability criterion is applied to each vortex separately. Global stability analysis of vortex polygons has been made by Reinaud (2019). Based on the quasi-geostrophic (QG) and inviscid assumptions, the author has shown that three-dimensional like-sign uniform PV (potential vorticity) vortices can remain stable in a m-fold axially symmetric configuration (with a moderate strength central vortex) for m < 7. The stable polygon patterns obtained by our simulations are in line with his conclusion. In contrast with the QG theory, our simulation includes viscous, nonadiabatic, turbulent effects, and the vortices are generated by the rotating turbulent convection.
Under the more complicated situation, the maintenance of closely packed vortices is more challenging because of the following two factors. First, the vortices could be destroyed by energy arising from be created and self-organized to a polygonal distribution at the polar region of a fast rotating, turbulent convection zone. It offers a viable mechanism to explain the long-lived closely-packed vortices observed in the polar regions of Jupiter.  (Iess et al. 2019). Therefore, the Ro c of Saturn can be estimated as ∼ 0.2 times that of Jupiter, which means that the relative influence of rotation on
This is an important question relevant to deciphering the vortex generation mechanism. Based on models with thin, stable layer settings, lower Rossby number should create more cyclones. But in deep convection models, the thickness of the CZ is non-negligible and relevant. Besides the convective Rossby number, the depth of the convection zone also plays an important role in vortex formation (correspondingly the aspect ratio in an experimental setup) . If the depth of a γ-box increases, the aspect ratio decreases (the horizontal extent of the box has to stay within a reasonably range of colatitude from the pole). A box with small aspect ratio cannot accommodate many vortices.
Though dependent on the Rossby number, the intrinsic aspect ratio (diameter over depth) of a deep vortex in the polar region, at least for the range of Ro c reachable so far, is on the order of 3 (see for example Fig. 3). As the depth of Saturn's convection zone is estimated to be three times deeper than that of Jupiter, we expect that the aspect ratio of Saturn's polar convection zone is about one third of that of Jupiter. Its polar vortex is thus much larger and the number is limited to 1. This argument ties the number of long-lived vortices directly with the depth of the convection zone relative to the planetary radius (the parameter H/R).
To examine the argument above, we have performed some preliminary numerical experiments for Saturn-like configurations. Table 2 lists the parameters of three additional simulation cases (labeled by C, D, and E). The particular characteristics of these cases are the three-fold increase of H/R, the doubling of θ max , and the successive decrease of Ro c (the minimum value is four times lower than that of Table 1). The number of horizontal grids is smaller, as the aspect ratio of the computational box is smaller (by a factor of 9/16 relative to that of Cases A and B; also see Appendix B). The left-side panels of Fig. 11 show examples of horizontal cuts of the vertical vorticity at z = 0.85 (the stability interface is now at z = 0.75) for the three cases C, D, and E in sequential order; the right-  (2016)); they will be very useful for closer comparison with observational results.

A. DIMENSIONAL QUANTITIES
As the model is based on similarity of dimensionless parameters, the quantities are made dimensionless. We scale all the physical variables by the pressure, density, temperature, and height at the top of the box. In Table 3 p scl , ρ scl , T scl , and L scl represent the scales of pressure, density, tempera-ture, and length, respectively. Consequently, the velocity scale is v scl = (p scl /ρ scl ) 1/2 ; the time scale is t scl = L scl /v scl ; and the flux scale is F scl = ρ scl v 3 scl . For comparison with dimensional quantities of Jupiter and Saturn, the values of the essential scale factors are listed. A dimensional quantity can be obtained by multiplying the dimensionless value with the corresponding scale factor. large-scale cyclone is found in the coarser-grid case. Besides illustrating the importance of resolution, this comparison also shows that the horizontal resolution of 900 2 is adequate to accommodate more than one roaming cyclone in a polar f-plane box with moderate aspect ratio, when there is no polar beta effect to push them to merge at the pole (as in Cases C-E).
where ψ = exp(ωt)Ψ(r, z) is the meridional streamfunction of the perturbed velocity, ω is the time frequency, r is the radius away from the vortex center, z is the height, N 2 is the Brunt Väisälä frequency, Φ = (f + 2V c /r)(f + V c /r + dV c /dr) is the Rayleigh discriminant, and V is the volume of the domain of Ψ. With some rearrangement, the time frequency square can be expressed as where The flow is stable to axisymmetric modes when ω 2 < 0. As I 2 is positive definite, a stable vortex requires I 1 > 0. The sign of I 1 depends on N 2 and Φ. If N 2 > 0, then Φ > 0 is a sufficient but not a necessary condition for inertial stability (under axisymmetric perturbations). If N 2 = 0, then Φ > 0 is a sufficient and necessary condition. If N 2 < 0, then Φ > 0 is a necessary but not sufficient condition for inertial stability. Though our fluid is non-Boussinesq, this criterion provides a helpful analytical framework for theoretical understanding.
D. SUPPLEMENTARY MOVIES MovieS1 (Fig. 13 is the static figure) shows the generation of small vortices in a small f-box. The aspect ratio of the box is 4. The time period of this movie is 1332 (about 98 planetary rotation periods). MovieS2 (Fig. 14 is the static figure) shows the generation and formation of the hexagon pattern by the closely-packed cyclones. The aspect ratio of the box is 16. The time period of this movie is over 33000 (over 2600 planetary rotation periods). MovieS3 (Fig. 15 is the static figure) shows the time evolution of the pentagon pattern. The time period of this movie is 1200 (around 96 planetary rotation periods). In all the movies, the bright (dark) color represents higher (