# Computer simulations of Jupiter’s deep internal dynamics help interpret what Juno sees

See allHide authors and affiliations

Contributed by Gary A. Glatzmaier, May 21, 2018 (sent for review February 13, 2018; reviewed by Chris Jones and Peter L. Olson)

## Significance

NASA’s Juno spacecraft is currently orbiting Jupiter, collecting high-fidelity measurements of its near surface. The resulting patterns of magnetic and gravity fields above Jupiter’s surface have already started providing insight into the type of winds in Jupiter’s deep interior. In addition, by monitoring cloud motions on Jupiter’s surface, Juno may detect a particular type of wave with properties that could tell us the size of Jupiter’s small rocky core. Here, we describe 3D time-dependent computer simulations of thermal convection and magnetic field generation in Jupiter’s deep interior that will help interpret and explain the maintenance of these global near-surface patterns that Juno continues to measure.

## Abstract

We describe computer simulations of thermal convection and magnetic field generation in Jupiter’s deep interior: that is, its convective dynamo. Results from three different simulations highlight the importance of including the dynamics in the very deep interior, although much of the convection and field generation seems to be confined to the upper part of the interior. A long-debated question is to what depth do Jupiter’s zonal winds extend below its surface. Our simulations suggest that, if global latitudinally banded patterns in Jupiter’s near-surface magnetic and gravity fields were detected by Juno, NASA’s orbiting spacecraft at Jupiter [Bolton S, et al. (2017) *Science* 356:821–825], they would provide evidence for Jupiter’s zonal winds extending deep below the surface. One of our simulations has also maintained, for a couple simulated years, a deep axisymmetric inertial wave, with properties at the surface that depend on the size of the model’s small rocky core. If such a wave was detected on Jupiter’s surface, its latitudes and oscillation frequency would provide evidence for the existence and size of Jupiter’s rocky core.

Jupiter’s latitudinally banded cloud pattern (Fig. 1) has intrigued people for hundreds of years. Ground-based and spacecraft observations have shown that these banded clouds are being sheared by latitudinally banded zonal winds (1) [that is, east–west jet streams alternating in latitude (2)]. NASA’s Galileo Probe, which was dropped from the Galileo Orbiter into Jupiter’s atmosphere near its equator in 1995, measured zonal wind speeds that increased from about 80 m / s at a pressure of about 1 bar (roughly Jupiter’s cloud top radius,

Jupiter has the most intense planetary magnetic field in our solar system. It is dominantly dipolar when measured far above its surface. However, small perturbations exist in this field (4) and in Jupiter’s gravity field (5) due to the planet’s internal dynamics (6). These perturbations, which are more prominent the closer the measurements are made to Jupiter’s surface, can reveal much about Jupiter’s deep interior.

Although many theoretical studies of Jupiter’s interior have been conducted using a variety of methods, many questions still remain about the structures of its deep fluid flows and magnetic fields and about the physical mechanisms that maintain them. Do Jupiter’s banded zonal winds at the surface extend deep enough below the surface to where electrical conductivity and zonal wind shear are high enough to generate a latitudinally banded magnetic field that could be detected by Juno, NASA’s orbiting spacecraft (7⇓–9)? Would strong zonal winds that extend to such a depth perturb enough mass to produce latitudinally banded perturbations in the gravity field that Juno could detect? The recent analyses (10⇓–12) of the gravity field from two of Juno’s perijove passes (i.e., its closest orbital approach to Jupiter, 5% of

## Model

Jupiter is a giant rotating planet of mainly hydrogen, a smaller amount of helium, and minor amounts of heavier elements (13). Fluid in the shallow atmospheric layer is approximately an ideal gas, but it compresses with depth into a dense liquid. However, as pressure increases with depth, the fluid becomes partially electron degenerate; that is, it becomes a “quantum fluid” with density and pressure being relatively insensitive to temperature. Experimental (14) and computational (15) studies indicate that Jupiter’s electrical conductivity increases rapidly with depth. The gas atmosphere at Jupiter’s surface is a poor conductor; however, as pressure increases with depth, diatomic hydrogen electrons become conducting electrons, which increase electrical conductivity until, at a radius of roughly 0.9

Computer models simulate thermal convection and magnetic field generation in Jupiter’s deep interior (16⇓–18) in attempts to explain the flows and fields observed at Jupiter’s surface and to predict their intensities and structures well below the surface. These models are defined by a coupled set of nonlinear partial differential equations based on the classical laws of conservation of mass, momentum, and energy that describe fluid flow, heat flow, and magnetic field generation in a rotating, self-gravitating, density-stratified fluid sphere as an analogue of Jupiter’s dynamic interior. Massively parallel supercomputers produce the 3D time-dependent numerical solutions to these equations. We have published (16, 19) detailed descriptions of our basic computer model used to produce the computer simulations presented here. The set of equations and the model prescriptions are presented in *Methods*. However, before discussing the results, we briefly mention some of the model details that define these particular Jupiter simulations.

We fit our 1D reference-state density (Fig. 2, *Left*), pressure, and temperature to a 1D evolutionary interior model (5, 15) using a polytropic relationship of pressure proportional to density squared (20). The 1D evolutionary model solves only for the spherically averaged radial profiles of the thermodynamic variables and simply parametrizes the effects of 3D convection. However, this 1D model can afford to incorporate a more realistic equation of state with a more detailed treatment of chemical composition and radiation transfer than can current 3D global convective dynamo models. It can also afford the much higher resolution in radius required to model the very shallow atmospheric layer, which current dynamo models ignore. Our code then numerically solves the system of equations (*Methods*) that describes the 3D time-dependent perturbations relative to our 1D time-independent rotating reference state [that is, the thermodynamic and gravitational perturbations and the fluid velocity and magnetic field vectors (*Methods*) (16)].

Our model prescriptions for the (dimensional) planetary mass, radius, and rotation rate are set to Jupiter’s. Based on the calculations of French et al. (15), we choose constant values for the viscous and thermal diffusivities. However, the estimated values of these molecular diffusivities would result in a cascade of energy down to eddies too small to resolve with any 3D global simulation of convection in Jupiter’s interior. Moreover, this unresolved small-scale turbulence transports much more heat and momentum than molecular diffusion. Therefore, we numerically resolve as much of the energy spectrum as we practically can and then parametrize the transport of heat and momentum by the remaining unresolved eddies as a diffusion process using enhanced values for the (turbulent) viscous and thermal diffusivities. However, to obtain peak amplitudes of the winds and fields at the surface that are also Jupiter like, we have to compensate for these enhanced diffusivities by driving thermal convection with a luminosity that is larger than Jupiter’s.

Jupiter’s magnetic diffusivity is inversely proportional to its electrical conductivity and is larger than its viscous and thermal diffusivities. We prescribe an electrical conductivity for the model that increases (approximately) exponentially (21) by four orders of magnitude from the model’s surface down to 0.9 *Right*) it transitions to a constant in the deeper interior. It is important for the simulated dynamics to have the amplitudes of these three model diffusivities be in the same order as they are for Jupiter. Therefore, since we are forced to prescribe enhanced values for the viscous and thermal diffusivities, the maximum electrical conductivity that we prescribe below 0.9 *Methods* and Fig. 2, *Right*) is several orders of magnitude less than predicted values (14, 15).

A required feature of our study is the self-consistent calculation (in 3D and at every numerical time step) of the perturbation in the gravitational potential by solving a Poisson equation (Eq. **3**) forced by the local density perturbation. We also monitor the local mass flux in and out of the model’s (fixed but permeable) spherical upper boundary and approximate the resulting local time-dependent accumulation or depletion of mass above this boundary as a surface mass density perturbation in hydrostatic balance on this upper boundary surface. Then, assuming no other density perturbations above this boundary, this procedure allows us to easily calculate the gravity perturbations anywhere above the upper boundary as a function of time. This provides a connection between centrifugal accelerations due to, for example, the surface zonal winds and gravity perturbations at Jupiter’s perijove. A similar method (16) allows us to calculate the magnetic field at perijove, assuming that electric currents are negligible above the surface.

The simulations presented here were initially run with a spatial resolution in grid space (radial × latitudinal × longitudinal levels) of 241 × 384 × 768 and then after reaching a statistical equilibrium, continued at 241 × 768 × 1,536. In spectral space, all spherical harmonics are used up to degree and order 255 initially and later 511 (for the latitudinal and longitudinal expansions) and in Chebyshev polynomials up to degree 241 (for the radial expansions). The typical numerical time step is 100 (simulated) seconds based on a 10-h rotation period. Each of the simulations presented here has run at least 20 simulated years (i.e., about 6 million numerical time steps). The details of this computational procedure are described in ref. 16.

We discuss three cases that illustrate how the prescribed depth of the model’s fluid interior affects the patterns and intensities of the winds and fields at the surface, which are compared with those observed on Jupiter to estimate how realistic the simulated interior dynamics may be. In addition, different mean temperature profiles in radius are forced in these three cases via thermal boundary conditions and prescribed internal heating that determine if and where thermal convection occurs. Case 1 simulates thermal convection and magnetic field generation in a deep rotating density-stratified fluid shell with an upper boundary radius of 0.98

## Case 1: Convective Dynamo in a Deep Shell

Rotating thermal convection produces bands of longitudinal flow that alternate eastward and westward with latitude relative to the rotating frame of reference (Fig. 3, *Left*). After roughly 3 million numerical time steps, the average in longitude of this flow, the zonal wind, becomes virtually steady in time (Fig. 3, *Right*). This is also called differential rotation, since in the inertia frame of reference, it is a radially and latitudinally dependent rotation rate. It is maintained in our simulations by the redistribution of angular momentum by (nonlinear) advective transport, which naturally occurs in a rotating convection zone when vortices are generated by Coriolis forces as they rise (expand) or sink (contract) through the density-stratified interior (16, 22). Some recent studies (11, 12), however, are based on the assumption that Jupiter’s differential rotation is maintained as a “thermal wind.” That is, instead of solving the full set of equations (*Methods*), they consider only two terms, the buoyancy and Coriolis forces in Eq. **5**, and assume that the curls of these exactly balance. This is not a bad approximation for the Earth’s shallow atmosphere driven by solar insolation, but our deep convective dynamo simulations of Jupiter do not indicate such a balance. However, it is encouraging that these two studies (11, 12), which fit different types of models to the early Juno gravity data, agree that strong fluid flows exist down to at least 4% of

The zonal wind for Case 1 (Fig. 3, *Right*) is (very nearly) only a function of the distance from the planetary rotation axis; that is, the amplitude of this zonal wind is “constant on cylinders” that are coaxial with the rotation axis. The total kinetic energy of the zonal wind, integrated over the entire convection zone, is at least an order of magnitude greater than the total kinetic energy of the 3D convection, which maintains the zonal wind. The local kinetic energy density of convection peaks in the upper part of the interior where the temperature gradient is the most unstable (superadiabatic) and the relative drop in density with radius is greatest. However, for this Case 1, convection is significant enough throughout the entire fluid shell to maintain zonal winds constant on cylinders down to the model’s rocky core (Fig. 3, *Right*).

The surface manifestation of this deep zonal wind is an eastward equatorial jet, with a maximum amplitude of about 120 m / s, slightly greater than Jupiter’s maximum eastward wind at cloud tops (about 100 m / s). The maximum convective velocity near the surface is about 10 m / s. However, instead of multiple, alternating, narrow zonal jets at higher latitude as observed on Jupiter’s surface (Fig. 1), only one latitudinally broad westward jet exists at middle latitude and one broad eastward jet exists at high latitude in each hemisphere.

Note that, long before this steady zonal wind structure sets in, banded zonal winds first develop in a shallow layer just below the surface, where the steep temperature gradient there drives vigorous convection. A strong steady eastward jet develops in the equatorial region. At higher latitude, multiple narrow alternating jets develop in this shallow layer that are not constant on cylinders and do not extend from the northern to southern hemispheres. However, these shallow high-latitude jets are relatively weak and very intermittent, unlike those observed on Jupiter. Eventually, this pattern evolves into a deep, slowly changing, approximately constant on cylinders pattern of differential rotation that then takes a couple million more time steps to evolve into the steady, cylindrically symmetric structure that extends throughout the deep interior (Fig. 3).

These convective flows and zonal winds also generate magnetic field in the electrically conducting fluid deep below the surface, and as seen in Fig. 4, this field extends up through the electrically insulating surface layer, where the field is unaffected by fluid flows. For example, Fig. 5, *Right* shows the latitudinally banded magnetic field seen in Fig. 4 in a spherical surface at a distance above the model’s upper boundary equal to Juno’s closest orbital approach (perijove) to Jupiter’s cloud surface. At a depth where electrical conductivity is large, a pair of deep zonal winds in opposite directions would shear poloidal (vertical and north–south directed) field into toroidal (east–west directed) field. Also, convective vortices twist toroidal field into latitudinal bands of poloidal magnetic field loops. This “dynamo” process converts kinetic energy of the fluid flow into magnetic energy. The greatest magnetic energy density occurs where the product of the fluid flow amplitude (which decreases with depth) and the electrical conductivity (which rapidly increases with depth) is maximum; this occurs in Case 1 near 0.9

To limit the distortion of the poloidal field and the resulting ohmic heating within the deep electrically conducting interior, poloidal field tends to align itself on surfaces of constant angular velocity (in this case, cylindrical surfaces of constant zonal wind speed) as suggested by Ferraro’s isorotation law (23). The resulting banded (poloidal) magnetic structure, formed deep below the surface, is retained in the magnetic field that extends up through the insulating upper layer and out to Juno’s perijove. A similar magnetic effect has been produced (24) with a prescribed differential rotation using a 2D time-independent kinematic mean-field model of the semiconducting region above 0.9 *Right* would be virtually nonexistent along most of Juno’s highly elliptic (53-d) polar orbit far from Jupiter, where the field appears dominantly dipolar.

Other 3D (nearly full-sphere) convective dynamo simulation studies of Jupiter (17, 18) have not developed latitudinally banded magnetic structures like those in Fig. 5, *Right*. Instead, their surface magnetic fields are dominantly dipolar with differential rotation in the deep interior suppressed by magnetic drag. The differences between their results and the banded results presented here could be due to the different choices for the model parameters, boundary conditions, internal heating, centrifugal force, or our relatively weak electrical conductivity below 0.9

These same fluid flows determine the 3D distribution of local low and high mass densities, which slightly perturb the nearly spherically symmetric (reference state) gravity field below and above the surface. Consider the radial component of the gravity perturbation at Juno’s perijove in Fig. 5, *Upper Left*. Blue bands (downward-directed gravity perturbations) in this image correspond to faster rotating (eastward) jets at the surface (Fig. 3), which due to greater centrifugal forces, slightly elevate mass through our permeable upper boundary, enhancing the (downward-directed) gravitational acceleration at perijove. Juno data from two perijove passes (10) already show evidence of local gravity variations banded in latitude a few times smaller in amplitude compared with those in Fig. 5, *Upper Left*. When Juno data from many more than two orbital passes are analyzed, we will see if Jupiter’s near-surface gravity field is globally banded.

Spherical harmonic spectra of the simulated magnetic field patterns (measured at perijove) provide additional information to be compared with Juno’s data. Fig. 6, *Upper* is a plot of the amplitudes of the axisymmetric spherical harmonic (Gauss) coefficients,

Fig. 6, *Lower* is a plot of the amplitudes of the axisymmetric spherical harmonic coefficients,

Juno is able to detect these spherical harmonic gravity coefficients down to an amplitude of roughly *Lower*) are slightly smaller than those reported in ref. 10.

## Case 2: Convective Dynamo in a Shallow Shell

Here, we examine an approximation similar to what has been made in many previously published global simulations of Jupiter (25). Instead of simulating the fluid dynamics down to a small rocky core, an impermeable lower boundary of the convection zone is prescribed at a radius of 0.8 *Right*) is not severely suppressed in the deep interior, likely because (as mentioned above) the magnetic field lines there tend to be in surfaces of constant angular velocity, which reduces magnetic drag and ohmic heating.

The radial profiles for Case 2 reference state density, temperature, pressure, and electrical conductivity are the same as those in the upper region of Case 1. Thermal boundary conditions are prescribed for Case 2 that drive convection throughout this (shallower) fluid shell. However, since Jupiter’s rocky core would be much deeper than this artificial lower boundary, here (unlike for Case 1) we apply viscous stress-free and magnetic stress-free conditions on this lower boundary.

As for Case 1, a shallow zonal wind layer initially develops in the upper part of this fluid shell, with multiple narrow weak intermittent jets at high latitude. However, after a couple million more numerical time steps, the zonal wind evolves into a steady constant on cylinders structure throughout the fluid shell (Fig. 7). Several high-latitude jets are maintained on the surface, a pattern slightly more like Jupiter’s (Fig. 1) than that for Case 1. The maximum surface zonal wind speed occurs in the eastward equatorial jet with an amplitude of 150 m / s, which is 50% higher than the maximum measured at Jupiter’s cloud tops but less than what the Galileo probe measured on Jupiter below the cloud tops. As for Case 1, the total kinetic energy in these Case 2 zonal winds is typically an order of magnitude greater than that of the convection. The magnetic and gravity fields at perijove for this case also have latitudinally banded patterns, although not as well-defined as those for Case 1. An earlier simulation of Saturn’s dynamics (16) in its outer region produced latitudinally banded magnetic fields qualitatively similar to Case 2. Although Case 2 maintains surface winds that better match those of Jupiter, the flows, fields, and thermodynamics that should exist below the (artificial) impermeable lower boundary at 0.8

## Case 3: Convective Dynamo in a Shallow Shell Above a Deep Convectively Stable Interior

As discussed for Case 2, a more Jupiter-like surface zonal wind pattern can be simulated when an impermeable lower boundary to the convection zone is imposed at radius 0.8

Case 3 is a simple test simulation of such a scenario. Convectively stable fluid regions are maintained from the lower boundary at radius 0.10

The model maintains these three regimes using a time-dependent and radially dependent heat source that continually “nudges” the mean entropy perturbation toward a prescribed target profile (16). This heat source is proportional to the difference between the current (local) mean entropy and the target value and inversely proportional to a time constant that determines how closely the entropy perturbation tracks the target profile. For Case 3, the target entropy increases linearly with radius in the lower (stable) region by a total of 250 J / (kg K), it decreases linearly with radius in the middle (convectively unstable) region by 22 J / (kg K), and it increases linearly with radius in the upper (stable) region by 22 J / (kg K). These changes in entropy are relatively small compared with the model’s specific heat capacity at constant pressure (

Downwelling plumes in the convecting region (between 0.90 and 0.96

As in the two previous cases, the eastward zonal jet in the equatorial region for Case 3, seen in Fig. 8, *Top*, is maintained by the nonlinear convergence of angular momentum flux that is driven by 3D rotating convection within the density-stratified unstable region. Weak, relatively shallow zonal wind jets, mainly westward directed relative to the rotating frame of reference, also exist at high latitude. In addition, for Case 3, an axisymmetric inertial wave developed (Fig. 8) at midlatitude in both hemispheres throughout the interior. This snapshot illustrates the ray-like pattern of shear layers in the axisymmetric longitudinal and radial components of fluid velocity on conical surfaces. Although several waves like this but propagating at different angles were initially excited, this one survived, because it is continually reinforced by closing on itself after “reflecting” off the upper boundary at four latitudes and being tangent to the lower (rocky core) boundary at two latitudes. A dominant inertial wave would likely be difficult to maintain in a full sphere with no rocky core to force nodes in the transverse oscillation. All legs of this wave oscillate at the same frequency and propagate at the same angle relative to the planetary rotation axis (16, 26, 27). The latitudinal separation between the two rays at the surface in each hemisphere equals the diameter of the model’s rocky core. If a similar oscillating axisymmetric pattern of zonal and radial flows was observed on Jupiter’s surface by, for example, time-lapse imagery from Junocam (28), its frequency and surface latitudes could determine the existence and size of the planet’s rocky core.

This inertial wave (Fig. 8) was a dominant feature for at least 2 simulated years. Its period was comparable with the 10-h rotation period, roughly six times less than periods of long-wavelength internal gravity waves in this simulation, and much less than characteristic periods of the simulated Alfvén waves. However, maintaining such a large-scale dominant inertial wave in Jupiter could be difficult if Jupiter’s rocky core boundary was not well-defined. Likewise, reflections at Jupiter’s surface would be more complicated and less efficient than those off our model’s upper boundary. Even without these complications, the inertial wave in the Case 3 simulation eventually disappeared and has not returned after an additional 13 simulated years, likely because of interference with convective flows or short period gravity waves.

The shearing and twisting fluid flows in Case 3 self-consistently generate somewhat of a globally banded magnetic field (Fig. 8, *Bottom*). The convective dynamo mechanism is most efficient near 0.9

Note that the dipolar part of the magnetic field in Case 3 (check the polar regions in Fig. 8, *Bottom*) has the reversed polarity of that in Case 1 (Figs. 4 and 5) and likewise, the reverse of Jupiter’s dipole polarity. This is not an issue because the resulting polarity in these simulations depends on the arbitrary initial conditions and because the exactly reversed polarity would also satisfy the coupled set of differential equations for exactly the same flow and thermal patterns. Note also that the magnetic dipole part of the field in Case 1 is nearly an axial dipole; that is, the magnetic dipole axis is nearly the same as the planetary rotation axis. However, the dipole axis for Case 3 is a little more tilted (somewhat more like Jupiter’s) and is a little more time dependent than that for Case 1. This greater variability in space and time for Case 3 occurs because the zonal wind structure is not constant on cylinders throughout the entire fluid interior; the dominant dynamic length scales are more confined to the shallow convecting region and therefore have less global influence. No dipole reversals have occurred in these Jupiter simulations; however, as mentioned, the simulations span only about 20 simulated years.

## Discussion

This project has two main objectives. One is to investigate the importance of including the entire fluid interior in a model when studying fluid flows and magnetic field generation in Jupiter. The three cases presented here illustrate how the dynamics in Jupiter’s deep interior can affect the flows and fields near Jupiter’s surface. The simulated zonal winds in all three cases are considered “deep” zonal winds maintained by the dynamics of the deep fluid interior, not zonal winds maintained by solar insolation in Jupiter’s very shallow gaseous atmosphere. The early results from Juno (10) at least support the existence of local winds that are deep and strong. Cases 2 and 3 maintain somewhat more Jupiter-like zonal winds at the surface than does Case 1 by completely ignoring the very deep interior or by suppressing convection there, respectively. Case 3 crudely explores such a situation by forcing, via internal heating, a subadiabatic (convectively stable) fluid interior below a relatively shallow upper convecting region that excites internal gravity waves in the lower stable region and Alfvén waves and inertial waves throughout the entire fluid interior.

However, as discussed, our computer model is far from being perfect. No one has been able to produce a 3D global convective dynamo simulation with all model prescriptions being Jupiter-like because of the huge computational resources that would be required. Although we set the planet size, mass, and rotation rate to Jupiter values, we are forced to use enhanced values for the diffusivities, which require a luminosity greater than Jupiter’s to evolve the surface amplitudes of zonal winds and magnetic fields to Jupiter-like values. However, encouraged by Case 3 results, the next step in this investigation could be to continue Case 3, gradually decreasing these diffusivities while increasing spatial and temporal resolutions. Although one could not afford to extend the simulation very far in time, it might be enough to provide useful insight.

Another step in this investigation should be the development of a more physical model of the deep interior below 0.9

Our other objective for this project is to provide interpretations of certain global patterns that the Juno mission might discover at Jupiter. The simulations presented here suggest that a detection by Juno of a global latitudinally banded pattern in the magnetic field would provide evidence for Jupiter’s banded zonal winds extending deep below the cloud surface to where electrical conductivity is sufficiently high for these winds to shear the deep magnetic field enough to be detectable out to Juno’s perijove. Similarly, a detection by Juno of a global latitudinally banded pattern in the gravity field perturbations would argue for deep banded zonal winds that shear a sufficient amount of mass to produce detectable gravity variations at Juno’s perijove. The recent analyses of Juno’s near-surface gravity from two perijove passes (10⇓–12) provide encouraging evidence for strong fluid flow well below Jupiter’s cloud surface. Soon, Juno will have gathered sufficient global coverage of Jupiter’s near surface to reveal how globally banded its magnetic and gravity fields are. In addition, although much less likely, a detection by Junocam of an axisymmetric inertial wave, manifested at the surface as a pair of zonal wind jets oscillating in sync with radial flows and occurring at the same latitudes in both hemispheres, could provide exciting evidence for the existence and size of Jupiter’s rocky core.

## Methods

The details of the anelastic approximation, boundary conditions, numerical method, parallel programing, and computer graphics are described in refs. 16 and 19. Here, we briefly described the set of equations that are solved and the model details for the cases presented in this paper:

This coupled system of nonlinear partial differential equations describes mass conservation (Eq. **1**), magnetic flux conservation (Eq. **2**), the gravitational potential perturbation (Eq. **3**), the perturbation equation of state (Eq. **4**), momentum conservation (Eq. **5**), magnetic induction (Eq. **6**), and energy conservation (Eq. **7**). The numerical solution to these equations updates, at each numerical time step, the 3D time-dependent fluid flow v, magnetic field B, and perturbations in density ρ, pressure p, specific entropy S, and gravitational potential U. The rate of strain tensor is

For Cases 1 and 2, the radial gradient of the spherically symmetric part of the entropy perturbation is determined by applying a zero radial gradient lower boundary condition on the entropy perturbation and a fixed entropy perturbation condition on the upper boundary. Convection is driven by a prescribed internal heating source (**7**), which represents the slow cooling rate of the planet (17, 32). For these two cases, Q is a constant in space and time. For Case 3, both the lower and upper thermal boundary conditions are fixed entropy, equal to the target entropy values at these boundaries. The internal heating source for this case constantly nudges the spherically averaged perturbation entropy toward the target profile as described above.

The pressure perturbation part of the buoyancy term in the momentum equation (Eq. **5**) has been combined with the gradient of the pressure perturbation using the Lantz–Braginsky–Roberts formulation within the anelastic approximation (33, 34). Here,

For the simulations presented here, the total mass of the model planet is ^{6} and *Right*), is

## Acknowledgments

We thank C. Jones and P. Olson for their helpful suggestions for improving this manuscript. Support for this project was provided by NASA Grant OPR NNX13AK94G. Computational resources were provided by the NASA Ames Research Center and by National Science Foundation funds that supported parallel computing at the University of California Santa Cruz.

## Footnotes

- ↵
^{1}Email: glatz{at}es.ucsc.edu.

This contribution is part of the special series of Inaugural Articles by members of the National Academy of Sciences elected in 2010.

Author contributions: G.A.G. designed research, performed research, analyzed data, wrote the computer code, and wrote the paper.

Reviewers: C.J., University of Leeds; and P.L.O., Johns Hopkins University.

Conflict of interest statement: G.A.G. was one of 37 coauthors, including reviewers Chris Jones and Peter L. Olson, on a 2016 report of a computational performance benchmark exercise. They were not research collaborators on that report.

See QnAs on page 6880.

- Copyright © 2018 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- Porco C, et al.

- ↵
- Ingersoll A

- ↵
- Atkinson D,
- Pollack J,
- Seiff A

- ↵
- Connerney J,
- Acuna M,
- Ness N,
- Satoh T

- ↵
- ↵
- ↵
- ↵
- Jones C,
- Holme R

- ↵
- Bolton S, et al.

- ↵
- ↵
- ↵
- ↵
- Bagenal F,
- McKinnon W,
- Dowling T

- Guillot T,
- Stevenson D,
- Hubbard W,
- Saumon D

- ↵
- ↵
- French M, et al.

- ↵
- Glatzmaier G

- ↵
- Jones C

- ↵
- ↵
- ↵
- Hubbard W

- ↵
- ↵
- Glatzmaier G,
- Gilman P

- ↵
- ↵
- ↵
- Gastine T,
- Heimpel M,
- Wicht J

- ↵
- ↵
- Seelig T,
- Harlander U

- ↵
- Hansen C, et al.

- ↵
- Stevenson D

- ↵
- Wahl S, et al.

- ↵
- Moll R,
- Garaud P,
- Stellmach S

- ↵
- ↵
- Lantz S

- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Earth, Atmospheric, and Planetary Sciences

### See related content:

- QnAs with Gary A. Glatzmaier- Jun 25, 2018