Deep convection-driven vortex formation on Jupiter and Saturn

Keeping the general features of the giant planet interiors in mind (Figure 1), we explore two cases which will help us to better understand vortex formation in rotating spherical shells. We first consider the dynamics in a thin rotating spherical shell – a generic representation of the outer convective layers in giant planets that couple only weakly with the interior magnetic field [22, 37] – which spans a region from 0.97$r_{o}$ to $r_{o}$, where $r_{o}$ is the outermost radius of the shell. Such a shell thickness would be equivalent to about 2000 km and 1000 km if scaled to Jupiter and Saturn, respectively. Since contemporary 3D simulations can not attain the highly turbulent geostrophic regime present in the interiors of giant planets, we need to set up the simulation such that it can promote relevant physics. The choice of assuming a relatively thin layer is motivated by the fact that 2D and quasi-2D rotating flows are well known for vortex formation in a broad parameter regime [38]. The thin shell nature of the simulation will allow it to excite vortices but still retain the 3D nature of the convective flows. We ignore the fluid compressibility in this simulation for the sake of simplicity and to reduce the computational costs associated with the model. Furthermore, a four-fold symmetry is imposed in the azimuthal direction to limit the computational requirements. The crucial control parameters governing the employed non-dimensional Boussinesq equations are: the Ekman number defining the ratio of viscous and Coriolis forces is $10^{-5}$, the Rayleigh number defining the convective driving is $3\times 10^{7}$, and the thermal Prandtl number is 0.1. Additional model details can be found in the Materials & Methods section.

This thin layer model self-consistently generates about seven alternating zonal jets in each hemisphere with an eastward jet in the equatorial region (Supplementary Figure S1), as well as a number of compact vortices. Alternating zonal jet streams have already been reported in earlier simulations of spherical shells representing the deep convection in the outer atmospheres of giant planets (e.g. see [20, 21, 22, 23]). Furthermore, anticyclonic vortices, along with zonal flows, were also reported earlier [33] in a similar but thicker shell setup with an overlying stably stratified layer.

The simulation generates about 40 well defined, plus additional smaller and shorter-lived, vortices in each hemisphere (Figure 2). The vortices are formed in the high to mid-latitude regions (Supplementary Figure S2) and the lower latitudes are dominated by stronger zonal jets. This may be expected since the topographic $\beta$-effect (i.e., changing height of an axially aligned fluid column with latitude) is larger at low latitudes which promotes the development of zonal jets [39]. Around mid latitude a clear wavy pattern is present in the vicinity of the peak of an eastward jet (Figure 2) that roughly separates the low-latitude jet regime from the high latitude vortex regime. This feature is reminiscent of the ‘Ribbon wave’ observed on Saturn [40] at similar latitudes. A persistent cyclonic vortex was present throughout the simulation at or close to both rotational poles.

The dynamics present in the simulation is very rich and can be better appreciated by visualizing a time evolution of the flow; the Supplementary Animation 1 (https://youtu.be/bmBdI9sJbqc) shows the time evolution of mid-depth radial vorticity over 80 rotations. Well defined vortices can be easily tracked as they are advected by the background zonal flow. The vortex life time appears to be mainly dependent on the size, with smaller ones lasting for several days while some large ones, e.g. the anticyclone highlighted by the second black arrow from the north pole in Figure 2, persisting for the entire duration of the animation. During this time evolution, this stable anticyclone experienced a slow drift to lower latitudes. The animation shows instances where storms with the same sign of vorticity merge with each other, leaving a single storm, as well as where storms with opposite vorticity sign merge and vanish. In some instances, storms come very close to each other but are advected away quickly by the background zonal flow. Generally, both cyclones and anticyclones form all over the mid-to-high latitudes, however, (anti)cyclones dominate in the region where the local zonal flow shear direction is (anti)cyclonic. This tendency is evident in the 65N to 85N and in the 50N to 65N latitude ranges where a preference for anti-cyclones and cyclones, respectively, is visible. Such a vortex dynamics is reminiscent of the interactions of storms on both Saturn and Jupiter [27].

We contrast the results from this thin shell case with another study by Heimpel et al. (2016) [33] that also reported vortices in a rotating spherical shell with convective forcing. There, all the vortices were anticyclonic due to the underlying formation mechanism: a stably stratified thin layer in the outermost part of the shell makes radially outward convective plumes to diverge, which, later, form shallow anti-cyclonic circulations in the outer layers due to the Coriolis force. This vortex formation mechanism is different from the one in our simulation where the vortices are self-organization phenomenon promoted by the inverse cascade of energy as seen in earlier Cartesian box studies [34, 35, 36]. Storms in our simulation are deep geostrophic structures and change little along the direction of the planetary rotation (Supplementary Figure S3). Furthermore, the shallow stably-stratified-layer mechanism [33] only produces anti-cyclones, whereas, in our simulation both cyclones and anticyclones form depending on the local zonal flow shear.

Comparing this case with other simulations we performed with thicker spherical shells (not reported here) leads us to believe that the thin shell nature of the shell is helping to promote the self-organization of compact vortices. As mentioned earlier, the main idea behind choosing a thin layer is to promote vortex formation in the parameter regime we can access. We know that 2D layers can promote vortices in a broad parameter regime [38] but lack the complex convection dynamics. In a deeper spherical shell simulation setup (akin to Saturn’s deep atmospheric layer) we can not achieve the highly turbulent and geostrophic regime that is likely needed to excite self-organization of vortices. Therefore, our thin shell setup, retaining properties of both 2D and 3D setups, is a good compromise to simulate relevant physics.

Next, we investigate what happens when a hydrodynamic layer directly interacts with another layer where the electrical conductivity is large enough to allow generation of magnetic fields (Figure 1). We consider a rotating spherical shell that spans a region from 0.1$r_{o}$ to $r_{o}$, within which the fluid density decreases by a factor of about 150 along the radius in the shell. Similar to Dietrich & Jones (2018) [41], the electrical conductivity of the fluid changes with radius as a hyperbolic tangent function: it is constant until $\approx 0.84r_{o}$, then sharply decreases by seven orders of magnitude within $\approx 0.07r_{o}$, and stays constant afterwards (see ‘Profile 1’ in Supplementary Figure S4). Such an electrical conductivity profile defines a deep dynamo layer and an overlaying atmospheric layer where Lorentz forces will be negligible. Therefore, we label the radial level at $0.85r_{o}$ as the ‘transition’ radius separating the dynamo and the atmospheric region. Note that our chosen electrical conductivity profile differs from several earlier studies: the hyperbolic tangent function introduces sharper conductivity gradients in the transition radius as compared to a combination of constant for some radius and exponentially decreasing afterwards [42, 43, 44, 45]. The electrical conductivity in Jupiter starts to decay with a super-exponential rate after about 0.9$r_{o}$ [8] which leads to decoupling of the atmosphere from the dynamo in a relatively small radius range. In order to capture such a decoupling, it is crucial to include a strong electrical conductivity decay offered by a hyperbolic tangent function with a large step change. We then solve the convective magnetohydrodynamics in this setup. The crucial parameters governing the anelastic system of equations are: the Ekman number is $10^{-6}$, the Rayleigh number is $4\times 10^{9}$, the thermal Prandtl number is 0.1, and the magnetic Prandtl number is 1.5. We refer the reader to the Materials & Methods section for more details about the model.

In this two-layer setup, the dynamics and the emergent features are very different from the earlier thin shell case. The simulation generates a largely dipole dominant magnetic field in the interior dynamo region (Supplementary Figure S5). Unlike the previous simulation exhibiting multiple alternating zonal jets, here we only get one equatorial eastward jet extending from about 30 degree north to 30 degree south (Supplementary Figure S6a). Investigating the zonal flow profile throughout the interior helps to explain the surface behavior. The strong Lorentz forces associated with the deeper dynamo layer enforce a nearly rigid body rotation in the fluid below $\approx 0.85r_{o}$. Above this radial level, a hypothetical magnetic tangent cylinder (MTC; Figure 1), that is tangent to the transition radius and parallel to the rotation axis, demarcates two regions [41]. For cylindrical radii larger than MTC, i.e., outside the MTC, the flow sharply decouples from the magnetic field due to much lower conductivity and creates a strong eastward zonal jet. For smaller cylindrical radii, i.e., within the MTC, the zonal flow is weak despite the flow-field decoupling. This property can be understood as following: an axially aligned fluid column in region outside the MTC experiences very little friction from the free slip mechanical boundaries imposed on its two ends, allowing it to develop strong zonal flows; inside the MTC, however, although the fluid column has a free-slip mechanical boundary on the upper end, it has a dynamo-imposed nearly rigid rotation on the lower end that mimics strong boundary friction, prohibiting the development of zonal flows [21, 43]. We note that Jupiter and Saturn have alternating zonal flows at all latitudes. This is not a contradiction. It is generally believed that, unlike in the contemporary high-viscosity simulations, the zonal jets on Jupiter and Saturn experience much smaller drag due to the extremely small fluid viscosity [8]. This allows the jets to readily develop at all latitudes in the molecular envelope of these planets despite having a dynamo-imposed nearly rigid body rotation in the deep [46]. Recent laboratory experiments on jet formation with low viscosity fluids lend support to these ideas [47].

Analyzing the flow pattern on the outer boundary of the simulation reveals interesting dynamics. The most dominant feature on the surface is the eastward jet in the equatorial region. In mid-to-high latitudes energetic vortical structures are present, which stand out in the otherwise uniformly distributed background of convective plumes (Figure 3). The energetic vortices in both the north and the south latitudes are anticyclonic in nature. A closer inspection of the vortices reveals that they have a quiet core with small velocities, surrounded by a large fast-rotating rim (Figure 3c, top right). Some of the large vortices also show a more complex behavior where a clear center and a rim are less well defined, but overall the structure constitutes an anticyclonic vortex (Figure 3c, top left). These vortices form on scales ranging from 5 to 20 degrees in latitudinal span (equivalent to about 6000 to 24000 km in diameter if present on Jupiter’s surface). The life time also varies with small storms persisting for only one planetary rotation while some larger ones lasting for 20 or so (see Supplementary Animation 2; https://youtu.be/Azlfva56jb8). Some of the long lasting storms at low latitudes also experience a westward drift of about 0.5 deg/day.

Investigating the structure of a vortex in deeper layers reveals the underlying formation mechanism: as a convective plume rises from deeper depths it impinges on the outer boundary and diverges, the background Coriolis force then transforms the horizontal divergence to an anticyclonic circulation (Figure 3d). Therefore, although the anticyclones are relatively shallow and have a pancake shape, they have a very deep origin.

Presence of such large anticyclonic vortices has not been reported in conventional spherical shell simulations of the outer convective layers of giant planets [19, 20, 21, 22, 23], therefore, the deeper dynamo layer is the key for triggering them. The electrical conductivity jump at the transition radius essentially couples two dynamically distinct layers. The convection in the deep dynamo continuously launches plumes into the overlying layer. Some of these plumes manage to either cross the entirely of the atmospheric layer and impinge on the outer layer or interact with the upper layer to create plumes that go on to hit the outer boundary and create anticyclones. Since the plumes are coming from deeper denser regions, when they deposit their energy in the outer less-dense layers, they can create substantially stronger flows than the surrounding local flows. The importance of the dynamo layer is further underscored by the fact that some of the low-latitude long lasting storms experience a slow westward drift which is driven by the weak westward zonal flow produced in the dynamo region that permeates the atmospheric layer (Supplementary Figure S6b). Such westward drifts are commonly found in dipole-dominant dynamos [48, 49]. Furthermore, the sharp conductivity gradient in our model is also an important ingredient to excite the anticyclones, another simulation with a more conventional profile of constant conductivity till $0.85r_{o}$ and an exponential drop afterwards (‘Profile 2’ in Supplementary Figure S4) did not produce similar anticyclones. This is likely due to the presence of a gradual transition of the fluid dynamics from the dynamo layer to the hydrodynamic layer which makes it one system rather two layers interacting with each other. Since the anticyclones are formed in the outer low-conductivity layer, they do not have a corresponding magnetic signal associated with them (Supplementary Figure S7). However, the simulation does have a tendency to generate more anticyclones in the regions containing higher magnetic energy. This hints that magnetic Lorentz forces are likely playing a role in promoting the formation of the anticylones.

The first simulation discussed in the paper consists of a thin spherical shell bounded by a lower boundary at $0.97r_{o}$ and an outer boundary at $r_{o}$. The spherical shell rotates about a fixed axis, whose direction vector is represented by $\hat{z}$, with an angular velocity $\Omega$. The fluid dynamics in this shell is governed by thermal convection set up by a temperature contract of $\Delta T$ between the inner and the outer boundary. The fluid is assumed to be incompressible, i.e., $\nabla\cdot{\vec{u}}=0$, where $\vec{u}$ is the velocity. We use the classical Boussinesq approximation to model the fluid dynamics. Furthermore, we work with non-dimensional system of equations where the distances are scaled by the shell thickness $d=r_{o}-r_{i}$, time is scaled by $\Omega^{-1}$, and temperature is scaled by $\Delta T$. With these assumptions, the equations governing the velocity and temperature $T$ are:

	$$\displaystyle\frac{\partial\vec{u}}{\partial t}+\vec{u}\cdot\vec{\nabla}\vec{u}+2\hat{z}\times\vec{u}+\vec{\nabla}P=\frac{Ra\,E^{2}}{Pr}g(r)\,T\,\hat{r}+E\nabla^{2}\vec{u},$$		(1)
	$$\displaystyle\frac{\partial T}{\partial t}+{\vec{u}\cdot\vec{\nabla}}T=\frac{E}{Pr}\nabla^{2}T,$$		(2)

where $P$ is pressure. The non-dimensional gravity $g(r)$ varies as $(r_{o}/r)^{2}$ assuming that most of the planetary mass is concentrated below the inner boundary at $0.97r_{o}$. The relevant non-dimensional parameters are: Rayleigh number $Ra=\alpha g_{o}\Delta Td^{3}/(\kappa\nu)$, Ekman number $E=\nu/(\Omega D^{2})$, and Prandtl number $Pr=\nu/\kappa$, where $\alpha$ is the thermal expansion coefficient, $g_{o}$ is gravity at outer boundary, $\nu$ is viscosity, and $\kappa$ is thermal diffusivity. These control parameters are set to: $E=10^{-5}$, $Ra=3\times 10^{7}$, and $Pr=0.1$. The velocity matches a stress-free condition on both boundaries. The temperature is assumed to be constant on each boundary.

The above system of equations is solved using the open source MagIC code (https://magic-sph.github.io/; [56]) which employs the pseudo-spectral approach. The latitudinal and longitudinal quantities are expanded using the Legendre polynomials and the radial ones are expanded using the Chebyshev polynomials. The code also uses the toroidal-poloidal decomposition to maintain strict divergence-less nature of the relevant quantities. The code uses the fast spherical harmonic library SHTns (https://nschaeff.bitbucket.io/shtns/; [57]). The equations are time-advanced using an explicit second-order Adams-Bashforth scheme for the Coriolis and the non-linear terms, and an implicit Crank-Nicolson scheme for other terms [58].

We simulate the system for about 1500 rotations which was enough to show a statistically stationary behavior. To limit the computational requirements, we used four fold symmetry in the azimuthal direction [59, 22], which effectively makes it a quarter-spherical-wedge simulation with periodic boundary conditions on the meridional planes of the wedge. The first 930 rotations were carried out on a simulation grid of [1056, 2112, 160] where the numbers represent resolution in azimuthal, latitudinal and the radial direction, respectively. The corresponding maximum spherical harmonic degree is 1408. The remaining time period was simulated on a larger grid of [1152, 2304, 200]. Apart from sharpening the fluid dynamic features, the higher resolution did not change the simulation behavior. Despite the large grid we used, a hyperdiffusion scheme, which suppresses energies at small length scales, had to be applied to run the simulation. Following several earlier studies [23, 22, 33], the viscosity in our setup becomes a function of spherical harmonic degree after a certain cutoff. MagIC multiplies the following function to the main viscous diffusion operator:

$$\displaystyle d(\ell)=1+D\left[\frac{\ell+1-\ell_{hd}}{\ell_{max}+1-\ell_{hd}}\right]^{\beta}$$

(3)

where $D$ defines the amplitude of the function, $\ell_{max}$ is the maximum spherical harmonic degree utilized in the simulation, $\ell_{hd}$ is the degree after which the hyerdiffusion starts, and $\beta$ defines the rise of the function for degrees higher than $\ell_{hd}$. For our simulation, we use $D=10$, $\beta=6$, $\ell_{hd}=450$.

In this simulation we use the anelastic system of equations [60, 61] to model compressible, subsonic flows present in the interior of giant planets. In this approximation, the thermodynamic quantities are assumed to be a combination of a static background (tilde) and small fluctuations (prime), i.e., $x=\tilde{x}+x^{\prime}$. The equations used in our simulation are in the entropy variable form with entropy contrast between the inner and the outer boundary driving the convection [37]. The magnetic field is scaled by $\sqrt{\tilde{\rho_{o}}\mu_{o}\lambda_{i}\Omega}$, where $\tilde{\rho_{o}}$ is density on the outer boundary, $\mu_{o}$ is magnetic permeability, and $\lambda_{i}$ is magnetic diffusivity on the inner boundary. The scales for distance and time are the same as above. With these assumptions and rescaling, the evolution for velocity is governed by:

	$$\displaystyle\nabla\cdot(\tilde{\rho}{\vec{u}})=0$$		(4)
	$$\displaystyle\left(\dfrac{\partial\vec{u}}{\partial t}+\vec{u}\cdot\vec{\nabla}\vec{u}\right)=-\vec{\nabla}{\frac{p^{\prime}}{\tilde{\rho}}}-\dfrac{2}{E}\hat{z}\times\vec{u}+\dfrac{Ra}{Pr}\tilde{g}\,s^{\prime}\,\hat{r}$$
	$$\displaystyle+\dfrac{1}{Pm_{i}\,E\,\tilde{\rho}}\left(\vec{\nabla}\times\vec{B}\right)\times\vec{B}+\dfrac{1}{\tilde{\rho}}\vec{\nabla}\cdot\mathsf{S},$$		(5)

where $p^{\prime}$ is pressure perturbation, $\tilde{g}(r)$ is gravity (assumed to vary linearly), $s^{\prime}$ is entropy perturbations, $\vec{B}$ is magnetic field, and $Pm_{i}$ is $\nu/\lambda_{i}$. The Rayleigh number $Ra$ in this formulation is given by $\alpha_{o}g_{o}T_{o}d^{3}\Delta s(c_{p}\kappa\nu)^{-1}$. The traceless rate-of-strain tensor $\mathsf{S}$ is defined by

$$\displaystyle S_{ij}=2\tilde{\rho}\left(e_{ij}-\frac{1}{3}\delta_{ij}\vec{\nabla}\cdot\vec{u}\right)\text{with}\,\,e_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right),$$

where $\delta_{ij}$ is the identity matrix. The energy conservation gives time evolution equation for entropy perturbation $s^{\prime}$:

	$$\displaystyle\tilde{\rho}\tilde{T}\left(\dfrac{\partial s^{\prime}}{\partial t}+\vec{u}\cdot\vec{\nabla}s^{\prime}\right)=\dfrac{1}{Pr}\vec{\nabla}\cdot\left(\tilde{\rho}\tilde{T}\vec{\nabla}s^{\prime}\right)+\dfrac{Pr\,Di}{Ra}\Phi_{\nu}+$$
	$$\displaystyle\dfrac{Pr\,Di\,\lambda_{norm}}{Pm_{i}^{2}\,E\,Ra}\left(\vec{\nabla}\times\vec{B}\right)^{2},$$		(6)

where $\tilde{T}$ is background temperature and the viscous heating contribution is given by

$$\displaystyle\Phi_{\nu}=2\tilde{\rho}\left[e_{ij}e_{ji}-\frac{1}{3}(\vec{\nabla}\cdot\vec{u})^{2}\right].$$

(7)

The dissipation number $Di$ is ${\alpha_{o}{g_{o}}d}/{c_{p}}$, where $\alpha_{o},g_{o}$ are thermal expansivity and gravity on outer boundary and $c_{p}$ is specific heat at constant pressure. The evolution of the magnetic field is governed by

$$\displaystyle\dfrac{\partial\vec{B}}{\partial t}=\vec{\nabla}\times\left(\vec{u}\times\vec{B}\right)-\dfrac{1}{Pm_{i}}\vec{\nabla}\times\left(\lambda_{norm}\vec{\nabla}\times\vec{B}\right)$$

(8)

where $\lambda_{norm}$ is local magnetic diffusivity normalized by its value at the inner boundary $r_{i}$. The magnetic field also follows the divergence-less condition $\vec{\nabla}\cdot\vec{B}=0$. We use the anelastic version [62] of MagIC code to simulate this system of equations.

The inner boundary is at 0.1$r_{o}$ and the outer is at $r_{o}$. The control parameters for this simulation are: $E=10^{-6}$, $Ra=4\times 10^{9}$, $Pr=0.1$, $Pm_{i}=1.5$. The density of the fluid changes by five density scale heights (a factor of about 150) along the radius.

The electrical conductivity follows a step change using a hyperbolic tangent function [41]: it remains constant until about 0.85$r_{o}$, decreases by factors of $10^{7}$ in a small radius, and stays constant afterwards (See ‘Profile 1’ is Supplementary Figure S4). The velocity matches a stress-free condition, the entropy is held constant, and the magnetic field matches a potential field on the shell boundaries. Instead of starting the simulation from small flow and magnetic field perturbations, we used a saturated dynamo simulation (without any conductivity change in its interior) from our earlier study [63] and lowered its Ekman number from $10^{-5}$ to $10^{-6}$ in incremental stages. Once the $10^{-6}$ simulation showed statistical saturation in its kinetic and magnetic energy time series, we then introduced a conductivity jump in stages, reaching a decrease of $10^{7}$ in the outer layer. The simulation was run for about 2000 rotations in the final form. About 50 rotations were simulated on a grid of size [2112, 1054, 400] and the rest were on a grid with [1152, 576, 400]. Here, too, increasing the resolution did not change the simulation nature substantially. This simulation also required hyperdiffusion with parameters $D=20$, $\beta=3$, $\ell_{hd}=450$ (see earlier section) for the higher resolution grid and $D=5$, $\beta=2$, $\ell_{hd}=120$ in the lower resolution grid.