Deep convection–driven vortex formation on Jupiter and Saturn

See allHide authors and affiliations

Science Advances  13 Nov 2020:
Vol. 6, no. 46, eabb9298
DOI: 10.1126/sciadv.abb9298


The surfaces of Jupiter and Saturn have magnificent vortical storms that help shape the dynamic nature of their atmospheres. Land- and space-based observational campaigns have established several properties of these vortices, with some being similar between the two planets, while others are different. Shallow-water hydrodynamics, where the vortices are treated as shallow weather-layer phenomenon, is commonly evoked for explaining their formation and properties. Here, we report novel formation mechanisms for vortices where the primary driving mechanism is the deep planetary convection occurring in these planets. Using three-dimensional simulations of turbulent convection in rotating spherical shells, we propose two ideas: (i) Rotating turbulent convection generates deep axially aligned cyclones and anticyclones; (ii) a deep planetary dynamo acts to promote additional anticyclones, some as large as Jupiter’s Great Red Spot, in an overlying atmospheric layer. We use these ideas to interpret several observational properties of vortices on Jupiter and Saturn.


By modeling the dominant dynamical features present on the surface of Jupiter and Saturn, namely, zonal jet streams and storms/vortices, we can learn about what drives them and their connection to the planetary deep interior. Apart from the well-known and persistent hexagonal storm on Saturn and the Great Red Spot (GRS) on Jupiter, there are numerous other compact vortical storms with various sizes and lifetimes found at various latitudes on these planets. By analyzing the images taken by the Cassini spacecraft through its lifetime, Trammel et al. (1, 2) provide a detailed outlook on the properties of compact vortices on Jupiter and Saturn. In our study, we focus only on the fluid dynamical properties of the vortices where the following general features stand out. First, with the exception of the equatorial regions containing a broad eastward zonal jet, vortices are generally found at all latitudes but tend to occur preferentially in regions of westward zonal flow (1, 3). There is a stark contrast in the number of vortices between the two planets: about 200 with 1000 km or larger diameter on Jupiter, while only 10 to 50 on Saturn (1). Note that while both planets also have smaller vortices, the disparity in the numbers still persists (3, 4). Vortices are also present at or very close to the rotational poles: Saturn has a cyclonic polar vortex at each pole (5), while Jupiter has a cluster of cyclones at each pole (6).

The rotation sense of vortices is also an important property but much harder to infer through direct observations. Li et al. (4) analyzed 500 vortices with a diameter of at least 700 km on Jupiter using Cassini images. They inferred the vorticity of the largest 100 vortices, all of them being anticyclonic. An earlier study (7) using Voyager images reached a similar conclusion. Because of the more hazy atmosphere of Saturn, the contrast of the different fluid dynamic features is much lower compared with Jupiter. Paired with the lack of continuous images of Saturn’s surface, determining the vorticity is not possible for most of the storms. Instead, the vorticity of a compact vortex is assumed to be determined by the vorticity of the local zonal wind shear (13). On the basis of this assumption, Saturn’s atmosphere does not appear to have a strong bias toward anticyclones. We must, however, keep in mind that the number of vortices on Saturn is small and undergoes large variations over time. For example, during the 7-year period from 2008 to 2015, the number of vortices of at least 1000-km diameter changed from about 5 to 20 in the northern hemisphere of Saturn (2). Therefore, a robust statistical trend cannot be inferred for Saturnian vortices.

Much like in the study of zonal jets on gas giant planets [see (810) for shallow jets and (1115) for deep jets], vortices can also be understood as either existing only in the outermost weather layers or extending much deeper into the planetary interior. In the case of zonal jets, exciting progress was recently made by modeling the gravity data from the Juno mission and the Cassini Grand Finale. It indicates that the zonal jet streams on both planets are likely thousands of kilometers deep, about 9000 km in Saturn (16) and about 3000 km in Jupiter (17). These results lend support to the deep zonal jet scenario where planetary convection is the driving force. Given these results, it is certainly worth investigating whether some of the vortices visible on the surface extend deep into the interior.

So far, most vortex formation models consider only the thin outer weather layer, with several using the shallow-water equations with a prescribed forcing in the form of moist convection or solar heating (18, 19). These studies have successfully modeled several properties of the vortices, for example, the bias toward anticyclones (20, 21) and polar vortex formation (22, 23). However, it is also well known that coherent vortices can spontaneously form in three-dimensional (3D) fluid turbulence if rotational effects are substantial [e.g., see (24)]. Could deep planetary convection also generate coherent vortices along with zonal jets? There have been several attempts to study this. Coherent, primarily anticyclonic, vortices have previously been shown to coexist with zonal jets using an elastic, spherical deep convection model with a stably stratified layer (25). Others have investigated vortex formation in Cartesian geometry (2628), where the curvature effects due to the spherical shell are not present. Here, we report several physically motivated simulations of planetary convection in deep spherical shells to better understand how vortices form and behave under such conditions.


Thin shell case

Keeping the general features of the giant planet interiors in mind (Fig. 1), we explore two cases that 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 (14, 29)—which spans a region from 0.97ro to ro, where ro is the outermost radius of the shell. Such a shell thickness would be equivalent to about 2000 and 1000 km if scaled to Jupiter and Saturn, respectively. Since contemporary 3D simulations cannot 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 (30). 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 fourfold symmetry is imposed in the azimuthal direction to limit the computational requirements. The crucial control parameters governing the used nondimensional Boussinesq equations are as follows: the Ekman number defining the ratio of viscous and Coriolis forces is 10−5, the Rayleigh number defining the convective driving is 3 × 107, and the thermal Prandtl number is 0.1. Additional model details can be found in Materials and Methods.

Fig. 1 Gas giant interiors: A schematic representation of the interior of Jupiter and Saturn showing a meridional quadrant in the northern hemisphere.

The partial circles with arrows represent overturning convection. The overturns with thick lines highlight energetic convection plumes that impinge on the adjacent layer. Jupiter may be crudely approximated by a two-layer structure with the deep interior consisting of metallic hydrogen with larger conductivity, and a lower conductivity molecular hydrogen outer layer (38). Saturn, on the other hand, likely has a three-layer structure: an innermost metallic hydrogen layer sustaining the dynamo, a molecular hydrogen-rich low-conductivity outermost layer, and a stably stratified layer due to helium rainout (58) between the other two layers. Although the exact location and nature of the intermediate helium-enriched layer are uncertain, its existence is highly favored (58, 59) due to the very high axisymmetry of Saturn’s magnetic field (60). Jupiter may also have a similar stable layer, but it is likely much thinner (61); furthermore, Jupiter’s magnetic field is not highly axisymmetric (62, 63).

This thin layer model self-consistently generates about seven alternating zonal jets in each hemisphere with an eastward jet in the equatorial region (fig. 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 (1215)]. Furthermore, anticyclonic vortices, along with zonal flows, were also reported earlier (25) 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 (Fig. 2). The vortices are formed in the high to mid-latitude regions (fig. S2), and the lower latitudes are dominated by stronger zonal jets. This may be expected since the topographic β-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 (31). Around mid-latitude, a clear wavy pattern is present in the vicinity of the peak of an eastward jet (Fig. 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 (32) at similar latitudes. A persistent cyclonic vortex was present throughout the simulation at or close to both rotational poles.

Fig. 2 Cyclones and anticyclones: The radial component of the vorticity, defined by (∇ × ur) where u is velocity, plotted at mid depth of the simulation.

The velocity is in the units of Rossby number uro, where Ω is the planetary rotation frequency, and ro is the shell radius. The (−ve) +ve values are represented by the (blue) red color, and the plotted range is ±0.006. The (yellow) green arrows highlight the peaks of (eastward) westward zonal jets. Some of the easily recognizable cyclonic (red) and anticyclonic (blue) vortices are pointed out with arrows, as well as a large-scale modulation labeled as “Ribbon wave feature.” Note that when a quarter wedge is repeated to form a full spherical surface, as is done here, a vortex close to the poles would appear as a cluster of four vortices. This is most likely an artifact since in our other simulations where we considered full azimuthal extent, only one cyclonic vortex was present at or near the poles (64).

The dynamics present in the simulation is very rich and can be better appreciated by visualizing a time evolution of the flow; movie S1 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 lifetime 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 Fig. 2, persisting for the entire duration of movie S1. During this time evolution, this stable anticyclone experienced a slow drift to lower latitudes. Movie S1 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. In general, 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 65°N to 85°N and in the 50°N to 65°N latitude ranges, where a preference for anticyclones and cyclones, respectively, is visible. Such a vortex dynamics is reminiscent of the interactions of storms on both Saturn and Jupiter (19).

We contrast the results from this thin shell case with another study by Heimpel et al. (25), which 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 anticyclonic 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 (2628). Storms in our simulation are deep geostrophic structures and change little along the direction of the planetary rotation (fig. S3). Furthermore, the shallow stably stratified layer mechanism (25) only produces anticyclones, 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 (30) but lack the complex convection dynamics. In a deeper spherical shell simulation setup (akin to Saturn’s deep atmospheric layer), we cannot 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.

Thick shell case

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 (Fig. 1). We consider a rotating spherical shell that spans a region from 0.1ro to ro, within which the fluid density decreases by a factor of about 150 along the radius in the shell. Similar to Dietrich and Jones (33), the electrical conductivity of the fluid changes with radius as a hyperbolic tangent function: It is constant until ≈0.84ro, and then sharply decreases by seven orders of magnitude within ≈0.07ro, and stays constant afterward (see “Profile 1” in fig. 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.85ro 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 with a combination of constant for some radius and exponentially decreasing afterward (3437). The electrical conductivity in Jupiter starts to decay with a superexponential rate after about 0.9ro (38), which leads to decoupling of the atmosphere from the dynamo in a relatively small radius range. 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 as follows: the Ekman number is 10−6, the Rayleigh number is 4 × 109, the thermal Prandtl number is 0.1, and the magnetic Prandtl number is 1.5. We refer the reader to Materials and Methods 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 (fig. S5). Unlike the previous simulation exhibiting multiple alternating zonal jets, here we only get one equatorial eastward jet extending from about 30° north to 30° south (fig. S6A). Investigating the zonal flow profile throughout the interior helps explain the surface behavior. The strong Lorentz forces associated with the deeper dynamo layer enforce a nearly rigid body rotation in the fluid below ≈0.85ro. Above this radial level, a hypothetical magnetic tangent cylinder (MTC; Fig. 1) that is tangent to the transition radius and parallel to the rotation axis demarcates two regions (33). For cylindrical radii larger than MTC, i.e., outside the MTC, the flow sharply decouples from the magnetic field because of 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 follows: An axially aligned fluid column in the 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 (13, 35). 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 (38). 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 (39). Recent laboratory experiments on jet formation with low-viscosity fluids lend support to these ideas (40).

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 (Fig. 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 (Fig. 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 (Fig. 3C, top left). These vortices form on scales ranging from 5° to 20° in latitudinal span (equivalent to about 6000 to 24,000 km in diameter if present on Jupiter’s surface). The lifetime also varies with small storms persisting for only one planetary rotation, while some larger ones lasting for 20 or so (see movie S2). Some of the long-lasting storms at low latitudes also experience a westward drift of about 0.5°/day.

Fig. 3 Deep rooted anticyclones: The velocity magnitude as color maps, flow vectors, and streamlines at different locations.

Lighter colors represent larger velocities with maximum being 0.05 in the units of Rossby number. Panels (A) and (B) show the northern and southern views of the simulation. (C) Velocity magnitude and flow vectors in a subset of the region highlighted with a box in (A). (D) Flow streamlines and flow vectors—white ones for the convective plume triggering the storm and red ones for the rim of the storm—for an anticyclone highlighted with a box in (B). The dark gray surface in (D) highlights the transition radius where the electrical conductivity starts to drop. Furthermore, the translucent cylinder in (D) is aligned with the rotation axis.

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, and the background Coriolis force then transforms the horizontal divergence to an anticyclonic circulation (Fig. 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 (1115); 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 entirety 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 (fig. S6B). These westward drifts are commonly found in dipole-dominant dynamos (41, 42). 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 until 0.85ro and an exponential drop afterward (“Profile 2” in fig. 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 than 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 (fig. 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.


We have presented two simulations of highly turbulent convection in rotating spherical shells. The first simulation demonstrates the spontaneous formation of zonal jets, cyclones, and anticyclones in a thin spherical shell. The second one shows a novel vortex formation mechanism where plumes originating in a deep dynamo layer exclusively excite anticyclones in an overlying low-conductivity atmospheric layer.

The results from these simulations can help us better understand the atmospheric dynamics visible on the surface of Jupiter and Saturn. Although vortex formation (through an inverse cascade of energy) is well known in rotating Cartesian boxes with either forced or convective turbulence (2628), our simulation is among the first to produce numerous simultaneous cyclones and anticyclones, as well as multiple zonal jet streams, in a global spherical geometry with convective driving. In recent work, Heimpel et al. (25) report anticyclones in a fluid shell with a stably stratified layer above a convective layer, Garcia et al. (43) show formation of shell-depth spanning polar vortices in Boussinesq convection, as also reported by Yadav and Bloxham (64) in anelastic convection. Compared with local Cartesian box simulations, our thin shell simulation can be considered a step further toward the goal of modeling the global atmospheric dynamics of gas giant planets. The global geometry allows us to simultaneously capture the vortex-free near-equatorial regions, alternating zonal jets, formation of storms with vorticity governed by local shear flows, and nonlinear interactions between storms (including mergers). The simulation shows that rotating convection can generate deep axially aligned storms that mimic many properties of the storms on Jupiter and Saturn.

The second case shows how a dynamo layer might interact with an overlaying atmospheric layer. The dynamo region is conventionally viewed as an inert layer proving a drag force for the atmospheric zonal jets (1215). Our simulation, instead, shows that it could play an important role dynamically and may provide a source for seeding and sustaining large anticyclones in the atmospheric layer. Given that, there is a possible connection between our results and Jupiter’s GRS. Jupiter’s deep dynamo is expected to be in the magnetostrophic regime where small-scale and planetary-scale giant convective cells are expected (4447). If Jupiter’s dynamo layer does seed anticylones in the atmospheric layer, as suggested by our simulation, then a GRS-type anticyclone with a quiet center can be readily formed. Our simulation forms several GRS-scale anticylones. Furthermore, long-lived vortices in this simulation also showed a GRS-like westward drift: about 0.5°/day in the simulation, while about 0.3°/day for GRS in recent times (48). In the simulation, the deep dynamo sets up a slow westward drift in the low to mid-latitudes. The convection plumes in the dynamo region accordingly drift westward, dragging the correspondingly large vortex in the atmospheric layer. The stability of the GRS, however, could not be captured in the simulation. The most stable storms in our simulation existed for about 20 days, while the GRS has a lifetime of hundreds of years (48). Possibly, Jupiter is able to sustain much longer-lived convective plumes, spawning stable GRS-like vortices, than our simulation. A broader systematic parameter study of this simulation will be useful to investigate the possibility of longer-lived large anticyclones.

On the basis of these two case studies, we can begin piecing together a global picture of the deep convection–driven fluid dynamics occurring in Jupiter and Saturn. An atmospheric layer with dynamics governed by rapidly rotating deep convection will tend to promote either cyclonic or anticyclonic vortices/storms depending on the direction of the vorticity set by the local zonal wind shear. This phenomenon is likely active in the atmospheric layers of both Saturn and Jupiter. On top of that, if the interior structure of a planet allows the dynamo layer to directly interact with the outer hydrodynamic layer, then energetic plumes from the deep dynamo layer will exclusively excite anticyclonic vortices in the atmospheric layer, tilting the balance toward more anticyclones than cyclones. This scenario is likely active in Jupiter, explaining why it has many more anticyclones than cyclones (4, 19). On Saturn, however, the dynamo layer is largely decoupled from the outer hydrodynamic layer due to the likely presence of a thick stably stratified layer. Even if its dynamo layer generates energetic plumes, they get trapped by the overlying stably stratified layer, explaining the absence of a clear and consistent preference for anticyclones on Saturn (2, 3, 19).

The vortex formation mechanisms discussed above, however, do not exclude the presence of other mechanisms that have been proposed. For example, moist convection–driven shallow storms might be forming in the outermost layers of both Saturn and Jupiter. Furthermore, the idea of anticyclone formation through convective pummeling of a thin stably stratified layer from below (25) might also be active on both planets. The atmospheres of Jupiter and Saturn exhibit extremely rich dynamics, suggesting the presence and complex interaction of various storm formation mechanisms. Here, we offer two possibilities that appear promising. The Juno spacecraft currently in orbit around Jupiter will likely shed more light on this issue. The on-board microwave radiometer (49) observes thermal radiation from levels as deep as 1000 bars. If large atmospheric vortices on Jupiter maintain some correlation with the temperature at depth, then we might be able to ascertain if some of the storms on the surface go deeper into the atmosphere.


Thin shell case

The first simulation discussed here consists of a thin spherical shell bounded by a lower boundary at 0.97ro and an outer boundary at ro. The spherical shell rotates about a fixed axis, whose direction vector is represented by ẑ, with an angular velocity Ω. The fluid dynamics in this shell is governed by thermal convection set up by a temperature contract of ΔT between the inner and the outer boundary. The fluid is assumed to be incompressible, i.e., ·u=0, where u is the velocity. We use the classical Boussinesq approximation to model the fluid dynamics. Furthermore, we work with nondimensional system of equations where the distances are scaled by the shell thickness d = rori, time is scaled by Ω−1, and temperature is scaled by ΔT. With these assumptions, the equations governing the velocity and temperature T areut+u·u+2ẑ×u+P=Ra E2Prg(r) T r̂+E2u(1)Tt+u·T=EPr2T(2)where P is pressure. The nondimensional gravity g(r) varies as (ro/r)2, assuming that most of the planetary mass is concentrated below the inner boundary at 0.97ro. The relevant nondimensional parameters are as follows: Rayleigh number Ra = αgoΔTd3/(κν), Ekman number E = ν/(ΩD2), and Prandtl number Pr = ν/κ, where α is the thermal expansion coefficient, go is gravity at the outer boundary, ν is viscosity, and κ is thermal diffusivity. These control parameters are set to E = 10−5, Ra = 3 × 107, 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 [; (50)], which uses the pseudospectral 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 divergenceless nature of the relevant quantities. The code uses the fast spherical harmonic library SHTns [; (51)]. The equations are time advanced using an explicit second-order Adams-Bashforth scheme for the Coriolis and the nonlinear terms, and an implicit Crank-Nicolson scheme for other terms (52).

We simulate the system for about 1500 rotations, which was enough to show a statistically stationary behavior. To limit the computational requirements, we used fourfold symmetry in the azimuthal direction (14, 53), 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 the azimuthal, latitudinal, and radial directions, 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 (14, 15, 25), 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 operatord()=1+D [+1hdmax+1hd]β(3)where D defines the amplitude of the function, ℓmax is the maximum spherical harmonic degree used in the simulation, ℓhd is the degree after which the hyerdiffusion starts, and β defines the rise of the function for degrees higher than ℓhd. For our simulation, we use D = 10, β = 6, and ℓhd = 450.

Thick shell case

In this simulation, we use the anelastic system of equations (54, 55) 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=x+x. The equations used in our simulation are in the entropy variable form with entropy contrast between the inner and the outer boundaries driving the convection (29). The magnetic field is scaled by ρoμoλiΩ, where ρo is the density on the outer boundary, μo is the magnetic permeability, and λi is the 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·(ρu)=0(4)(ut+u·u)=pρ2Eẑ×u+RaPrg s r̂+1Pmi E ρ(×B)×B+1ρ·S(5)where p′ is the pressure perturbation, g(r) is the gravity (assumed to vary linearly), s′ is the entropy perturbations, B is the magnetic field, and Pmi is ν/λi. The Rayleigh number Ra in this formulation is given by αogoTod3Δs(cpκν)−1. The traceless rate-of-strain tensor S is defined bySij=2ρ(eij13δij·u)with eij=12(uixj+ujxi)where δij is the identity matrix. The energy conservation gives time evolution equation for entropy perturbation sρT(st+u·s)=1Pr·(ρTs)+PrDiRaΦν+Pr Di λnormPmi2 E Ra(×B)2(6)where T is the background temperature and the viscous heating contribution is given byΦν=2ρ[eijeji13(·u)2](7)

The dissipation number Di is αogod/cp, where αo, go are the thermal expansivity and gravity on the outer boundary and cp is the specific heat at constant pressure. The evolution of the magnetic field is governed byBt=×(u×B)1Pmi×(λnorm×B)(8)where λnorm is the local magnetic diffusivity normalized by its value at the inner boundary ri. The magnetic field also follows the divergenceless condition ·B=0. We use the anelastic version (56) of the MagIC code to simulate this system of equations.

The inner boundary is at 0.1ro, and the outer is at ro. The control parameters for this simulation are as follows: E = 10−6, Ra = 4 × 109, Pr = 0.1, and Pmi = 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 (33): It remains constant until about 0.85ro, decreases by factors of 107 in a small radius, and stays constant afterward (see Profile 1 in fig. 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 (57) 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 107 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, β = 3, and ℓhd = 450 (see earlier section) for the higher-resolution grid and D = 5, β = 2, and ℓhd = 120 for the lower-resolution grid.

Correction (18 March 2021): The authors inadvertently forgot to cite a relevant study by Garcia et al. (2020) prior to publication. The second sentence in the second paragraph of the Discussion section has been modified to reference this work and make clear the relationship of the study to this article. The PDF and HTML (full text) versions have been updated.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: R.K.Y. thanks H. Cao for the interesting discussions. Funding: The work was supported by the NASA Juno project. The computing resources were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and by Research Computing, Faculty of Arts and Sciences, Harvard University. Author contributions: R.K.Y. designed and carried out the simulations. R.K.Y., M.H., and J.B. interpreted the results and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper. Additional data are available from the corresponding author R.K.Y. upon request.

Stay Connected to Science Advances

Navigate This Article