Computer simulations of Jupiter's deep internal dynamics help interpret what Juno sees.

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, m) to a maximum of about 170 m / s at a pressure of 4 bars, below which the wind speed remained approximately constant until the probe’s signal was lost at about 21 bars (3). Although this suggests that the zonal winds on Jupiter extend below the visible surface, the 21-bar depth is less than 1% of below the surface.

Fig. 1.

(Left) An image of Jupiter in natural color taken by the Hubble Space Telescope [Image courtesy of NASA, ESA, and A. Simon (Goddard Space Flight Center)]. (Right) Jupiter’s surface zonal wind speed in meters per second relative to its (rotating) magnetic field (1) (Reprinted from ref. 1, with permission from AAAS).

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 above the surface) have already produced exciting results suggesting that strong fluid flows exist down to at least 4% of below Jupiter’s cloud surface. Providing detailed physical descriptions of the structure and maintenance of Jupiter’s zonal winds and its magnetic and gravity fields is a major goal of computational studies of Jupiter’s convective dynamo.

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 , hydrogen becomes a liquid metal (14). At the center of Jupiter, there may be a small core of rocky material with a radius of roughly 0.1 (13).

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 . 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 () 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 () (16)].

Fig. 2.

The model’s prescribed reference-state density (Left) and electrical conductivity (Right).

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 , where (Fig. 2, 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 ( 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. ) 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 and a lower boundary radius of 0.10 . Below this model’s lower boundary is a solid, electrically conducting, rocky core that rotates in response to the viscous and magnetic stresses on it and according to its Jupiter-like moment of inertia. Internal heating maintains a mean temperature profile in the radius that is convectively unstable; that is, convection occurs throughout the spherical fluid shell for Case 1. The lower boundary for Case 2 is defined as an impermeable stress-free spherical interface at a radius of 0.8 . Thermal convection is also driven throughout this (artificially shallower) spherical fluid shell by internal heating. Case 3 uses the same rocky core as that of Case 1; that is, the spherical fluid shell extends from radius 0.10 to 0.98 . However, its mean temperature profile is constantly “nudged” by internal heating toward a target mean temperature profile that is convectively stable below 0.90 and above 0.96 and convectively unstable between these two radii.

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 (), they consider only two terms, the buoyancy and Coriolis forces in Eq. , 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 below Jupiter’s surface.

Fig. 3.

A snapshot of the longitudinal component of the simulated fluid velocity (i.e., winds) for Case 1. (Left) The longitudinal winds plotted in an equal area projection of the entire spherical upper surface of the fluid model. (Right) The “zonal wind” (i.e., longitudinal average of the longitudinal winds) plotted in a meridian plane. Yellows and reds represent eastward-directed winds; blues represent westward-directed winds (relative to the rotating frame of reference).

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 . Note also that, for this particular simulation, the magnetic dipole axis is tilted (relative to the planetary rotation axis) much less than it is on Jupiter. This is probably due to the constant on cylinders zonal wind structure that dominates the entire fluid interior.

Fig. 4.

A snapshot of the simulated magnetic field for Case 1 plotted in 3D as magnetic lines of force. The axis of planetary rotation is centered and vertical. The model’s insulating near-surface region starts at the depth where the east–west magnetic lines begin to bend upward. Lines are gold where the field is outward directed and blue where the field is inward directed.

Fig. 5.

A snapshot of the simulated gravity field perturbation (Left) and magnetic field (Right) for Case 1 plotted in a spherical surface with a radius equal to Juno’s perijove. Upper shows the radial components of the fields (yellows and reds are outward directed, and blues are inward); Lower shows the colatitudinal components (yellows and reds are southward directed, and blues are northward). The peak amplitudes of the gravity perturbation for the radial and colatitudinal components are 18 and 14 mgal, respectively. The peak amplitudes of the magnetic field for the radial and colatitudinal components are 5 and 4 gauss, respectively.

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 . Small magnetic length scales above the surface decay more rapidly with distance from Jupiter than the dipole mode (degree ); therefore, the banded magnetic patterns seen in Fig. 5, 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 . Analysis of Juno’s first perijove pass (7) shows two intense magnetic flux patches on opposite sides of the equator indicating that Jupiter’s near-surface magnetic field is more complicated than previously expected. It will be interesting to see, after more Juno perijove passes have been analyzed, if Jupiter’s near-surface magnetic field pattern is globally banded.

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, , as a function of spherical harmonic degree, , for the Fig. 5 snapshot. Our simulated data are calculated every numerical time step for all spherical harmonic degrees and orders up to 511 but plotted here (for order 0) out to degree 30. They show how the intensity of the longitudinally averaged part of the simulated magnetic field decreases with increasing spherical harmonic degree (that is, as the length scale of the mode decreases). The relatively broad spectrum at low degrees reflects the relatively large contribution from the banded (nondipolar) part of the field (Fig. 4). However, this simulated surface field is somewhat more intense than Jupiter’s; for example, the dipole part, , for Case 1 is four times larger than Jupiter’s (4). Of course, even if the model was near perfect, we would only expect qualitative agreement with Jupiter’s magnetic spectra, because Jupiter and the simulation are both time dependent.

Fig. 6.

A snapshot of the longitudinally averaged magnetic and gravity fields as a function of spherical harmonic degree measured at what would be perijove for Juno’s orbits around Jupiter. (Upper) The amplitude of the magnetic Gauss coefficients in units of one billionth of a Tesla (i.e., ). (Lower) The amplitude of the perturbation in the gravity field relative to the spherically averaged value. The relatively straight line outlines the (discrete) even degree values that represent just the effect of the planet’s oblateness (due to the planet’s rotation) (6). The jagged curve outlines the (discrete) degree values that represent the remaining effects due to the simulated fluid flows within the modeled planet.

Fig. 6, Lower is a plot of the amplitudes of the axisymmetric spherical harmonic coefficients, , for the nonspherically symmetric part of the simulated gravity field (measured at perijove) as a function of spherical harmonic degree . It shows that this small aspherical perturbation in the simulated gravity field is dominated, for even degrees up to 10, by the planetary oblateness (due to centrifugal force) (6). For degrees higher than 10 (i.e., smaller length scales), fluid flows, especially the strong surface zonal winds, have a greater effect (than oblateness) on the structure of the aspherical part of the gravity field near the surface. Note that, because oblateness is symmetric across the equator, it makes no odd degree contribution to the gravity field. Therefore, the contributions from odd degrees, including those less than degree 10, are totally due to deep axisymmetric fluid flows that are asymmetric with respect to the equator.

Juno is able to detect these spherical harmonic gravity coefficients down to an amplitude of roughly (6, 10). An exciting recent analysis (10) of two of Juno’s near-surface perijove passes has already produced estimates of the axisymmetric spherical harmonic coefficients of Jupiter’s gravity field out to degree 12. In addition to the even degrees due to Jupiter’s oblateness, odd degrees have been measured indicating that strong wind jets exist on Jupiter well below the depth to which the Galileo probe measured. The amplitudes of the corresponding odd degrees produced in Case 1 (Fig. 6, 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 . The main reason for ignoring convection below this artificial lower boundary is to focus all of one’s computational resources in the region where magnetic field generation is expected to be the most efficient, which as in Case 1, is the upper part of the convection zone. Another reason has been based on the assumption that differential rotation in Jupiter is severely suppressed by magnetic drag below a depth where the electrical conductivity is sufficiently large. This is supported by 2D models of Jupiter’s interior fitted to the recent Juno gravity data (11, 12). It is also supported by some deep 3D dynamo simulations of Jupiter (17, 18) that maintain sufficiently vigorous convection and consequently intense magnetic fields. However, the evolved differential rotation in Case 1 (Fig. 3, 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 have been ignored.

Fig. 7.

Snapshots of the longitudinal wind for Case 2 in the same projections as those in Fig. 3 for Case 1.

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 (or above) instead of at 0.1 , which is what 1D evolutionary models predict. There may be important physical mechanisms operating in Jupiter, but missing in our model (in addition to higher electrical conductivity), that tend to confine convection to a relatively shallow upper region, while allowing only less energetic motions in the deeper part of the fluid interior. Therefore, we ask what modifications to the model could produce a more Jupiter-like latitudinally banded pattern at the surface with multiple alternating zonal wind jets up to high latitude while simulating the internal magnetohydrodynamics down to radius 0.1 .

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 up to 0.90 and from 0.96 up to the upper boundary at 0.98 ; a convectively unstable region is maintained between these two stable regions. Internal gravity waves are excited in the lower and upper (stable) regions, where the mean temperature decreases more slowly with radius than an adiabatic temperature profile would; that is, in these two regions, the mean temperature profiles are “subadiabatic,” and the mean (specific) entropy increases with radius. Thermal convection occurs in the middle (unstable) region, where the mean temperature decreases faster with radius than an adiabatic temperature profile would; that is, the mean temperature profile is “superadiabatic,” and the mean entropy decreases with radius. The shallow upper layer serves as a crude stratosphere, which provides some resistance to convective penetration from below, albeit without realistic atmospheric physics. Our tests show that this upper layer has little effect on the large-scale flow patterns at the surface, likely because this layer is not stable enough.

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 ( J / (kg K)) as required for codes like ours that use the “anelastic approximation” to the fully compressible system of equations (16).

Downwelling plumes in the convecting region (between 0.90 and 0.96 ) penetrate slightly into the lower stable region exciting 3D “internal gravity waves” that propagate through the interior due to buoyancy restoring forces. Magnetic fields, distorted by fluid flows, provide magnetic restoring forces that produce “Alfvén waves.” Coriolis forces, which occur (in the rotating frame of reference) when fluid flows perpendicular to the axis of rotation, serve as yet another restoring force that drive “inertial waves.” Alfvén waves and inertial waves exist in both the convectively stable and unstable regions; internal gravity waves get excited only in the convectively stable regions. All three types of waves are transverse waves; that is, the oscillating fluid motion is perpendicular to the direction of phase propagation. Internal gravity wave frequency is maximum when the wave propagates horizontally (i.e., when oscillating fluid motions are parallel to gravity). Alfvén wave frequency is maximum when phase propagation is parallel to the local magnetic field direction (i.e., when oscillating fluid motions are perpendicular to the local magnetic field). Inertial wave frequency is maximum when the wave is propagating parallel to the rotation axis (i.e., when oscillating fluid motions are perpendicular to the rotation axis). When all three types of waves are excited with comparable frequencies, the resulting wave motion can be complicated.

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.

Fig. 8.

A snapshot of the longitudinal and radial components of fluid velocity and the radial component of magnetic field for Case 3. (Left) Plotted in a constant radius surface slightly below the upper boundary. (Right) Longitudinally averaged and plotted in a meridian plane. Reds and yellows represent eastward (Top) or upward (Middle and Bottom); blues represent westward or downward.

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 , where convective eddies and zonal winds are still significant and electrical conductivity is at its maximum. The peak magnetic field intensity, about 500 gauss, occurs near this radius. The wave motion in the deep (stable) interior, which is excited by the convection above, is less energetic and less helical, and therefore, it is less efficient in generating magnetic field. We see no significant “pumping” of the magnetic field from the convecting region into the lower stable region by convective downwellings. However, as mentioned, the model’s prescribed electrical conductivity is much less than predicted for the deep interior of Jupiter (14, 15).

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 . For example, a plausible thermal mechanism for maintaining a convectively stable layer deep within Jupiter could be absorption of latent heat by molecular hydrogen when, while sinking, it dissociates into monatomic hydrogen. Likewise, latent heat is emitted in rising monatomic hydrogen as it combines into molecular hydrogen. This phase change, which is thought to continually occur in Jupiter between roughly 0.8 and 0.9 (14), should keep the temperature nearly constant in this layer and possibly subadiabatic. It has also been proposed that “helium rainout” (29, 30) may be occurring in this same region. That is, in this region, helium becomes insoluble in hydrogen and gravitationally falls through the fluid, which produces a stable layer with heavy fluid below light fluid. There may also be double-diffusive instabilities occurring below this layer due to very different molecular diffusion rates of temperature and heavy element composition in a region with a destabilizing (superadiabatic) thermal stratification and a stabilizing mean molecular weight stratification (31).

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:Each dependent variable in these equations is written as the sum of a prescribed time-independent (hydrostatic and nonmagnetic) reference state that only depends on radius (the “barred” variables in Eqs. 1–7) plus a small (unbarred) perturbation relative to its reference state value that depends on time and 3D space.

This coupled system of nonlinear partial differential equations describes mass conservation (Eq. ), magnetic flux conservation (Eq. ), the gravitational potential perturbation (Eq. ), the perturbation equation of state (Eq. ), momentum conservation (Eq. ), magnetic induction (Eq. ), and energy conservation (Eq. ). The numerical solution to these equations updates, at each numerical time step, the 3D time-dependent fluid flow , magnetic field , and perturbations in density , pressure , specific entropy , and gravitational potential . The rate of strain tensor is ; is the magnetic permeability; and are the turbulent viscous and thermal diffusivities, respectively; is the magnetic diffusivity; is the unit tensor; and and are unit vectors in the radial and colatitudinal directions, respectively.

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 ( in Eq. ), which represents the slow cooling rate of the planet (17, 32). For these two cases, 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. ) has been combined with the gradient of the pressure perturbation using the Lantz–Braginsky–Roberts formulation within the anelastic approximation (33, 34). Here, is the reference state gravitational acceleration, and the term is centrifugal acceleration. Near Jupiter’s surface at low latitude, the latter is nearly one-tenth of the former and therefore, is not neglected.

For the simulations presented here, the total mass of the model planet is , and the average planetary rotation rate is (that is, a 10-h rotation period). The reference state density, temperature, and pressure at the lower boundary are , , and , respectively. The viscous and thermal diffusivities are 106 and , respectively. Magnetic diffusivity, , where is electrical conductivity (Fig. 2, Right), is at the lower boundary and at the upper boundary.