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

NASA’s Juno spacecraft has discovered closely-packed cyclones on both the North and South polar regions of Jupiter (Bolton et al., 2017; Adriani et al., 2018; Tabataba-Vakili et al., 2020). The Northern Polar Cyclone, about 4000 km in diameter, is surrounded by eight circumpolar cyclones with similar sizes. The Southern Polar Cyclone was surrounded by five circumpolar cyclones, with diameters 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 (Adriani et al., 2018). 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 $\gamma$-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 $\beta$-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 $\gamma$-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 $\gamma$-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.

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

where $\rho$ is the density, $\bm{v}$ is the velocity, $p$ is the pressure, $\bm{\Sigma}$ is the viscous stress tensor, $\bm{g}$ is the gravitational acceleration, $\bm{\Omega}$ is the angular velocity, $E=e+\rho\mathbf{\it{v}}^{\mathbf{2}}/2$ is the total energy density, $e$ is the internal energy, $\bm{F_{\rm d}}$ is the diffusive energy flux, $T$ is the temperature, $c_{p}$ is the heat capacity under constant pressure, $\tau$ 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=\rho 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 $\gamma_{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 $\nu$ according to the Smagorinsky scheme (Smagorinsky, 1963)

where $\Delta x$, $\Delta z$ are local sizes of the horizontal and vertical grids, $\bm{\sigma}$ is the strain rate tensor, the double dot denotes tensor contraction, and the coefficient $c_{\nu}$ 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 ‘$\gamma$-box’ is considered to represent a piece of spherical shell with polar axis through its center, along the vertical direction. The local gravitational acceleration $\bm{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 $\bm{\Omega}$ can be tilted with respect to $\bm{g}$. The inclined angle between $\bm{\Omega}$ and $-\bm{g}$ is the colatitude $\theta$, and the latitude is $\phi=90^{o}-\theta$. The form of the vector equations remain unchanged, but $\bm{\Omega}$ now depends on location. As the box is three-dimensional, the effects of both the vertical Coriolis parameter $f=2\mathrm{\Omega\sin\phi}$ ($=2\mathrm{\Omega\cos\theta}$ in terms of the colatitude) and the horizontal Coriolis parameter $f^{\prime}=2\mathrm{\Omega\cos\phi}$ ($=2\mathrm{\Omega\sin\theta}$) are included. Here $\mathrm{\Omega}=\|\mathbf{\Omega}\|$, $\xi=((x-x_{c})^{2}+(y-y_{c})^{2})^{1/2}$ and $\theta=\xi/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^{\prime}$ 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 $\theta_{max}$ (the colatitude angle spanned by the half-width of the box) need to be much less than 1. We use the term ’$\gamma$-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 $\gamma$-plane approximation ($f=2\Omega-\Omega(\xi/R)^{2}$ and $f^{\prime}=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. Time scale is given by the amount of time for this isothermal sound speed to travel a distance equal to $H$. To obtain time in terms of the rotation period of the planet, one can multiply the dimensionless time with the dimensionless $\Omega/(2\pi)$. Dimensional quantities can be recovered once the values of the scaling quantities at the top level are specified (see Appendix A).

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 $\kappa$ 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. $\rho\propto 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/(\gamma_{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; $\lambda$ 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; $\lambda$ 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 $\sim 0.25$ for long-lived cyclones to be generated (Chan & Mayr, 2013). Another necessary condition is that the Reynolds number $Re$ should be $\gg 100$ (Guervilly et al., 2014). As a result, simulation of multiple cyclones needs large $\lambda$, 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 ${\Omega}=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 into the scaling factors. The surface temperature, density, and pressure of Jupiter are about $166$K, $0.167\rm kg/m^{3}$, and $10^{5}\rm Pa$, respectively. With these as dimensional values of the scaling quantities at the top of the box, the velocity scale is then $\approx 773.8\rm m/s$. The radius of Jupiter is about $70000$ km. The box thickness of case A, being 0.0524$R$, is about 3600km (the thickness of the stable layer is about 180 km), and the width is about 58000km. The scale for time in this case is about $4740\rm s$, with which and the rotation period can be found to be approximately $16.5$ hours. For Case B, the box thickness is about 1800km (the thickness of the stable layer is about 90km), and the width is about 29000km. The planetary rotation period is about 8.3 hours; it is not equal but closer to Jupiter’s rotation period of about 10 hours (Bagenal et al., 2007).

Both cases deal with deep stratification. Simulation with realistic values of internal heat flux and stratification of Jupiter is extremely difficult. Even for the shallower case $H/R=2.6\,\%$, the stratification spans 10 pressure scale heights. The observed internal heat flux of Jupiter is about $5.4-7.5\rm Wm^{-2}$ (Hanel et al., 1981; Li et al., 2018). The thermal relaxation time can be estimated by the Kelvin-Helmholtz timescale $t_{KH}=\int edz/F_{bot}$, which is over $O(10^{5})$ years (and over $O(10^{8})$ rotation periods). So far no numerical simulation of Jupiter’s deep convection zone can get close to such numbers. In practice, larger values of heat flux and smaller stratification are adopted (e.g. 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 $\lambda=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 $\lambda=16$ $\gamma$-box. The thermal structure, turbulence, and organized flows evolved in a mutually consistent way.

Although the focus of this paper is on Jupiter, the question about comparison with Saturn, which also has a deep convection zone, invariably arises. For deep convection models, the primary factor for the appearance of large-scale vortices is the convective Rossby number. The internal heat flux of Saturn is about one third of that of Jupiter. As convective velocity in a non-rotating CZ is approximately proportional to the cubic root of the convective energy flux (Böhm-Vitense, 1958; Chan & Sofia, 1989), the velocity scale for calculating Saturn’s $Ro^{c}$ would be roughly 0.7 times that of Jupiter. While the rotation period of Saturn is almost the same as that of Jupiter, Saturn’s CZ is estimated to be about three times deeper (Iess et al., 2019). Therefore, the $Ro^{c}$ of Saturn can be estimated as $\sim 0.2$ times that of Jupiter, which means that the relative influence of rotation on Saturn is even stronger. Why did observations (Godfrey, 1988; Vasavada et al., 2006; Sánchez-Lavega et al., 2014) find only single long-lived cyclones in Saturn’s polar regions?

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 $\gamma$-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 $\theta_{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-side panels of the same figure show corresponding horizontal velocity streamlines. The streamlines of Case C (Fig. 11B) indicate some similarity with the polar vortex obtained recently by the model of Yadav & Bloxham (2020), which was an anelastic global simulation of deep rotating convection ($H/R=0.1$, 5 density scale heights). They found a single large cyclone at the pole surrounded by a band of three large anticyclones. Our cases show that as $Ro^{c}$ decreases, more anticyclones appear in the surrounding band and polygonal features may arise in flows around the polar vortex, within which the distribution of positive PV broadens. The associated phenomena are intriguing but we postpone an investigation to later study.

In summary, we have shown that numerical computation solving the full set of hydrodynamics equation for rapidly rotating deep convection could produce both Jupiter-like closely-packed circumpolar cyclones and Saturn-like solitary polar cyclone. The convective Rossby number and the lateral-to-height aspect ratio are confirmed to play critical roles in the formation and the selection of the number of these polar cyclones. For each of the cyclones in the Jovian cases, we analysed the axisymmetric component of the tangential velocity field around the vortex and found that a large enough Coriolis parameter is a necessary condition for allowing closely-packed, multiple cyclones to coexist.

Due to the limitation of computer resources, we only give a few examples to illustrate the physical processes of how deep polar vortices on Jupiter and Saturn may be formed and maintained. Quantitative comparisons with observations have shown that these models are generally inadequate. One major reason could be the very large energy fluxes introduced for fast relaxation (so that velocities become higher and changes are faster). Simulations with more physical realism need more powerful codes (e.g. Cai (2016)); they will be very useful for closer comparison with observational results.