Geodynamo theory and simulations
Paul H. Roberts
Institute of Geophysics and Planetary Physics, University of California, Los Angeles, California 90095
Gary A. Glatzmaier
Department of Earth Sciences, University of California, Santa Cruz, California 95064
80 years ago, Joseph Larmor planted the seed that grew into today’s imposing body of knowledge about how the Earth’s magnetic ﬁeld is created. His simple idea, that the geomagnetic ﬁeld is the result of dynamo action in the Earth’s electrically conducting, ﬂuid core, encountered many difﬁculties, but these have by now been largely overcome, while alternative proposals have been found to be untenable. The development of the theory and its current status are reviewed below. The basic electrodynamics are summarized, but the main focus is on dynamical questions. A special study is made of the energy and entropy requirements of the dynamo and in particular of how efﬁcient it is, considered as a heat engine. Particular attention is paid to modeling core magnetohydrodynamics in a way that is tractable but nevertheless incorporates the dynamical effects of core turbulence in an approximate way. This theory has been tested by numerical integrations, some results from which are presented. The success of these simulations seems to be considerable, when measured against the known geomagnetic facts summarized here. Obstacles that still remain to be overcome are discussed, and some other future challenges are described.

CONTENTS
I. Groundwork A. Objectives B. A rudimentary model of the Earth C. What needs to be explained
II. Core convection A. Energy balance of the core B. Sources of energy C. Convection in a compressible ﬂuid; CO density D. Thermodynamic efﬁciency
III. Basic Dynamo Theory A. The induction equation B. Kinematic dynamos, Cowling’s theorem C. Turbulent helicity and the ␣ effect D. Large-scale helicity and ␣ effect
IV. Dynamical Theory A. The full (primitive) equations B. The reference state C. The anelastic equations D. Core turbulence E. Working equations and boundary conditions F. The Boussinesq approximation
V. RMHD A. Orders of magnitude B. Classical rotating ﬂows C. Magnetic effects D. The Taylor state and model z
VI. MHD Dynamo Simulations A. The development of models B. Some results
VII. The Future A. Turbulence, diffusion, and hyperdiffusion B. Boundary conditions C. The road ahead
Acknowledgments References

I. GROUNDWORK

A. Objectives 1081

1081

The Chinese were the ﬁrst to discover, probably sev-

1081 eral millennia ago, that the Earth is magnetic, but the

1083 reason for this remained a puzzle until quite recently.

1087 Theory ﬁrst started along the right lines early in the 20th

1087 century and has now reached the stage at which the fa-

1088 1090

vored mechanism of ﬁeld production has been simulated

1092 with considerable success. It is our pleasant task to try to

1093 explain the status of the subject as it stands today to

1093 physicists having no specialized knowledge of geomag-

1094 netism, and who may well believe that a subject based

1096 on 19th-century physics cannot possibly be interesting.

1098 In trying to convince them otherwise, we shall simplify

1098 the subject in two ways that may disturb the professional

1098 geophysicist. On the one hand, we shall ignore known

1099 geophysical facts that would complicate our arguments

1099 without adding anything essential to the physics. (The

1100 arguments are already complicated enough.) On the

1101 other hand, we shall try to avoid prevarication by mak-

1103 ing dogmatic assertions even when certainty is geophysi-

1103 cally unjustiﬁed.
1103

1104 B. A rudimentary model of the Earth 1107

1109

The Earth is an almost spherical body of mean radius

1110 1110 1111

rEϭ20 000/␲ km, consisting concentrically of a solid inner core, a ﬂuid outer core, a rocky mantle, and the thin

1116 crust on which we live. It is sketched in Fig. 1, where the

1116 following abbreviations, used throughout this review,

1118 appear:

1118

SICϭsolid inner core;

1120

1120

FOCϭﬂuid outer core;

Reviews of Modern Physics, Vol. 72, No. 4, October 2000 0034-6861/2000/72(4)/1081(43)/$23.60 ©2000 The American Physical Society 1081

1082

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

FIG. 1. Rudimentary sketch of the Earth’s internal structure: SIC, solid inner core; ICB, inner core boundary; FOC, ﬂuid outer core; CMB, core-mantle boundary.

ICBϭinner core boundary ͑ mean radius ¯rICB Ϸ1221.5 km);

CMBϭcore-mantle boundary ͑ mean radius rCMB Ϸ3480 km).

The Earth is in approximate hydrostatic equilibrium in an effective gravitational ﬁeld ¯g that is mostly its own (Newtonian) self-gravity but is also partly the centrifugal acceleration created by its rotation ⍀, which makes it slightly oblate. We have

“P¯ ϭ¯␳¯g,

(1.1)

¯gϭϪ“U¯ ,

(1.2)

where P is pressure, ␳ is density, and U is the effective gravitational potential. The overbar is used here to indicate horizontally averaged quantities (or more precisely quantities averaged over surfaces of constant U; see Sec. IV.B). Later, the overbar will also refer to the largescale component of the quantity under it. Table I gives some idea of how P¯ and ¯␳ increase with depth. (The sufﬁx M on the diffusivities listed there indicates that they are the ‘‘molecular’’ values, to distinguish them from turbulent values that will appear later.) There is a discontinuity ⌬␳ϭ600 kg mϪ3 in density at the inner core boundary.
The core has a density comparable to, but less than, what iron would have at core pressures, and it is probably an iron-rich alloy, made of an ‘‘uncertain mixture of all the elements’’ (Birch, 1952). The percentage by weight of iron is perhaps 85–90 %. It is not even known what the principal alloying element is, silicon, sulfur, oxygen, and hydrogen all having been proposed. To remove the ambiguity, it is supposed here that the core is made of FeX, where X stands for the principal, but unknown, light constituent. The last four entries of Table I are about right if XϭS or Si; see Braginsky and Roberts (1995), from which we draw most of our data. The inclusion of further light elements would greatly complicate the theory without adding signiﬁcant enlightenment.
To model seismic data successfully, it is necessary to

assume that the core is isentropic and chemically homogeneous:

“ ¯S ϭ 0,

(1.3)

“¯␰ ϭ 0.

(1.4)

Here S is the speciﬁc entropy and ␰ is the mass fraction of element X. From elementary thermodynamics and Eqs. (1.1)–(1.4) it follows that

“¯␳ ϭ¯␳¯g/ ¯u S2 ,

(1.5)

“T¯ ϭ¯␣S¯g,

(1.6)

“␮¯ ϭ¯␣␰¯g.

(1.7)

Here T is temperature, ␮ is chemical potential, uS is the speed of sound,

ͩ ͪ ͩ ͪ ␣

Sϭ

Ϫ

1 ␳

‫␳ץ‬ ‫ץ‬S

ϭ␳
P,␰

‫ץ‬T ‫ץ‬P

␣T

ϭ
S,␰

cP

(1.8)

is the entropic expansion coefﬁcient, and

ͩ ͪ ͩ ͪ ␣

␰

ϭ

Ϫ

1 ␳

‫␳ץ‬ ‫␰ץ‬

ϭ␳
P,S

‫␮ץ‬ ‫ץ‬P

S,␰

(1.9)

is the compositional expansion coefﬁcient. The coefﬁcient ␣S is more useful than the more familiar thermal expansion coefﬁcient ␣, to which it is related by the last term in Eq. (1.8), cP being the speciﬁc heat at constant pressure. The coefﬁcients ¯␣S and ¯␣␰ will occur frequently below; values are given in Table I.
In the general case, which is brieﬂy touched on in Sec.
IV, Eq. (1.1) shows that ¯␳ and P¯ and therefore all thermodynamic variables are constant over each equipoten-
tial surface, U¯ ϭconstant. We shall usually ignore centrifugal forces. The model Earth then has a spherically symmetric structure, and ¯␳, P¯ , T¯ , ␮¯ , ¯␣S, ¯␣␰ . . . are functions only of radial distance r from the geocenter O. (We shall denote ‫ץ‬/‫ץ‬r by ‫ץ‬r , and similarly for other differentials.) Since “•¯gϭϪ4␲G¯␳ where G is the gravitational constant, Eq. (1.2) gives

ٌ2U¯ ϭ4␲G¯␳.

(1.10)

In principle, a model of the Earth’s interior can be con-
structed from the seismically determined ¯uS by solving Eqs. (1.5) and (1.10) self-consistently for ¯␳ and U¯ . Then
P¯ is obtained by solving Eq. (1.1) subject to P¯ (rE)ϭ0. Similarly, T¯ and ␮¯ are obtained, apart from integration
constants, from Eqs. (1.6) and (1.7). For example, if we assume that ¯␣S and ¯␣␰ are constant, we obtain

T¯ ϪT¯ CMBϭϪ¯␣S͑ U¯ ϪU¯ CMB͒,

(1.11)

␮¯ Ϫ␮¯ ICBϭ¯␣␰͑ U¯ ICBϪU¯ ͒.

(1.12)

Equation (1.6) deﬁnes the adiabatic gradient Ϫ‫ץ‬rT¯ , which is less than 1 K kmϪ1; see Table I. The implied

outward heat ﬂux, mal conductivity),

i¯IsradaϭboϪutK0M.0‫ץ‬1r5T¯ W(wmhϪer2eatKtMheisintnheertchoerre-

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1083

TABLE I. Notation and magnitudes.

⍀ϭ7.292ϫ10Ϫ5 sϪ1

P¯ CMBϭ135.75 GPa

P¯ ICBϭ328.85 GPa

␳0ϭ10.9ϫ103 kg mϪ3

¯␳CMBϭ9.9ϫ103 kg mϪ3

¯␳IϩCBϭ12.166ϫ103 kg mϪ3

¯␳IϪCBϭ12.764ϫ103 kg mϪ3

⌬¯␳ ϭ (¯␳

ICB)

Ϫ ϩ

ϭ

0.6ϫ

103

kg

mϪ3

¯gCMBϭ10.68 m sϪ2

¯gICBϭ4.40 m sϪ2

U¯ CMBϪU¯ ICBϭ2.55ϫ107 m2 sϪ2 ¯uS,CMBϭ8.065ϫ103 m sϪ1 ¯uS,ICBϭ10.356ϫ103 m sϪ1

T¯ CMBϭ4000 K T¯ ICBϭ5300 K T¯ 0ϭ4590 K Ϫ(‫ץ‬rT¯ )ICBϭ0.28 K kmϪ1 Ϫ(‫ץ‬rT¯ )CMBϭ0.89 K kmϪ1 ¯␣CMBϭ1.8ϫ10Ϫ5 KϪ1 ¯␣ICBϭ1.0ϫ10Ϫ5 KϪ1 ¯cP,CMBϭ866 J kgϪ1 KϪ1 ¯cP,ICBϭ842 J kgϪ1 KϪ1 ¯␣CS MBϭ8.5ϫ10Ϫ5 kg JϪ1 K ¯␣ISCBϭ6.3ϫ10Ϫ5 kg JϪ1 K
¯␩ϭ2 m2 sϪ1
(K¯ M)CMBϭ40 W mϪ1 KϪ1 (K¯ M)ICBϭ53 W mϪ1 KϪ1 ¯␬Mϭ5ϫ10Ϫ6 m2 sϪ1 D¯ Mϭ10Ϫ9 m2 sϪ1 ¯␯Mϭ10Ϫ6 m2 sϪ1 ¯␣ ␰ ϭ 0.7 hϭ106 J kgϪ1
¯␰ ϭ 0.16
¯␰ SICϭ 0.1 ⌬T¯ mϭ700 K

Well-determined parameters

Angular velocity of the Earth Pressure at the CMB

Pressure at the ICB

Mean density of the FOC Density of the FOC at the CMB Density of the FOC at the ICB Density of the SIC at the ICB Density jump at the ICB Acceleration due to gravity at the CMB Acceleration due to gravity at the ICB Gravitational potential difference across FOC

Thermodynamic properties

Seismic velocity in FOC at the CMB Seismic velocity in FOC at the ICB
Temperature of the CMB

Temperature of the ICB

Mean temperature of the FOC

Adiabatic gradient at ICB

Adiabatic gradient at CMB

Thermal expansion coefﬁcient at the CMB Thermal expansion coefﬁcient at the ICB Speciﬁc heat at constant P at the CMB Speciﬁc heat at constant P at the ICB Entropic expansion coefﬁcient at CMB Entropic expansion coefﬁcient at ICB Relevant but less well-determined quantities Magnetic diffusivity Thermal conductivity at CMB

Thermal conductivity at ICB

Thermal diffusivity of FOC Compositional diffusivity of FOC

Kinematic viscosity of FOC Compositional expansion coefﬁcient Latent heat of crystallization Mass fraction of light constituents in FOC

Mass fraction of light constituents in SIC

Depression of melting point through alloying

boundary and 0.036 W mϪ2 at the core-mantle boundary. Multiplying by the area of these surfaces, we ﬁnd that the adiabatic heat ﬂow is

HICBϭ0.3 TW,

(1.13)

HCMBϭ5.4 TW.

(1.14)

Being highly metallic, the electrical conductivity ␴M of the core is large. Taking ␴Mϭ4ϫ105 S mϪ1 in the ﬂuid outer core, we see that the more useful magnetic diffusivity ¯␩ϭ1/␮0␴M is 2 m2 sϪ1, a value we adopt for

the solid inner core also. Because T¯ is so large, we take the magnetic permeability to be that of free space, ␮0 ϭ4␲ϫ10Ϫ7 H mϪ1.

C. What needs to be explained

Sufﬁciently far from the Earth, its external magnetic ﬁeld Bˆ is dipolar. In mathematical terms,

Bˆ ϭϪ“V,

(1.15)

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1084

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

“•Bˆ ϭ0,

(1.16)

where

VϳϪm•“rϪ1, for r→ϱ.

(1.17)

Here m(t) is the strength of the ‘‘centered dipole’’ at time t; currently mϷ7.835ϫ1022 A m2. Its direction is in-

clined to the geographic axis Oz by only about 11° at the

present time, or perhaps we should say ‘‘169°,’’ since the

dipole axis is almost antiparallel to the Earth’s angular velocity ⍀, the South magnetic pole being close to the

North geographic pole. It is this proximity of the poles

that makes the magnetic compass needle such a boon to

seafarers. The fact that Bˆ is represented by a ﬁeld that

diminishes with distance from the Earth lends support to

the conjecture made by William Gilbert in 1600 that the

origin of the Earth’s magnetism lies within it.

The dipole (1.17) is merely the ﬁrst term in an expan-

sion of V that includes quadrupoles, octupoles, etc. and

is more economically represented by an expansion in

spherical

harmonics,

P

m l

(

␪

)

cos

m␾

and

P

m l

(

␪

)

sin

m␾,

where

the

P

m l

are

Legendre

functions,

␪

is

colatitude,

and ␾ is longitude. (Colatitude differs from latitude by

1 2

␲

,

the

North

pole

being

␪ϭ0

and

the

South

pole

␪

ϭ␲.) The full expansion of V(r,␪,␾,t) is

͚ ͚ ͩ ͪ ϱ l
VϭrE
lϭ0 mϭ0

rE r

lϩ1

P

m l

͑

␪

͒

͓

g

lm͑

t

͒

cos

m

␾

ϩ

h

m l

͑

t

͒

sin

m

␾

͔

.

(1.18)

The

coefﬁcients

g

m l

and

h

m l

are

named

after

Gauss,

who

in 1839 was the ﬁrst to analyze the geomagnetic data in

this way. In addition to the interior harmonics proportional to (rE /r)lϩ1 appearing in Eq. (1.18), Gauss also
included the exterior harmonics, which are proportional to (r/rE)l. He showed that these are so weak as to be
effectively nonexistent. Equations (1.15) and (1.18) are

tenable only because the electric current density at and
above rϭrE is essentially zero. Gauss’s analysis supported Gilbert’s conjecture, but

also suggested more. Let the mean magnetic energy den-

sity on the spherical surface of radius r be W(r,t). For
rϾrE , this may be expressed as the sum of contributions from each spherical harmonic:

͚ W͑r,t ͒ϭ Wl͑r,t ͒. l

(1.19)

The set of values Wl is called the power spectrum of the ﬁeld. Two approximate relations are found for the Earth’s present surface ﬁeld (Langel and Estes, 1982):
ͭ Ϫ3.270Ϫ0.569l, for 2рlр12,
log10 Wl͑ rE͒Ϸ Ϫ10.83Ϫ0.0114l, for 16рlр23. (1.20)

For 13рlр15 the spectrum makes a transition between these two linear relations. The largeness of the constant 0.569 in the ﬁrst relation indicates (but does not prove) that the sources of these terms lie far below the Earth’s surface. They create the main geomagnetic ﬁeld, or sim-

ply (as we shall call it) the ‘‘main ﬁeld.’’ This is the object of our study. The smallness of the constant 0.0114 in the second relation indicates (but does not prove) that the sources of these terms lie close to the Earth’s surface; they are attributed to the permanent magnetism of the crust. The crustal sources in the sum (1.18) swamp the main ﬁeld sources for lϾ12, thus destroying all the short-wavelength information about the main ﬁeld; for our purposes, the sum can be truncated at lϭ12. This sets a limit of about 5 nT on the accuracy to which we can know the main ﬁeld at the Earth’s surface.1
Since the lth harmonic of Bˆ is proportional to (rE /r)lϩ2, we have

Wl͑r͒ϭWl͑rE͒͑rE /r͒2lϩ4,

(1.21)

so that, if applied at the core surface, the ﬁrst relation (1.20) would give

log10 Wl͑ rCMB͒ϷϪ2.221Ϫ0.044l, for 2рlр12. (1.22)
The fact that the coefﬁcient of l is negative means that, if taken to lϭϱ (the crustal sources having been removed), the series (1.20) would converge at the core surface. This suggests that there are no further sources of magnetism in the mantle and in particular no electric currents ﬂowing there. The mantle is certainly a poor conductor compared with the core, and we shall usually assume it is insulating. Obviously, the expansion (1.20) has no meaning below the core surface where Eq. (1.15) is drastically violated, but the fact that the constant 0.044 in Eq. (1.22) is so small is unexpected and suggests that there is abundant small-scale structure on the core surface.2
Figure 2 shows typical magnetic energy spectra at the core surface ͓Wl(rCMB)͔ through degree 12 for the Earth in 1980 according to Langel and Estes (1985) and through degree 95 for a snapshot from our simulation (Sec. VI below). In both cases, the energy in the dipole (lϭ1) exceeds that in any other harmonic. We shall call this dipole dominance. Since the dipole part of the ﬁeld decreases the least rapidly with distance from the core, the dipole dominates even more strongly at the Earth’s surface. This is illustrated in Fig. 3, where the radial component of the ﬁeld is plotted in equal area projections. The three panels on the left are projections on the Earth’s surface; the three on the right are for the coremantle boundary. The top panels depict the real geomagnetic ﬁeld truncated at lϭ12; the center panels show

1Since crustal sources, if they change at all, do so on time scales much longer than the core’s, one could in principle derive information about the time rate of change of the Gauss coefﬁcients of the main ﬁeld for lϾ12.
2The fact that the power spectrum of the ﬁeld is nearly ﬂat at the core surface led Hide (1978) to propose that the core radii of other planets and satellites in the solar system that generate their own magnetic ﬁelds could be determined by analyzing their ﬁelds (as determined in ﬂybys) into spherical harmonics and then extrapolating these to the depth at which the power spectrum becomes ﬂat. See also Hide and Malin (1979, 1981).

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1085

FIG. 2. Magnetic energy density Wl(rCMB) at the core-mantle boundary as a function of spherical harmonic degree l for
(solid symbols) the Earth in 1980 and (open symbols) a snap-
shot from the Glatzmaier-Roberts simulation. Values of Wl have been multiplied by 2␮0 .

our simulation, also truncated at lϭ12, while the bottom panels show the same ﬁelds at the maximum truncation level of the simulation, lϭ95. One can see that, although including degrees 13–95 produces no detectable differ-
ence in the surface ﬁeld, it makes a signiﬁcant difference in the structure of the ﬁeld at the core-mantle boundary. The intense concentrations of magnetic ﬂux (core spots) seen in the simulated ﬁeld at this boundary, which are essentially undetectable from the (ﬁltered) geomagneticﬁeld measurements at the Earth’s surface, are somewhat similar to sun spots on the solar surface. These ﬁgures graphically illustrate the difﬁculty in inferring the ﬁeld at the core-mantle boundary from the observed ﬁeld at the Earth’s surface.
Although the power spectrum at the core surface is unknown for the real Earth beyond lϭ12, it cannot be absolutely white, as that would imply an inﬁnite W and an inﬁnite magnetic energy EˆB outside the core. If we adopt law (1.22) for all lу2, and add in the contribution from the dipole (which is more than 70% of the whole), we ﬁnd that EˆBϷ6.9ϫ1018 J. If there were no ﬁeld sources, this energy would be carried by the Poynting ﬂux back into the core and be dissipated there through ohmic heating in a time of order (rCMB /␲)2/¯␩Ϸ2 ϫ104 yr. This is the longest efolding time for the decay modes of ﬁeld in a stationary sphere of radius rCMB ; it belongs to the lϭ1 harmonic. We may write it as

␶␩ϭL¯ 2/¯␩,

(1.23)

with L¯ ϭrCMB /␲Ϸ103 km. The efolding times of all other harmonics are shorter and, to make an allowance

for this, we shall take as our estimate of the characteristic length scale of ﬁelds in the core

L¯ ϭ500 km.

(1.24)

The overbar indicates that this refers to the large scales of core magnetohydrodynamics (MHD). The electromagnetic time constant (1.23) from the estimate (1.24) is 4ϫ103 yr and corresponds roughly to lϭ3. To maintain the external magnetic energy Eˆ, the power delivered by the sources must exceed EˆB/␶␩Ϸ45 MW. This considerably underestimates the total energy requirements; we shall ﬁnd in Sec. II.A that the magnetic ﬁeld in the core (which we shall denote by B to distinguish it from the ﬁeld Bˆ in and above the mantle) makes larger energy demands than 45 MW.
In reality the main ﬁeld is not constant, but varies slowly in time, t, a phenomenon called the secular variation. It has many time scales. These range3 from secular variation impulses or ‘‘jerks’’ (Le Moue¨ l and Courtillot, 1981; Le Moue¨ l et al., 1982) that are completed in about one year to total reversals of the polarity of the ﬁeld, the last of which occurred about 7.8ϫ105 yr ago. Knowledge of the ﬁeld over such long times is derived from archeomagnetic and paleomagnetic data. Paleomagnetism is the study of the main ﬁeld from Precambrian times to the present, as revealed by the ﬁeld recorded by rocks and sediments at their birth. Archeomagnetism provides similar information from historic and prehistoric times from the ﬁeld imprinted on manmade artifacts such as shards of pottery and the bricks from the kilns that made it. We shall not describe here the processes by which ancient ﬁelds are trapped, nor the techniques and pitfalls of extracting reliable geomagnetic information from samples. We simply report conclusions that are relevant to our theme. More detailed information can be found in Merrill et al. (1996).
Perhaps the most striking facts that have emerged from paleomagnetism are that the Earth has possessed a ﬁeld for more than 3ϫ109 yr and that, except during a polarity reversal (a process that is typically completed in a few thousand years), its intensity has not varied by a factor of more than about 3 over most of geological time (Kono and Tanaka, 1995). The age of the ﬁeld therefore greatly exceeds ␶␩ , and it is necessary to ﬁnd an energy source for it. The ﬁeld and currents cannot be relics of the Earth’s creation.
We have little to say about secular variation impulses. They are small but have not yet been convincingly explained. They seem to be worldwide phenomena and to have their origin in the core. Two observations suggest this. First, the jerks seem to be correlated to changes in the westward drift rate (Le Moue¨ l et al., 1981; Gire et al.,

3There are probably also much shorter time scales, but these are shielded from observation by electrical conduction in the mantle, which has an electromagnetic time constant of order 1 yr. See, for example, Gubbins and Roberts (1987).

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1086

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

FIG. 3. Radial component of the magnetic ﬁeld (reds for outward directed and blues for inward) plotted at the Earth’s surface and at the core-mantle boundary. The surface ﬁelds are multiplied by 10 to obtain comparable color contrast. The Earth’s ﬁeld is plotted out to spherical harmonic degree 12. A snapshot from the Glatzmaier-Roberts simulation is plotted out to degree 12 (for comparison) and out to degree 95 [Color].

1983). Second, they also seem to be correlated with equally rapid variations in the length of the day. It has long been known that the Earth is not a perfect time keeper and that the length of the day can change by as much as 5 ms in 10 yr during so-called decade variations. The magnitude of these changes is too large to be explained by a transfer of angular momentum between the mantle and the atmosphere. The global wind system would have to more than reverse or more than double to explain it. The rapidity of the changes shows that the ﬂuctuating unbalanced torque acting on the mantle is

typically of order4 1018 N m, and this also is too great to be provided by the atmosphere or oceans. The exchange of angular momentum must be between mantle and
4The axial moment of inertia of the mantle is 7.12 ϫ1037 kg m2, so that a ﬂuctuation of 1 ms, corresponding to a change in angular velocity ⍀ of the mantle of ␦⍀Ϸ8.4 ϫ10Ϫ13 sϪ1, is associated with a change in angular momentum of 5.07ϫ1025 kg m2 sϪ1. The torque required to bring this about in a year is 1.6ϫ1018 N m.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1087

core. Correlations between westward drift and mantle rotation have long been suspected (see, for example Le Moue¨ l et al., 1981), and there are indications (Jackson et al., 1993; Jault et al., 1988) that changes in the angular momentum of the mantle are accompanied by equal and opposite changes in the angular momentum of the core associated with a ‘‘torsional wave’’; see Sec. V.C below.
The secular variation operates on many time scales intermediate between the two extremes described above. Because of the historical importance of the geomagnetic ﬁeld to mariners, the ﬁeld has been increasingly well studied since the ﬁfteenth century, and many maps have been constructed that describe its structure at the Earth’s surface and its variation over the last 400 years (see Langel, 1987). One striking phenomenon is the westward drift of the ﬁeld patterns, a phenomenon discovered by Halley, who presciently attributed it to ﬂuid motions inside the Earth. The drift is latitude dependent and irregular. Its typical value in midlatitudes is commonly used to estimate U¯ , the characteristic largescale velocity in the core:

U¯ ϭ5ϫ10Ϫ4 m sϪ1.

(1.25)

Three other paleomagnetic facts will be of interest to us later. First, consider Fig. 4, which shows the location of the South geomagnetic pole in the northern hemisphere at various times during the past 5 Myr. (More precisely, Fig. 4 shows the virtual geomagnetic pole, or VGP, which is computed at a site from a measurement of the direction of the prevailing magnetic ﬁeld. This pole represents the position that the south magnetic pole would have were the ﬁeld precisely dipolar at the site.) Figure 4 suggests that the present angle of about 10° between the geographic and geomagnetic axes is not uncharacteristic of the past, but archeomagnetic data (not shown here) from the past 1000 years indicates that the average angle is roughly 5°. When averaged over a time of order 104 yr, the geographic and geomagnetic axes coincide.5 Several points shown in Fig. 4 lie far from the North geographic pole, and very often in the past the South geomagnetic pole has actually moved to the southern hemisphere and remained there until later another ﬂuctuation has returned it to the northern hemisphere. These are geomagnetic ﬁeld reversals and they bring us to our second interesting fact: there is no ﬁrm evidence that the Earth prefers to be in one polarity state rather than the other. If there is a bias, it is too small to be convincing. Our third interesting fact is the irregularity of reversals. Figure 5 shows the reversal frequency during the past 165 Myr. It may be seen that

5Over time scales longer than 5 Myr, the mean geomagnetic poles appear to move over the Earth’s surface in ‘‘polar wander’’ paths. This provides the main tool used to infer how continents have moved relative to one another over geological time. Further discussion of this fascinating topic, and the related subject of sea ﬂoor spreading, would take us too far from our present objectives; see, for example, McElhinny (1973) and Merrill et al. (1996).

during a period of about 35 Myr during the Cretaceous (144-66 Myr ago) the ﬁeld remained of one polarity state. It is also known that, during 50 Myr in the PermoCarboniferous, the ﬁeld did not reverse. At the present time, however, the ﬁeld is reversing roughly 4–5 times per Myr. When a constant-polarity interval endures for more than 105 yr, it deﬁnes a polarity epoch, otherwise it is called a polarity event. When a magnetic pole moves further from the nearest geographic pole than 45°, but then returns, a polarity excursion has occurred. There is evidence that 14 of these have so far occurred in the present polarity epoch (Lund et al., 1998). On occasions, called ‘‘cryptochrons,’’ a pair of back-to-back reversals occur in quick succession, restoring the initial polarity state. Sometimes the magnetic ﬁeld may appear to reverse at one site where the data was collected but not at another, suggesting that the nondipole ﬁeld dominates the dipole ﬁeld during this time.
The geomagnetic facts just summarized raise many challenging questions, and we aim to answer in at least a qualitative way virtually all of them.
II. CORE CONVECTION
A. Energy balance of the core
The temperature T beneath the crust exceeds the Curie point of all known materials; permanent magnetism does not exist anywhere in the Earth, except the crust. The obvious explanation of the Earth’s magnetism is therefore untenable; the main ﬁeld is created by electric currents ﬂowing mainly in the core.6 Several possible origins of these currents have been proposed,7 but all except one have been found wanting. The favored idea today is that they are generated by self-excited dynamo action associated with the motion of core ﬂuid, a suggestion ﬁrst made by Larmor (1919). We discuss the foundations of dynamo theory in Sec. III.
What drives core motions? It is obvious that, by Lenz’s law, the Lorentz force created by currents induced by core motion oppose that motion and bring it to rest, unless an energy source is available to drive the motion. Looked at another way, the ﬂuid motions must supply the Joule heat losses of the electric currents. If they fail to do so, the ﬁeld will diminish until dynamo action ceases. How large are these losses? Let us estimate the Joule losses QJ from B¯ as QJϭ(J¯ 2/
6We ignore the by now discredited idea that all rotating bodies produce a ﬁeld because of their rotation (Blackett, 1947).
7For example, some geophysicists have argued that the currents are driven by thermoelectric potential differences between rising and falling convection currents or between the core and the mantle. According to others, the currents have an electrochemical origin, as in a battery. The NernstEttinghauser effect has also been suggested. Like the Blackett proposal, all these theories ﬁnd it impossible to explain polarity reversals and the apparent indifference of the ﬁeld to its polarity state, i.e., the sign of Bˆ .

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1088

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

FIG. 5. Estimated reversal rate of the geomagnetic ﬁeld over the past 160 Myr. From Merrill et al., 1996.

FIG. 4. Paleomagnetic poles north of 40° N during the past 5 Myr in polar stereographic projection; the white diagonal cross is the present South magnetic pole. From McElhinny, 1973.

␴M)Vcore , where VcoreϷ1.77ϫ1024 m3 is the volume of the core. Using Ampe` re’s law [see Eq. (3.4) below], we estimate the characteristic current density J¯ for a typical ﬁeld strength of B¯ as J¯ ϷB¯ /␮0L¯ . In Sec. V.A we shall assess B¯ by equating the Coriolis force 2⍀¯␳VϷ2⍀␳0U¯ , where V is ﬂuid velocity, to the Lorentz force JÃBϷJ¯ B¯ , so that B¯ Ϸͱ(2⍀␳0␮0U¯ L¯ )Ϸ20 mT, J¯ Ϸ0.04 A mϪ2, and QJϷ1 TW. These values are prob-
ably overestimates. In the simulations described in Sec.
VI, the maximum ﬁeld strength is of order 20 mT and
the typical ﬁeld strength is closer to 5 mT, and it is found
that

QJϷ0.3 TW.

(2.1)

It is hard to be conﬁdent about estimate (2.1), which is clearly sensitive to the assumed B¯ and L¯ . It is even harder to assess the energy dissipation rate Qt of core turbulence, which may be at least as large as QJ, though the energy dissipation rate Q␯ from the large-scale ﬂow is likely to be much smaller; see Secs. IV and V. The turbulent energy loss is mainly through Joule dissipation of the small-scale currents, and only slightly from the viscous dissipation of the small-scale motions; see Sec. IV.D. The problem of estimating the total dissipation rate,

QD ϭ QJ ϩ Q␯ ϩ Qt ,

(2.2)

will be the principal concern of Sec. II.C. The source of QJ has been argued about for most of
the last 50 years, and by now the main consensus is that the ﬂow is convectively driven. The only alternative explanation that has not been ruled out is that the motions are powered by the luni-solar precession of the Earth’s rotation axis; see Bullard (1949), Malkus (1963, 1968), and Vanyo (1991). In this review we go with the major-

ity and assume that core motions are driven by buoyancy. The rate of working QC of the buoyancy forces
must balance the dissipative losses:

QC ϭ QD ϵ QJ ϩ Q␯ ϩ Qt .

(2.3)

This statement, like Eqs. (2.9) and (2.33) below, should

be interpreted as an average, since the quantities in-

volved ﬂuctuate over time; see Sec. 7 of Braginsky and

Roberts (1995) for a detailed derivation. An expression for QC is given in Eq. (2.15) below.

If we use B¯ ϭ5 mT and the estimate (1.25) for U¯ to

calculate the magnetic energy density as B¯ 2/2␮0 and the kinetic energy density (of core motions relative to the

mantle)

as

1 2

␳

0

U¯

2

,

we

obtain

10

J mϪ3

and

1.4

mJ mϪ3,

respectively. The ratio of the magnetic energy EBϷ1.8

ϫ1021 J of the Earth to the kinetic energy of core mo-

tions EVϷ2.4ϫ1017 J is therefore large: EB/EVϷ7500; see

Sec. V.

B. Sources of energy
What is the origin of the buoyancy? Most geochemists argue that there is no signiﬁcant radioactivity in the core to heat it; see, for example, Stacey (1992).8 It is now thought that the buoyancy arises from processes driven by the slow cooling of the Earth over geological time. The Earth continually radiates heat into space. This heat comes partly from radioactive sources in the crust and mantle, but the Earth loses more heat than that, i.e., the temperature of the Earth is dropping. Within the Earth, the heat is carried outwards mainly by thermal conduction down the adiabatic temperature gradient (1.6). But this state is convectionally unstable; see Sec. II.C. Both mantle and (probably) core also carry heat outward by convection, though, because their viscosities are so very different, mantle motions are measured in cm/yr,

8For a recent reassessment of the radiogenic heat release in the core, see Chabot and Drake (1999).

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1089

FIG. 6. Phase diagram (left) for the binary alloy representing the core and its relationship (right) to core structure. The arrows from the solid inner core (SIC) represent the outward ﬂuxes of heat and light component of the alloy from the inner core boundary as it freezes; the arrows in the mantle represent the heat ﬂux from core to mantle, which acts as a ‘‘valve’’ controlling core convection.
whereas those in the ﬂuid core are about six orders of magnitude greater, i.e., a few mm/s. Correspondingly, the overturning time ␶c of core convection is measured in hundreds of years and that of the mantle ␶m in hundreds of millions of years.
The downward temperature gradient created in a ﬂuid layer cooled sufﬁciently strongly from the top may cause it to lose convective stability, but this is not an efﬁcient way to make it convect. For given buoyancy sources, stronger motions arise when the sources are deep in the layer rather than distributed volumetrically throughout it. Best of all, they should be at the base of the layer. The core achieves this in an interesting way (Braginsky, 1963). For simplicity, let us (as is commonly done; see Sec. I.B) model the core as a binary, iron-rich alloy. A phase diagram usually displays the solidus and liquidus of the alloy as curves in ␰T-space, where ␰ is the mass fraction of the (unknown) light constituent. This is because freezing is usually discussed in contexts where variations in pressure P are unimportant. In reality, the liquidus and solidus depend on the thermodynamic state of the material. They should therefore appear as surfaces in, say, ␰PT-space. The traditional diagrams are merely projections of these surfaces onto the appropriate constant-P plane. For the core it makes more sense to plot the solidus and liquidus in ␰PS-space and to project them onto the appropriate constant-S plane. The result is a ␰P plot, such as that shown in Fig. 6.
On descending through the ﬂuid core from the coremantle boundary, we eventually encounter the inner core boundary (ICB), where solid freezes from the overlying melt. In fact, the core is a body of ﬂuid that is cooled from the top but freezes from the bottom! This may seem paradoxical, but only because we are not familiar with such processes in everyday life on Earth’s surface. The paradox is inevitable whenever the melting point of a material increases more rapidly with increasing pressure than the ambient temperature does. Although the ICB is a freezing front, the material lying below it is unlikely to be completely solid. Rather, it will be a mixed-phase region, where liquid and solid coexist

(Loper and Roberts, 1981); this is often called a ‘‘mush’’

in the metallurgical literature. This idea is corroborated

by seismically determined estimates of the Q of the solid

inner core (SIC), which are fairly low; Vidale and Earle

(2000) give values between 244 and 450. The mush be-

neath the ICB is called the ‘‘solid’’ inner core, but the

adjective is appropriate because the mass fraction ␸ of

solid is likely to increase rapidly with depth and become

close to 1 within a few meters of the ICB. There are two

reasons for this. First, the ﬂuid that percolates through

the channels of the mush deposits solid on them and

gradually closes them up. Second, a mush is of low me-

chanical strength and will easily compact in the prevail-

ing pressure gradient. The upper layers of the SIC, in

which ␸ differs signiﬁcantly from 1, impart a ‘‘fuzziness’’

to the ICB that we shall disregard in the remainder of

this review. The mechanical weakness of the SIC sug-

gests that it cannot sustain signiﬁcant internal nonhydro-

static stresses. It is therefore close to being in hydro-

static equilibrium and the ICB is close to being an

equipotential surface, a fact used in Sec. IV.

The idea that the SIC is only a solidiﬁed form of the

ﬂuid lying above it, and that the inner core surface is a

freezing front that advances slowly into the ﬂuid as the

Earth cools, was ﬁrst suggested by Jacobs (1953) and is

now generally accepted. If the SIC has added to its mass MSICϷ9.7ϫ1022 kg at a uniform rate over most of its
age, which (for consistency with other estimates made

below) we take to be ¯␶ϭ1.2 Gyr, the rate of increase in

its mass is 2.5ϫ106 kg sϪ1, corresponding to a rate of

amdmvaynrcϪe1.oVf ethrheoIoCgBencu(1rr9e6n1t)lypooifnrថteICdBϷou1t0thpamt

sϪ1 the

or 0.3 latent

heat released at the inner surface as it freezes provides a

thermal buoyancy source to stir the core. The heat re-

lease at the ICB is

QLϭT¯ ICB⌬¯SM˙ SIC ,

(2.4)

where ⌬¯Sϭ¯SϩϪ¯SϪ is the jump in ¯S at the SIC. This is
closely related to the latent heat of crystallization h,
which is usually deﬁned as the discontinuity in enthalpy: hϭT⌬Sϩ␮⌬␰. Since h is rather uncertain, there is little
point in distinguishing between T¯ ICB⌬¯S and h. From our estimate of rថ ICB above and Table I, we ﬁnd that

QLϷ2.5 TW.

(2.5)

Braginsky (1963) remarked that the difference in the

compositions ␰ of the liquidus and solidus (see Fig. 6)

means that the inner core is richer in iron than the outer

core, and that this largely accounts for the observed density jump ⌬¯␳ϭ¯␳ϪϪ¯␳ϩ at the ICB, i.e., the density jump

is due less to the contraction of a material as it freezes

(an effect we do not include) than to a discontinuity

⌬¯␰

ϭ¯␰

Ϫ¯␰

Ϫ ICB

in

¯␰ .

This

also

implies,

he

noted,

that

some

of the light component of the alloy is released as the

inner core freezes and that this provides a second, com-

positional buoyancy source to stir the core. Suppose that

the core ﬂuid is 16% by weight of light component and

that it retains 40% of this when it freezes onto the sur-

face of the SIC. The remainder provides an average

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1090

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

mass ﬂux Ir␰ϭ¯␳IϪCB⌬¯␰rថ ICB of light component away from the ICB of approximately 8.7ϫ10Ϫ9 kg mϪ2 sϪ1 and a

rate of increase of ¯␰ in the ﬂuid outer core (FOC) of ¯␰

Ϸ 8.9ϫ 10Ϫ 20

sϪ1.

(In

reality,

rថ ICB ,

¯I

␰ ICB

,

and

¯I

S ICB

vary

over the ICB but, since each depends on the rate of

freezing, they are everywhere proportional to one an-

other; see Sec. IV.E.) If ␰ were ever to approach the

eutectic, compositional forcing would cease.

The fact that the solid core is growing and is richer in

iron than the ﬂuid core means that the Earth is continu-

ally becoming more centrally condensed as it cools. In

addition, therefore, to the release of internal energy, of

which the latent heat is the most signiﬁcant part, there is a release of gravitational energy QG, which we may es-

timate by the following argument. A change ␦¯␰ in ¯␰ cre-

ates a fractional change in outer core density of

Ϫ¯␳¯␣␰␦¯␰, which is the same everywhere in the FOC if (as we shall temporarily assume for simplicity) ¯␣␰ is con-

stant. We compute the change in gravitational energy by

imagining that the mass (¯␣␰␦¯␰)¯␳d3x in a volume element d3x is moved directly to the ICB, to form part of

the mass 4␲¯rI2CB␦¯rICB⌬␳ϭMFOC¯␣␰␦¯␰ added to the SIC through freezing. The gravitational energy released is

¯␣␰␦¯␰(U¯ ϪU¯ ICB)¯␳d3x. Summing over the FOC and supposing that the increase ␦¯␰ in ¯␰ happens in time ␦t, we

obtain in the limit ␦¯␰/␦t→␰ថ

͵ QG ϭ Ϫ ¯␣ ␰ ␰ថ

͑ U¯ ICBϪU¯ ͒¯␳d3x.

FOC

(2.6)

Using ␰ថ Ϸ8.9ϫ10Ϫ20 sϪ1, we ﬁnd that

QGϷ1.3 TW.

(2.7)

Equation (1.12) allows us to express QG differently:

͵ QGϭϪ␰ថ

͑ ␮¯ ICBϪ␮¯ ͒¯␳d3x,

FOC

(2.8)

and this is the appropriate form for QG if the assumption of constant ¯␣␰ is lifted.
The average heat balance of the FOC is expressed by9

QCMBϭ HICBϩ QL ϩ QG ϩ QS ϩ QR ,

(2.9)

where HICB has been estimated in Eq. (1.13), the subscript CMB refers to the core-mantle boundary, QR is

the rate of energy input from core radioactivity (if any), and QS is the heat loss through the cooling of the FOC:

͵ QSϭϪSថ

T¯ ¯␳d3x.

FOC

(2.10)

9The form (2.9) differs from Eq. (7.32) of Braginsky and Roberts (1995), who give QCMBϭQLϩQGϩQtSotalϩQtRotal , where QtSotal and QtRotal also include the solid inner core. The form (2.9) excludes the SIC by supposing that the heat ﬂux from SIC
to FOC arises only from the adiabatic gradient in the SIC,
which is assumed to be continuous across the ICB. Again, Eq.
(2.9) should be interpreted as an average over the convective
time scale.

Equation (2.10) may also be written as QS

ϭϪT¯ 0Sថ MFOC where T¯ 0 is the mass-weighted average of T¯ in the FOC; see Table I. We give later, in Eq.

(4.66), a means of computing S˙ (which is proportional to

rថ ICB), and this gives S˙ ϷϪ3.7ϫ10Ϫ16 W kgϪ1 KϪ1. Since MFOCϭ1.8367ϫ1024 kg, we now have

QSϷ3.1 TW.

(2.11)

Assuming that QRϭ0, we ﬁnd from Eq. (2.9) that QCMBϭ7.2 TW.
Since QL, QG, and QS are proportional to rថ ICB , we can easily explore other possibilities, such as the follow-
ing: (a) if QCMBϭHCMB and QRϭ0, then ¯␶ϭ1.6 Gyr; (b) if ¯␶ϭ4 Gyr and QRϭ0, then QCMBϭ2.4 TW ϽHCMB , so that convection pumps heat downward; see Sec. II.C; (c) if QRϾQCMBϪHICB , then the SIC would start to remelt (rថ ICBϽ0).

C. Convection in a compressible ﬂuid; CO density

Convection experiments in the laboratory usually in-

volve thin layers of ﬂuid that are almost uniform in den-

sity, but the compressibility of a thick shell of ﬂuid like

the ﬂuid outer core cannot be ignored, since (see Sec.

I.B) the adiabatic increase of T¯ with depth results in a heat loss HCMB from the core far exceeding the energy demands QD of the dynamo.

Let us consider ﬁrst the case of a chemically homoge-

neous ﬂuid in which ¯␰ is constant so that compositional

buoyancy is absent. Suppose that the ﬂuid is in a hydro-

static equilibrium in which the speciﬁc entropy ¯S de-

pends on depth, and apply a ‘‘test parcel’’ argument (see

also Sec. 4 of Landau and Lifshitz, 1987). Imagine that a

small parcel of ﬂuid is inﬁnitesimally displaced down-

wards so fast that heat has no time to change its entropy

content but so slowly that its internal pressure adjusts to

that of its new surroundings. Its density will therefore

increase

by

␦

␳

ϭ

␦

P¯

/¯u

2 S

.

The

density

of

the

new

sur-

roundings is greater than that of the old surroundings by

␦¯␳

ϭ

␦

P¯

/

¯u

2 S

Ϫ¯␳

¯␣

S

␦

¯S

,

where

¯␣ S Ͼ 0.

If

␦ ␳ Ͼ ␦¯␳ ,

the

parcel

is denser than its new surroundings and will tend to sink

further. The equilibrium is then said to be convectively

unstable (although, strictly speaking, the argument ig-

nores thermal conduction and viscosity, both of which

tend to stabilize the equilibrium; see Sec. V). If ␦␳

Ͻ␦¯␳, the buoyancy force will drive it upwards towards

its original position; this deﬁnes convective stability; ␦␳

ϭ␦¯␳ signiﬁes neutral stability. Stated another way, the

ﬂuid is unstable if ¯S increases downwards anywhere, and

it is stable if ¯S decreases downwards everywhere (as, for

example, happens if the ﬂuid is isothermal). The down-

ward temperature gradient in the isentropic state of uni-

form ¯S and neutral stability is the adiabatic (tempera-

ture) gradient deﬁned in Sec. I.B. When the heat ﬂow QCMB from the top of the layer
exceeds HCMB given by Eq. (1.14), the layer convects

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1091

and the difference QCcoMnvB is the convective heat ﬂowing from core to mantle. To estimate QCcoMnvB , we suppose that the convective velocities Uϳ5ϫ10Ϫ4 m sϪ1 are of

the order inferred in Sec. I.C and that the temperature

difference between the rising and falling convective streams is of order only ␦Tϳ10Ϫ4 K. (We shall ﬁnd in

Secs. V.B and VI.B that this estimate of ␦T is realistic

for

the

core.)

The

outward

convective

heat

ﬂux

I

conv r

ϳ ␳ 0¯c P U␦ T

then

signiﬁcantly

exceeds

¯I

ad r

,

and

the

ﬂuid

may be said to be convecting strongly. Vigorous convec-

tion homogenizes all extensive properties of the ﬂuid, so that ¯S becomes almost uniform everywhere, except in

boundary layers. As we saw in Sec. I.B, the success of

models of the Earth’s interior, based on seismic and

other data, depends on the assumption that the core is

close to being in an isentropic state, and this success is a

strong indication that it is convecting vigorously.

Let us now ask if and how compositional buoyancy

changes any of these conclusions. Suppose that ¯S and ¯␰ depend on depth and that the parcel retains both its ¯S and ¯␰ contents when it sinks. The increase in its density

is

again

␦

␳

ϭ

␦

P¯

/¯u

2 S

,

but

that

of

the

surroundings

is

now

␦¯␳ϭ␦P¯ /¯uS2Ϫ¯␳(¯␣S␦¯Sϩ¯␣␰␦¯␰). Again the sign of ␦␳Ϫ␦¯␳

is all important. Stability means that Ϫ¯g•(¯␣S“¯S

ϩ¯␣␰“¯␰)Ͼ0 everywhere, but if

Ϫ¯g•͑ ¯␣S“¯Sϩ¯␣␰“¯␰ ͒Ͻ0

(2.12)

anywhere, the ﬂuid is unstable. If the resulting convection is sufﬁciently vigorous, it will homogenize both ¯S and ¯␰ everywhere, except in boundary layers.10
If ¯␣S and ¯␣␰ are constants, we may rewrite inequality (2.12) as

¯g•“C¯ Ͻ0, where
CϭϪ␣SSϪ␣␰␰

(2.13) (2.14)

is the CO density, which stands for the convection-
originating density (although strictly speaking it is a relative density).11 This name underscores its virtue: the

10The scenario envisaged here for the core is very different

from that of the oceans, where Ϫ␰ would correspond to the salt

content of sea water. The resulting thermohaline convection is

usually modeled using a Boussinesq approximation, so that in

place of Eq. (2.12), for instance, the criterion for instability is

Ϫ¯g•

(

¯␣

“

T¯

ϩ

¯␣

␰ T

“¯␰

)

Ͻ

0,

where

␣

␰ T

is

the

isothermal

coefﬁcient

of compositional expansion; see Sec. IV.F. Vertical mixing in

the oceans is not sufﬁciently effective to homogenize T¯ and ¯␰.

Similarly, a stably stratiﬁed layer may exist at the top of the

FOC (Braginsky, 1999). 11The CO density has been widely used for convection in

Boussinesq (i.e., almost incompressible) ﬂuids, but Braginsky

and Roberts (1995), who invented the acronym, seem to have

been the ﬁrst to show that the same idea works for compress-

ible ﬂuids. For a discussion of the history of the idea, see Bra-

ginsky and Roberts (2000).

buoyancy force arising from the density differences created by changes in S and ␰ produces convective circulations, but those created by a change in P do not; its buoyancy force can be absorbed into the gradient of a reduced pressure (see Sec. IV.C). The CO density neatly splits off the important from the unimportant. The rate at which buoyancy supplies energy to the ﬂow is

͵ QCϭ CV•¯g¯␳d3x. FOC

(2.15)

It is obvious from inequality (2.12) that, if d¯␰/dr is

sufﬁciently negative, the ﬂuid is unstable even when

d¯S/dr is positive. And, as Loper (1978) observed, if the

outward

compositional

ﬂux

I

␰ r

is

sufﬁciently

large

and

positive, convection will occur through compositional

buoyancy, irrespective of whether the outward heat ﬂux

I

q r

exceeds

the

adiabatic

ﬂux

¯I

ad r

or

not.

The

heat

QCMB

leaving the core may be greater or less than the heat

HCMB conducted down the adiabat; compositional buoy-

ancy can be so strong that vigorous mixing occurs that

homogenizes ¯S and ¯␰ even though heat is being pumped

downwards by hot descending ﬂuid and cold rising cur-

rents! Nevertheless, most simulations of core convection

and the geodynamo assume that the convective heat ﬂow QCcoMnvB at the core-mantle boundary is positive (but see Glatzmaier and Roberts, 1997). The reason is

simple: a downward entropy ﬂux reduces the vigor of

convection and is detrimental to dynamo action.

The mantle and core should ideally be studied as a

coupled system. This, however, is impractical; because of

their very different viscosities, the mantle and core op-

erate on vastly different time scales. Usually therefore

they are considered separately, the thermal state of the

one providing a boundary condition at the CMB for the

other. The mobility of the core allows it to carry heat

readily from one area of the CMB to another, so equal-

izing their temperatures; this boundary is therefore an

isothermal surface. Its temperature prescribes the bot-

tom boundary condition for mantle convection. The heat ﬂux ICq MB(␪,␾) per unit area is then found as a function of (␪,␾) by solving the mantle convection equations.12 This ICq MB(␪,␾) prescribes the upper
boundary condition on the solution of the core convec-

tion equations. In this way, convection in the mantle is

affected by heating from the core, and the core is cooled

by the overlying mantle. Convection in the mantle re-

sembles

a

valve

controlling

the

heat

ﬂux

I

q CMB

from

the

core. Core convection and the geodynamo are com-

pletely at the mercy of this valve. Unless the valve is set

to allow enough heat through (i.e., unless mantle con-

vection is sufﬁciently strong), the dynamo will fail and

12Cold subducting lithospheric plates may descend to a graveyard of plates at the base of the mantle where, until their temperature rises, they cause the heat ﬂux from the core to the mantle to increase locally. See, for example, Ricard et al. (1993).

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1092

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

core convection may even shut down. In the reference
frame moving with the mantle, the setting of the valve
changes only very slowly, on the convective time scale ␶mϳ108 yr of the mantle. The character of core convection and the geodynamo respond on the same ␶m time scale. This provides the most plausible explanation of
the variations in reversal frequency shown in Fig. 5.

D. Thermodynamic efﬁciency

Even if energy is available to satisfy the balance (2.3), is it obvious that the dynamo will function? A further complication was ﬁrst noticed by Braginsky (1964d), who, regarding the dynamo as a heat engine, questioned whether enough buoyant power would be available to maintain the ﬁeld. The issue of thermodynamic efﬁciency was taken further by several authors, most recently by Braginsky and Roberts (1995), whose article gives references to earlier work. This subsection aims to estimate the efﬁciency of the geodynamo.
There are signiﬁcant points of difference between the geodynamo and the heat engine considered in texts on thermodynamics. To understand the former better, let us start from the classic heat engine operating steadily in the well-known Carnot cycle. The heat input Qin is provided at a higher temperature Tin than the temperature Tout at which heat Qout is extracted. In a perfect machine, the difference QinϪQout is all useful work, but in practice a part QD of this energy is dissipated uselessly as heat. The rate at which the machine does useful work is therefore

Aϭ QinϪ QoutϪ QD .

(2.16)

The Efﬁciency of the engine ␩E is deﬁned as

␩Eϵ

A Qin

ϭ1

Ϫ

QQoiuntϪ

QD Qin

.

(2.17)

For simplicity, let us suppose that QD is expended at a
single moment during the cycle, when the temperature
TD is within the operating temperature range: Tout ϽTDϽTin of the machine. The entropy balance is then

Qin T in

ϭ

Qout T out

ϩ

QD TD

.

(2.18)

Equations (2.16)–(2.18) imply that

␩EϭfF␩C ,

(2.19)

␩Cϭ1ϪTout /Tin .

(2.20)

Here ␩C is the Carnot efﬁciency and fFϭ1 Ϫ(QD/Qin␩C)(1ϪTout /TD) is a frictional factor, lying between 0 and 1 and representing a reduction in the
efﬁciency of the machine below the theoretical maximum ␩C .
We now abandon the classical heat engine and at-
tempt to apply similar ideas to Earth’s core. There are
new and signiﬁcant features. First, the deﬁnition of the ‘‘useful work’’ A done by the engine is largely arbitrary. Let us suppose that it is the rate of production of large-
scale magnetic energy; we do this even though that en-

ergy is eventually ohmically degraded into heat at the same rate QJ. The engine must also make good the remaining frictional losses, QF (say), such as Q␯ and Qt in
Eq. (2.2). It follows that

Aϭ QJ ,

(2.21)

QD ϭ QJ ϩ QF .

(2.22)

A second difference is that, since both QJ and QF reappear within the ﬂuid, they must be regarded as part of the energy source driving the engine. The energy balance replacing Eq. (2.16) is therefore simply

QinϭQout .

(2.23)

The entropy balance is expressed by

Qin T in

ϩ

QD TD

ϭ

Qout T out

.

(2.24)

To be slightly more sophisticated, we may suppose that QJ is dissipated at one temperature TJ and QF at an-
other, TF ; then TD is a compromise between TJ and TF deﬁned by

QD QJ QF TD ϭ TJ ϩ TF .

(2.25)

By Eq. (2.23), we may rewrite Eq. (2.24) as

ͩ ͪ QDϭQinTD

11

T

Ϫ
out

T in

.

(2.26)

It may be seen from Eqs. (2.23) and (2.24) that heat is needed not to maintain the energy balance but to preserve the entropy balance. This limits the efﬁciency of the device in producing magnetic power. It suggests that, in analogy with the oft encountered phrase ‘‘the available energy,’’ it is useful to call QD ‘‘the available dissipation,’’ of which only the fraction fF goes into useful work, QJ.
The dynamo Efﬁciency is

QJ QJ QD ␩Eϵ Qin ϭ QD • Qin ϭfF␩I ,

(2.27)

where the ‘‘frictional factor,’’ fFϭQJ/(QJϩQF), is the fraction of QD that is ‘‘useful dissipation,’’ while the fac-

tor

ͩ ͪ 1 1

␩IϭTD

Ϫ T out T in

(2.28)

is the Ideal efﬁciency, which cannot be exceeded, even
when there is no internal friction. Since ToutрTD рTin , it follows from Eq. (2.28) that

␩Cр␩Iр␩B , where ␩Bϭ͑ Tin /Tout͒Ϫ1. (2.29)

The argument is now further modiﬁed by recognizing the third signiﬁcant difference between the classical heat engine and the Earth’s core: the former is a machine that works steadily, but the core is a slowly evolving system that, apart possibly from some radioactivity, receives all its energy from cooling and gravitational settling. We shall conﬁne this nonstationarity to the refer-

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1093

ence state and shall consider the superimposed
convection as cyclic, i.e., one that, when averaged over the time scale ␶c of the convection, varies only on the time scale of the slow evolution of the core ¯␶. We deﬁne the ideal dynamo efﬁciency ␩D as

␩DϭQJ/QCMB .

(2.30)

The reciprocal St of ␩D was introduced by Stevenson (1984), who pointed out that StϾ1 is a necessary condition for dynamo action.
In considering the average entropy balance of the FOC, we deﬁne [in the spirit of our earlier discussion; see Eq. (2.25)] effective temperatures TR and TD for radioactivity and dissipation:

͵ ͵ QD/TDϭ

͑ qD/T¯ ͒d3x, QR/TRϭ

͑ qR/T¯ ͒d3x,

FOC

FOC

(2.31)

where qD and qR are the volumetric sources that give rise to QD and QR on integration; see Sec. IV. We also encounter the term

͵ Sថ ¯␳ FOC

d

3

x

ϭ

MFOCSថ

ϭ

Ϫ

QS T¯ 0

,

(2.32)

by Eq. (2.10). The entropy balance then assumes the form

QCMB ϭ HICBϩQL ϩ⌺ϩ QD ϩ QS ϩ QR .

T CMB

T ICB

TD T¯ 0 TR

(2.33)

The left-hand side is the entropy ﬂow out of the ﬂuid outer core across the CMB and the ﬁrst term on the right-hand side is the ﬂow into the ﬂuid outer core across the ICB; ⌺ is the entropy source from conduction of heat down the adiabat:

͵⌺ϭ

K¯ M͑ “T¯ /T¯ ͒2d3x.

FOC

(2.34)

From the model described in Sec. IV.B below, we ﬁnd that

⌺Ϸ190 MW KϪ1.

(2.35)

According to Eqs. (2.9) and (2.33), we have QD ϭ QG ϩ QH ,

(2.36)

where

ͩ ͪ ͩ ͪ QHϭ͑HICBϩQL͒

1Ϫ TD T ICB

ϩ QCMB

TD Ϫ1 T CMB

ͩ ͪ ͩ ͪ ϩQS

1Ϫ TD T¯ 0

ϩ QR

1Ϫ TD TR

ϪTD⌺.

(2.37)

Taking T¯ 0 as a rough but plausible estimate of TD (see Table I and Sec. II.A) and assuming that QRϭ0, we obtain

QHϷ0.6 TW, QDϷ1.9 TW.

(2.38)

Thus all the gravitational energy released is available for dissipation but (as anticipated from the earlier discus-

sion) only a part, proportional to ⌬T/T for some ⌬T, of each thermal contribution to QD is available; in the case

of ⌺ this is evident from its deﬁnition (2.34). Recalling that Q␯ is probably an insigniﬁcant part of QD, we see

from Eqs. (2.3) and (2.38) that the turbulent dissipation

is about six times the large-scale Joule losses (2.1).

By Eq. (2.36), the efﬁciency (2.30) of the dynamo is ␩DϭfF␩G where ␩G is the ideal Geodynamo efﬁciency,

QD ␩Gϵ QCMB ,

(2.39)

and fFϭQJ/QD is the frictional factor. Our estimates of these are

␩GϷ28%, ␩DϷ4%, fFϷ0.15.

(2.40)

A similar conclusion follows from an expression derived

by Braginsky and Roberts (1995):

ͫ ͬ 1
␩GϷ QCMB

QGϩ

⌬T0 T CMB

͑

QCcoMnvBϩ

HICBϩ

QL

͒

,

(2.41)

where ⌬T0ϭT¯ 0ϪTCMB ; this is also about 28%. The heat conducted down the adiabat does not drive convec-
tion and is conspicuous by its absence from Eq. (2.41).
The two remaining thermal terms are also diminished in usefulness by the factor ⌬T0 /TCMBϷ0.13.
The importance of the issue of thermodynamic efﬁ-
ciency was ﬁrst raised by Braginsky (1964d), who argued
that the efﬁciency of a thermally driven dynamo would be proportional to ⌬T0 and therefore small. He also foresaw that the compositionally driven dynamo would
be 100% efﬁcient. This led him to the conclusion that
the geodynamo is primarily compositionally driven. This
is supported by our present estimates, which suggest that QG is about twice QH.
From an assumed QCMB , the efﬁciency ␩G determines QD and therefore an upper bound on QJ. Two other examples of this are the following: (a) If QCcoMnvBϭQRϭ0, then ␩GϷ24% and QDϷ1.3 TW; (b) If ¯␶ϭ4 Gyr and QRϭ0, thermal buoyancy opposes compositional buoyancy, and ␩GϷ2%, which implies that QDϷ0.04 TW. Example (b) touches on a perplexing question of how a
thermally driven geodynamo, which is the only type of
dynamo that can operate before the birth of the solid
inner core, can maintain a ﬁeld having an intensity much
the same as at present, as the paleomagnetic evidence
requires (Sec. II.C). Are the estimates made above so greatly in error? Is the assumption that QR was negligible in the remote past (or even now) unjustiﬁed?

III. BASIC DYNAMO THEORY
A. The induction equation
This section describes solutions of ‫ ץ‬t Bϭ “ Ã͑ VÃBϪ ¯␩ “ ÃB͒ , “ • Bϭ 0.

(3.1) (3.2)

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1094

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

These are the equations governing the magnetic ﬁeld B in the conducting volume V (ϭFOCϩSIC) of the core. Equation (3.1) is called the induction equation. In this subsection we sketch its origins. Further details can be found in most texts on MHD.
We start from the pre-Maxwell equations, i.e., the Maxwell equations with the displacement current neglected because U¯ is tiny compared with the speed of light. These include Eq. (3.2) together with

‫ ץ‬t Bϭ Ϫ “ ÃE,

(3.3)

Jϭ“ÃB/␮0 ,

(3.4)

where E is the electric ﬁeld. Implicit in Eq. (3.4) is the assumption that the permeability ␮0 is that of free space, so that the magnetizing force H is simply B/␮0 . The source of B and E is the electric current density J, which
for a dense ﬂuid such as the Earth’s core is given by
Ohm’s law, in the form appropriate for a moving, dense,
electrical conductor:

Jϭ␴M͑ EϩVÃB͒,

(3.5)

where ␴M is the electrical conductivity. Equations (3.4) and (3.5) imply13

Eϭ Ϫ VÃBϩ ¯␩ “ ÃB,

(3.6)

where (see Sec. I.B) ¯␩ϭ1/␮0␴M . Equation (3.1) follows at once from Eqs. (3.3) and (3.6) and, if ␴M is constant,
it may be written alternatively as

‫ץ‬tBϭ“Ã͑ VÃB͒ϩ¯␩ٌ2B.

(3.7)

Pre-Maxwell theory is Galilean invariant and not Lorentz invariant; B, and therefore J, are frame invariant, though plainly E is not. Also, although ‫ץ‬t and V depend on the reference frame, Eq. (3.1) does not.
For simplicity, let us ignore electrical conduction everywhere in the exterior Vˆ of V, including the mantle. The electromagnetic ﬁeld again obeys Eqs. (3.2)–(3.4), but Eq. (3.5) is replaced by Jˆ ϭ0, so that Bˆ is a potential ﬁeld [see Eqs. (1.15) and (1.16)]. And

BϭBˆ , on the core-mantle boundary.

(3.8)

13This form of Ohm’s law highlights the fact that, in contrast to full Maxwell theory, E and B are not on an equal footing in the pre-Maxwell approximation: B is the master and E the slave. If desired, E can be obtained from Eq. (3.6) after B has been determined. The energy density of the electric ﬁeld is negligible compared with that of the magnetic ﬁeld, and the electric stresses are negligible compared with the magnetic stresses. The free charge density is in general nonzero but exerts a force on the conductor that is negligible compared with the Lorentz force. If ¯␩ and/or rÃV is assumed to be discontinuous on the inner core boundary, E makes a ﬂeeting reappearance, since it is necessary to ensure that rÃE is continuous. Because of the surface charge density generally present on a surface of discontinuity of material properties, there is a jump in the normal component of E there, but (since ¯␩ 0) all components of B are continuous.

B. Kinematic dynamos, Cowling’s theorem

The dynamos we study in this review are often called ‘‘homogeneous dynamos,’’ to distinguish them from the manmade dynamos that are deliberately constructed to be inhomogeneous and that obviously work. It is not so evident that a continuous simply connected mass of ﬂuid can create ﬁeld efﬁciently (or at all).
The dynamo problem arises when we require that there be no sources of magnetic ﬁeld in Vˆ even ‘‘at inﬁnity.’’ This demand is crucial, since it is always possible to maintain a ﬁeld in a conductor by applying one from the outside. Sources in Vˆ at a ﬁnite distance from the conductor, arising from electric currents or permanent magnetism, are excluded by Eqs. (1.15) and (1.16), and sources at inﬁnity are eliminated by the requirement (1.17) that the ﬁeld be dipolar at great distances:

Bˆ ϭO͑ rϪ3͒, as r→ϱ.

(3.9)

We shall refer to this as the ‘‘dynamo condition.’’ A successful dynamo is one that actively maintains B for as long as energy sources exist to maintain V:

B→” 0, as t→ϱ.

(3.10)

In practical terms, t→ϱ can simply mean tӷ␶␩ (see below).

There are two kinds of dynamos: the kinematic dynamo and the MHD dynamo.14 The kinematic dynamo

problem is, Given V, ﬁnd B. The MHD dynamo prob-

lem is, Given an energy source such as buoyancy, ﬁnd B

and V. Evidently, the MHD problem contains the kine-

matic problem. The kinematic problem is linear in the

unknown B and is therefore considerably easier to solve

than the MHD problem, which is nonlinear in B and V.

The linearity of kinematic dynamo theory is also its prin-

cipal weakness. It predicts that B either grows without

limit for the given V, or dies inexorably, or rests on a

knife edge between. This unreality disappears when mo-

mentum balance is demanded. The Lorentz force JÃB regulates V so that the knife edge becomes the norm for

a successful dynamo, at least when averaged over time.

We may regard kinematic theory as being realistic for

studying the growth of a weak ‘‘seed ﬁeld’’ until (or if)

that seed ﬁeld becomes so strong that its Lorentz force

alters V.

When V is time independent, the kinematic dynamo

deﬁnes a linear eigenvalue problem for the growth rate ␭ of the ﬁeld ͓B(t)ϭB(0)exp(␭t)͔. Because the eigenvalue problem is not self-adjoint, ␭ is usually complex. If

Im(␭) 0, the solution oscillates, corresponding to ﬁeld
reversal, but on a time scale of order ␶␩ , which is too short to have a direct bearing on the geomagnetic ﬁeld

(although it is of interest in the modeling of magnetic

activity in the Sun and other lower main-sequence stars).

14Sometimes the MHD dynamo is called the ‘‘self-consistent dynamo.’’ We refrain from using this term, to avoid the unfortunate implication that the kinematic theory is inconsistently formulated.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1095

By condition (3.10), dynamo action requires that Re(␭) у0 for at least one eigenvalue, with Re(␭)ϭ0 for that
eigenvalue deﬁning the marginal state. There is no guarantee that, if V is a dynamo motion, ϪV is one also.
When V is periodic in time with period P, dynamo action requires that a solution exist in which B(tϩP) ϭB(t)exp(␭P) with Re(␭)у0. When V is time dependent but not periodic (as can happen when solving the MHD problem—see Sec. VI), the claim of dynamo action is not so easily made precise but, if B has persisted without noticeable diminution over a time long compared with ␶␩ (so that transients from the initial state have disappeared), it is reasonable to claim that a dynamo has been found. A less stringent deﬁnition of dynamo action has been proposed by Hughes (1993). But both deﬁnitions of a dynamo differ from what is usually meant by this term when describing the reversed-ﬁeld pinch in plasma research (see, for example, Fowler, 1999; Kabantsev et al., 1999).
In this section we consider only kinematic theory; the remainder of the article deals with MHD dynamos. By taking the scalar product of Eq. (3.1) with B/␮0 and integrating over the core V, we ﬁnd that

Ϫ ‫ ץ‬t EB ϭ QB ϩ QJ ,

(3.11)

where

͵ ͵ EBϭ

2

1 ␮

0

B2d3x, QBϭ
Vϩ Vˆ

VV•͑ JÃB͒d3x,

͵ · QJϭ␮0

¯␩ J2 d 3 x .
V

(3.12)

Equation (3.11) has a simple interpretation: the lefthand side is the rate of decrease of magnetic energy EB; the ﬁrst term on the right-hand side is the rate QB at
which the Lorentz force converts magnetic energy into
kinetic energy; the ﬁnal term is the Ohmic dissipation QJ, which is positive. In a successful dynamo, the ﬁrst
term on the right-hand side is negative and creates mag-
netic energy from kinetic energy at a rate greater than the rate QJ at which electrical resistance can transform it
into heat. By dimensional analysis of Eq. (3.1), we see that a necessary condition for this to happen is15 Rm
տO(1), where

RmϭU¯ L¯ /¯␩

(3.13)

is the magnetic Reynolds number. According to estimates (1.24) and (1.25), this is about 125 for the Earth. A long-standing challenge to experimenters has been to exhibit convincingly a homogeneous ﬂuid dynamo in the laboratory where, even with the use of liquid metals, O(1) values of Rm are hard to attain because of the comparatively small values of L¯ available; see Roberts

and Jensen (1993), Busse (2000). Very recent reports suggest that this challenge has now been met; see Gailitis et al. (2000), Stieglitz and Mu¨ ller (2000).
Unfortunately, the condition RmտO(1), though necessary, is not sufﬁcient for dynamo action, and it is easy to construct examples in which Rm is large and the magnetic energy EB grows strongly at ﬁrst, but ultimately tends to zero. Cowling (1933) was responsible for a major setback in the subject:
Cowling’s theorem: an axisymmetric B cannot be sustained by dynamo action. [A ﬁeld is said to be axisymmetric if it is the same in every meridional plane about the symmetry axis; this does not mean that its zonal component is zero. In our application, the polar axis Oz deﬁnes axisymmetry; see Eq. (3.20) below.] It took 25 more years before the existence of homogeneous dynamos was unequivocally established. By now there are several simple examples, based on the observation that, although Cowling’s theorem rules out axisymmetric B, it does not rule out kinematic16 dynamos with axisymmetric V. The simplest example of all (Ponomarenko, 1973) operates in a spaceﬁlling conductor which is stationary apart from an inﬁnite cylinder (sϽa) that is in helical motion:

Vϭ␻s1␾ϩU1z ,

(3.14)

where (s,␾,z) are cylindrical coordinates, unit vectors

in the direction of coordinate q increasing are denoted by 1q , and U and ␻ are nonzero constants. Ponomarenko’s model is obviously remote from geophysics, but it does dispel any gloomy thoughts Cowling’s theorem might provoke. And simple axisymmetric models have been devised (Dudley and James, 1989) that resemble the Ponomarenko model but that ﬁt into a sphere. The Ponomarenko model and the Dudley and James models maintain asymmetric B, thus evading Cowling’s theorem. The ﬁeld travels as a wave around the axis of symmetry, a wave that grows in amplitude if Rm exceeds a critical value Rmc but that diminishes to zero if Rm ϽRmc .
With the Earth in mind and with symmetry dictated by the polar axis (because of the importance of Coriolis forces; see Sec. V), we should focus on working dynamos in which B possesses a nonzero axisymmetric part B˜ . In what follows, we shall write

BϭB˜ ϩBЈ, VϭV˜ ϩVЈ,

(3.15)

and deﬁne the emf (electromotive force) due to the asymmetries by

EϵVЈÃBЈϭE˜ ϩEЈ.

(3.16)

Dividing Eq. (3.7) into axisymmetric and asymmetric parts, we obtain

15A fascinating branch of kinematic dynamo theory concerns the fast-dynamo limit, Rm→ϱ. The geodynamo operates on
the ‘‘slow’’ diffusive time scale ␶␩ , and we shall therefore not discuss fast dynamos, but see Childress and Gilbert (1995).

16It does, however, prevent MHD dynamos from being axisymmetric in V. Kinematic dynamos with axisymmetric V create asymmetric B. The concomitant asymmetry in the Lorentz force adds asymmetry to V.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1096

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

‫ץ‬tB˜ ϭ“Ã͑ V˜ ÃB˜ ϩE˜ ͒ϩ¯␩ٌ2B˜ ,

(3.17)

‫ץ‬tBЈϪ“Ã͑V˜ ÃBЈϩEЈ͒Ϫ¯␩ٌ2BЈϭ“Ã͑VЈÃB˜ ͒. (3.18)
A transparent way of establishing Cowling’s theorem was devised by Braginsky (1964a) and requires V˜ to be an incompressible ﬂow (“•V˜ ϭ0).17 First V˜ and B˜ are divided into their (axisymmetric) zonal and meridional parts,

V˜ ϭs␨͑ s,z,t ͒1␾ϩVp͑ s,z,t ͒,

(3.19)

B˜ ϭB͑ s,z,t ͒1␾ϩBp͑ s,z,t ͒,

(3.20)

where p is used for the ‘‘meridional’’ or ‘‘poloidal’’ component; this can be represented by a vector potential,
e.g.,

Bpϭ“Ã͓A͑ s,z,t ͒1␾͔.

Then Eq. (3.17) gives

ͫ ͬ 1
s

‫ץ‬ ‫ץ‬t ϩVp•“

sAϪ¯␩⌬Aϭ˜E␾ ,

ͫ ͬ ‫ץ‬
s ‫ץ‬t ϩVp•“

B s

Ϫ

¯␩

⌬

B

ϭ

s

Bp•

“

␨

ϩ

͑

“

ÃE˜ ͒

␾

,

(3.21) (3.22) (3.23)

where ⌬ϭٌ2ϪsϪ2. In Vˆ , we have

⌬Aˆ ϭBˆ ϭ0, Aˆ →0 as r→ϱ,

(3.24)

and Eq. (3.8) requires that

AϭAˆ ,

‫ץ‬A ‫ץ‬Aˆ ‫ץ‬r ϭ ‫ץ‬r ,

BϭBˆ ,

on

the

CMB.

(3.25)

If Eϵ0, it follows from Eqs. (3.22)–(3.25) that

͵ ͵ 1
2

‫ץ‬t

͑ sA ͒2d3xϭϪ¯␩
V

͓
Vϩ Vˆ

“

͑

s

A

͒

͔

2

d

3

x

Ͻ

0,

(3.26)

͵ ͩ ͪ ͵ ͫ ͩ ͪͬ 1
2

‫ץ‬t

V

B s

2
d 3 x ϭ Ϫ ¯␩

V“

B s

2
d3x

͵ϩ

V

B s2

͓“͑ sA ͒Ã“␨͔␾d3x.

(3.27)

By Eq. (3.26), A→0 as t→ϱ, so that the second term on the right-hand side of Eq. (3.27) ultimately disappears
too. The right-hand side of Eq. (3.27) is then negative, so that B→0 also. In short, Eϵ0 implies that B˜ →0 as t→ϱ. This is Cowling’s result: A working dynamo necessarily maintains a B having an asymmetric part BЈ. Perhaps this is why the geomagnetic and geographic
axes are persistently inclined to one another (see Sec.
I.C).
This proof of Cowling’s theorem shows that to establish dynamo action it is never sufﬁcient to show that EB

17This assumption was not made by Cowling in one of his original proofs, but this proof was not totally satisfactory. Ivers and James (1984) provided the mathematically most general demonstration of Cowling’s theorem.

initially increases in time. Equation (3.27) demonstrates
that, until A becomes small, the rate of increase of the
magnetic energy stored in the zonal ﬁeld may be positive if ␨ is large enough. The term sBp•“␨ in Eq. (3.23) responsible for this can be interpreted by
Alfve´ n’s frozen-ﬂux theorem: Flux tubes in a perfectly
conducting ﬂuid are carried by the ﬂuid in its motion, just
as though they were frozen to it.
The freezing is imperfect when Rm is large but not inﬁ-
nite. The theorem nevertheless provides a useful and
qualitatively correct way of visualizing induction pro-
cesses, except in regions where the ﬁeld gradients are
large and where a locally deﬁned Rm is uncharacteristi-
cally small, thus allowing the severing and reconnection of ﬁeld lines that would be forbidden when Rmϭϱ. In the present case, the zonal shear ␨ stretches the lines of force (sAϭconst) of the meridional ﬁeld Bp along lines of latitude to create B˜ ␾ . Often (though not here) V˜ ␾ /s is denoted by ␻ rather than ␨, and this mechanism for creating B became known as the ‘‘␻ effect’’.

C. Turbulent helicity and the ␣ effect

For given B˜ and V, Eq. (3.18) is an inhomogeneous linear equation for BЈ, its right-hand side being then a known source. Solving Eq. (3.18) for BЈ, we can then obtain E˜ from Eq. (3.16). Thus BЈ(B˜ ) and E˜ (B˜ ) are linear functionals of B˜ , and Eq. (3.17) becomes a closed equation for B˜ .
At ﬁrst sight, this method of solving the kinematic dynamo does not recommend itself. The functionals BЈ(B˜ ) and E˜ (B˜ ) depend on B˜ and V at all points x in the ﬂuid and at all earlier times t. To determine them seems to be a task every bit as daunting as solving Eq. (3.1) itself. Simpliﬁcations arise, however, when the conducting ﬂuid is in turbulent motion: we separate V and B, not as in Eq. (3.15) according to their symmetry, but into their large-scale (which we often call the ‘‘macroscale’’ below) and small-scale (‘‘microscale’’) parts, the former being the average over the turbulence (denoted by angle brackets) and the latter being the ﬂuctuating remnant, for example,

Vϭ͗V͘ϩVt, where ͗Vt͘ϭ0.

(3.28)

In this subsection (but not later) we shall have particu-
larly in mind situations of large magnetic Reynolds number, RmtϭLtV t/␩; these are relevant in many astrophysical contexts.
The emf due to the turbulence is similar to Eq. (3.16):

FϵVtÃBtϭ͗F͘ϩF t.

(3.29)

Separating Eq. (3.7) into its macroscale and microscale parts, we obtain, as in Eqs. (3.17) and (3.18),

‫ ץ‬t ͗ B͘ ϭ “ Ã͑ ͗ V͘ Ã͗ B͘ ϩ ͗ F͘ ͒ ϩ ¯␩ ٌ 2 ͗ B͘ ,

(3.30)

‫ץ‬tBtϪ“Ã͑͗V͘ÃBtϩF t͒Ϫ¯␩ٌ2Btϭ“Ã͑ VtÃ͗B͘ ͒. (3.31)

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1097

FIG. 7. Cartoons describing the microscopic ␣ effect. (a) rising (left) and sinking (right), eddies in the northern hemisphere of a rotating, turbulently convecting ﬂuid; (b) the effect of these on a horizontal ﬁeld line; (c) the situation after the ﬂux loop created shown in panel (b) has broken away from the parent ﬁeld line.

We would like to be able to solve Eq. (3.31) for Bt and thence to obtain ͗F͘ as a linear functional of ͗B͘ that
could be used to close Eq. (3.30). As before, this is a
daunting task. One may, however, argue heuristically that, if the microscale Lt is sufﬁciently small compared with the macroscale L¯ , the functional ͗F͘(͗B͘) at the point x should be well represented by the ﬁrst few terms
in its Taylor expansion about x, e.g.,

͗ F͘ i ϭ ␣ i j ͗ B j ͘ ϩ ␤ i j k ٌ j ͗ B k ͘ ,

(3.32)

where ␣ij and ␤ijk depend on ͗V͘ and on the statistical properties of the turbulence Vt. According to Eq. (3.32), ͗F͘ at x depends only on ͗B͘ and its derivatives at x. This is what is meant by a local turbulence theory; we shall meet it again in Sec. IV.D. The existence and importance of the ﬁrst term in the expansion (3.32) was ﬁrst noticed by Parker (1955). We brieﬂy review his picture of microscale induction by ‘‘cyclonic turbulence.’’
Consider (in the corotating frame) turbulent convection in the northern hemisphere of a rotating sphere of compressible ﬂuid rotating about an axis Oz. To sim-

plify the discussion, suppose that ͗V͘ϭ0. Visualize, as in

Fig. 7(a), the turbulence as consisting of tiny rising and

sinking ‘‘eddies,’’ each having a brief identity before dis-

solving back into its surroundings. The left-hand eddy in

Fig. 7(a) is rising (Vrt Ͼ0), and as it rises it expands,

tending to conserve its angular momentum about the

vertical as it does so. The vertical component (“ÃVt)r

of its vorticity (relative to the rotating frame) is there-

fore negative, as is Vrt (“ÃVt)r . The right-hand eddy in Fig. 7(a) is sinking (Vrt Ͻ0), and as it does so it com-

presses and tends to rotate more rapidly about the ver-

tical

͓ ( “ ÃVt ) r Ͼ 0 ͔ ,

so

that

again

V

t r

(

“

ÃVt

)

r

Ͻ

0.

Aver-

aging over rising and falling eddies, we see that the

helicity18 of the turbulence, deﬁned by

H ϭ ͗ Vt • “ ÃVt ͘ ,

(3.33)

is negative.
In understanding the inductive effects of an eddy, we suppose that its magnetic Reynolds number Rmt is large enough for Alfve´ n’s frozen ﬂux theorem to be useful.
The effect of the cyclonic motions on a large-scale horizontal ﬁeld ͗B͘ is sketched in Fig. 7(b). An upwardmoving eddy makes an ⍀-shaped indentation in a horizontal ﬁeld line and the vertical vorticity simultaneously twists that ⍀ out of the plane of the paper. Now recall that, though Rmt is large, it is not inﬁnite, so that diffu-
sion acts particularly strongly where the ﬁeld gradients are large, as at the base of the ⍀. This causes the ⍀ to be severed from its parent line, to form a loop of ﬂux in a
plane perpendicular to the paper, as in Fig. 7(c). A
downward-moving eddy creates a ﬂux loop of the same
type. These processes are repeated throughout the ﬂuid,
and their net effect is to impose on the initial ﬁeld a
right-handed helical structure, just as though an electric ﬁeld existed to drive current parallel to the initial ﬁeld19:

͗ F͘ ϭ ␣ ͗ B͘ .

(3.34)

Steenbeck and Krause (1966) christened this the ‘‘alpha
effect,’’ for no better reason than that they used the letter ␣ to describe the proportionality. But the name has stuck. It has spawned an entirely new ﬁeld of study,
mean-ﬁeld electrodynamics, which investigates in greater depth how ͗F͘ depends on the strength and statistical properties of the turbulence [see Moffatt (1978) and Krause and Ra¨ dler (1980)].
On comparing the Ansatz (3.34) with Eq. (3.32), we see that the latter implies that ␣ijϭ␣␦ij , which requires that the turbulence be pseudoisotropic, i.e., turbulence
having statistics that are independent of direction but

18This felicitous term is due to Moffatt (1969). The quantity itself was ﬁrst introduced by Steenbeck and Krause (1966), who gave it the more forbidding German name Schraubensinn, which might be translated as ‘‘sense of screws.’’
19From HϽ0, we therefore have ␣Ͼ0, and similarly, ␣Ͻ0 in the southern hemisphere where HϾ0. These pseudoscalars tend to have opposite signs. The use of ␣ elsewhere in this review to mean the coefﬁcient of volume expansion should not cause confusion.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1098

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

that are not invariant under coordinate reﬂection. (In
truth, not even the cyclonic turbulence described above is like that.) Isotropy requires that ␤ijkϭϪ¯␩t⑀ijk where ¯␩t is a turbulent magnetic diffusivity. When ¯␩tӷ¯␩, the ␣ effect Ansatz leads, not to Eq. (3.1), but to an induction equation for ͗B͘ of the form

‫ץ‬

t

͗

B͘

ϭ

“

Ã͑

͗

V͘

Ã͗

B͘

ϩ

␣

͗

B͘

͒

Ϫ

“

Ã͑

¯␩

t
“

Ã͗

B͘

͒

.

(3.35)

Cowling’s theorem does not apply to this equation, and in principle it can have nontrivial axisymmetric solutions when ␣ 0.

D. Large-scale helicity and ␣ effect

Although the Earth’s core is undoubtedly turbulent, it is unlikely, in view of the large magnetic diffusivity of the ﬂuid, that Rmt is large. This implies that microscale induction is not very signiﬁcant, and Eq. (3.1) is more realistic than Eq. (3.35). Our reason for discussing cyclonic turbulence was not geophysical realism, but to introduce in a simple way the concepts of helicity and the ␣ effect and to show how they arise naturally in rotating convection. In fact, convective motions of every scale are helical in a rotating ﬂuid, and they create E¯ by something like an ␣ effect, though one that is nonlocal and not simply related to helicity. We believe that the geodynamo (though not every dynamo in nature) is maintained principally by large-scale ﬂows. The successful Ponomorenko motion (3.14) is of large scale and helical, and the same is true of the Dudley and James models.20
The Ansatz

E˜ ϭ␣B˜

(3.36)

[cf. Eq. (3.34)] is not so easily defended when E˜ is produced nonlocally by VЈ and BЈ, but an interesting class of large-scale motions with helicity that generate a local ␣ effect was discovered by Braginsky (1964a, 1964b). He had a very attractive idea: axisymmetric B is impossible
(i.e., requires inﬁnite Rm), but perhaps nearly axisym-
metric dynamos can work if Rm is large enough. He
supposed that VЈ/V˜ ␾ϭO(RmϪ1/2), which produced BЈ/B˜ ␾ϭO(RmϪ1/2) by Eq. (3.18) and therefore an E˜ ␾ of order RmϪ1. He found, in analogy with the turbulent
␣ effect, that E˜ ϭ␣B˜ ␾1␾ , so that Cowling’s theorem is again evaded; nearly axisymmetric ﬂows can maintain a
dynamo in which B˜ p is of order B˜ ␾ /Rm. Braginsky computed ␣ explicitly in terms of the assumed VЈ/V˜ ␾ . In this way he provided the ﬁrst mathematical justiﬁca-
tion of Parker’s ideas and could also compute what
were, in a sense, the ﬁrst 3D kinematic geodynamos
(Braginsky, 1964c). Soward (1972) later made use of asymptotic techniques based on Rm→ϱ to create a
powerful, pseudo-Lagrangian alternative to Braginsky’s

20Strictly speaking, helicity is not essential for a dynamo to function, but it certainly helps; see Gailitis (1993).

method. He rederived Braginsky’s results and found a direct connection between Braginsky’s ␣ and the helicity of the ﬂow; see also Chap. 8 of Moffatt (1978). This approach also explains why Braginsky’s ␣ is local, in the sense deﬁned above, even though VЈ is a large-scale ﬂow; see also Soward (1990).
The Ansatz (3.36) provides a simple and popular way
of constructing 2D mean-ﬁeld dynamos (Roberts, 1972)
by a two-stage process: step 1: A is created from B by the ˜E␾ϭ␣B source in
Eq. (3.22);
step 2: B is created from A through the source sBp•“␨ϩ(“ÃE˜ )␾ϭ͓“(sA)Ã“␨ϩ“Ã(␣“ÃA1␾)͔␾ in Eq. (3.23). Clearly ␣ is essential in step 1, but there are three possibilities in step 2. If (“ÃE˜ )␾ is large compared with sBp•“␨, step 2 relies on ␣ for a second time, so that such models are called ‘‘␣2 dynamos.’’ If sBp•“␨ is large compared with (“ÃE˜ )␾ , then B is created by the ␻ effect, and the model is an ‘‘␣␻ dynamo.’’ When ␣ and ␻ effects are equally signiﬁcant in step 2, it is an ‘‘␣2␻ dynamo.’’ The terms ␣2 dynamo, ␣␻ dynamo, etc. are also used in Sec. VI to describe how a 3D dynamo
model maintains the axisymmetric part B˜ of its ﬁeld B.
Many useful kinematic dynamos operating with 3D
large-scale motions have been produced by numerical
integration, without any appeal to mean-ﬁeld electrody-
namics or the Braginsky-Soward theory; see, for ex-
ample, Love and Gubbins (1996). It is fashionable to say
that the kinematic dynamo problem is solved, despite
the fact that it is still a nontrivial operation to extract B
from an arbitrarily speciﬁed V.

IV. DYNAMICAL THEORY

A. The full (primitive) equations

The main motivation for this section is to explain how the full MHD theory of the core can be reduced to realistic but manageable proportions. A more complete analysis has been given by Braginsky and Roberts (1995). Readers interested only in a basis for numerical simulations may wish to move forward to Sec. IV.E.
The starting point is the set of primitive equations:

‫ ץ‬t␳ ϭ Ϫ “ • ͑ ␳ V͒ , ␳ d t Vϭ Ϫ “ P ϩ ␳ gϪ 2 ␳ ⍀ÃVϩ ␳ FB ϩ ␳ F␯ , ␳ d t S ϭ Ϫ “ • IS ϩ ␴ S , ␳ d t ␰ ϭ Ϫ “ • I␰ ,

(4.1) (4.2) (4.3) (4.4)

“ • Bϭ 0,

(4.5)

‫ ץ‬t Bϭ “ Ã͑ VÃB͒ Ϫ “ Ã͑ ␩ “ ÃB͒ ,

(4.6)

where dtϭ‫ץ‬tϩV•“ is the motional (or Lagrangian) derivative in the chosen reference frame, which is usually ﬁxed to the mantle. Equations (4.1) and (4.2) express conservation of total mass and momentum. The latter is no more than the usual Navier-Stokes equation in the

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1099

rotating frame, with the buoyancy force ␳g and Lorentz force ␳FBϭJÃB included, where B is determined by Eqs. (4.5) and (4.6), the topic of Sec. III. The rotation of
the frame has introduced not only a Coriolis force Ϫ2␳⍀ÃV but also a centripetal acceleration Ϫ⍀Ã(⍀Ãr), which is absorbed into the Newtonian gravitational acceleration to create the ‘‘effective’’ gravitational ﬁeld gϭϪ“U. The effective gravitational po-
tential U satisﬁes

ٌ2Uϭ4␲G␳Ϫ2⍀2,

(4.7)

where G is the constant of gravitation. The Poincare´ force Ϫ⍀˙ Ãr has been ignored (but would have to be restored if we were to include the luni-solar precession; see Sec. II.A). The viscous force is

␳ F␯ ϭ “ • ␲I ␯ ,

where

␲

␯ ij

ϭ

2

␳

␯

M͑

e

ij

Ϫ

1 3

ekk␦ij͒,

(4.8)

where

e

i

j

ϭ

1 2

(

ٌ

i

V

j

ϩ

ٌ

j

V

i

)

is

the

rate

of

strain

tensor.

Equation (4.3) governs the evolution of entropy, IS be-

ing the entropy ﬂux and ␴S the rate of entropy produc-

tion. Equation (4.4) when combined with Eq. (4.1) en-

sures mass conservation for the individual constituents of the alloy; I␰ is a mass ﬂux proportional to the small

difference between the velocity of the light constituent

and that of the heavy.

Energy conservation is expressed by

‫ ץ‬t u totalϩ “ • Itotalϭ q R ,

(4.9)

where utotal is the total energy density and Itotal is the total energy ﬂux, which includes the heat ﬂux Iq. The radioactive heat source qR is the only volumetric source
of energy. Equations (4.1)–(4.8) must imply Eq. (4.9),
which is the case only if

Iq ϭ T IS ϩ ␮ I␰ ,

(4.10)

T ␴ S ϭ q R ϩ q J ϩ q ␯ Ϫ IS • “ T Ϫ I␰ • “ ␮ ,

(4.11)

where

q J ϭ ␮ 0 ␩ J2 ,

q

␯

ϭ

e

i

j

␲

␯ ij

(4.12)

are the Joule heating and the viscous regeneration of
heat; they are both non-negative. The ﬁnal two terms in
Eq. (4.11) arise from diffusion of heat and composition; this combination must be non-negative too: qtу0 where

q t ϵ Ϫ IS • “ T Ϫ I␰ • “ ␮ .

(4.13)

We shall not need the general expressions for IS, I␰, and ␴S arising from molecular diffusion; they may be found
in, for example, Chap. VI of Landau and Lifshitz (1987).
We give below alternative expressions for turbulent dif-
fusion, and this has motivated the otherwise idiosyncratic notation qt in Eq. (4.13).

B. The reference state

The reference state is hydrostatic with uniform S and ␰; see Sec. II.C. Reference state variables carry an overbar. We have

“P¯ ϭϪ¯␳“U¯ ,

(4.14)

“ ¯S ϭ 0,

(4.15)

“¯␰ ϭ 0.

(4.16)

Applying “Ã to Eq. (4.14), we see that ¯␳ϭ¯␳(U¯ ) and that therefore P¯ and all other thermodynamic variables
are constant on each equipotential surface U¯ ϭconst. Since the inner core boundary is a phase boundary, it too is an equipotential surface.21 We shall brieﬂy return
in Sec. IV.D to the full generality of Eqs. (4.14)–(4.16),
but elsewhere shall assume spherical symmetry. All ref-
erence state variables then depend spatially on r alone, the ⍀2 term in Eq. (4.7) is disregarded, and Eqs. (4.14)– (4.16) become

‫ץ‬rP¯ ϭϪ¯g¯␳, ‫ץ‬r¯Sϭ0, ‫ץ‬r¯␰ϭ0,

(4.17)

where ¯gϭϪ¯grϾ0. These imply

‫ץ‬r¯␳ϭϪ¯␳¯g/¯uS2, ‫ץ‬rT¯ ϭϪ¯␣S¯g, ‫ץ‬r␮¯ ϭϪ¯␣␰¯g. (4.18)
Parameters such as ¯␣S and ¯␣␰ depend spatially on r. The reference state used in the simulations to be described in Sec. VI was based on the preliminary reference Earth model (PREM) of Dziewonski and Anderson (1981).
Because the core evolves in time as the Earth cools, the reference state must be continually updated. One should not lose sight of the fact that the reference state is merely a mathematical convenience and has no profound physical meaning. One reference state would not, even in an immobilized core, evolve into another reference state. To see this, we ﬁrst discard ¯I␰ as negligible. Heat conduction down the adiabat in the reference state would then give (see, for example, Landau and Lifshitz, 1987)

¯ISϭ¯Iq/T¯ , ¯IqϭϪK¯ M“T¯ , ¯␴SϭK¯ M“T¯ /T¯ )2ϩqR/T¯ , (4.19)
where radioactive sources qR have been included (in case any exist). The right-hand side of Eq. (4.3) would therefore be

¯␴

S Ϫ

ϭ

¯␴

T Ϫ

ϩ

q

R/

T¯

,

where

¯␴

ϪT ϭ

T¯

Ϫ

1
“

•

͑

K¯

M

“

T¯

͒

.

(4.20)

Since

¯␴

T Ϫ

depends

on

r,

the

assumed

r

independence

of

¯S in the initial reference state would, according to Eq.

(4.3), be lost. [It may be noted for future reference that

¯␴

T Ϫ

is

negative,

even

though

¯␴ S Ͼ 0 ;

this

is

what

we

ex-

pect for a cooling core (Sថ Ͻ0).]

C. The anelastic equations The convection equations are obtained by expanding
about the reference state by writing
21This statement is true, only insofar as the solid inner core can be treated hydrostatically: see Sec. II.B.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1100

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

Sϭ¯SϩSc , ISϭ¯ISϩIcS , ␳ϭ¯␳ϩ␳c ,... .

(4.21)

We have seen in Sec. II.C that Sc /¯S, ␰c /¯␰,... are of order 10Ϫ8. We may therefore conﬁdently linearize the thermodynamic variables about their reference values, as we essentially did in Sec. II.C. All other nonlinearities are retained. Since22 V¯ ϭB¯ ϭ0, we may omit c from V, B, and J.
The core is cooling and (in the absence of qR) it is driven into motion only by that cooling; it is imperative to incorporate this into the theory. Because the evolution of the core is slow, a ready-made mathematical technique can be used: the method of multiple time scales. See, for example, Bender and Orszag (1978), Chap. 11. All variables are functions of a reference time variable ¯t , and all convection variables are also functions of a convective time variable tc , where ¯t /tc ϭO(108). We write

‫ץ‬

t

ϭ

‫¯ץ‬t

ϩ

‫ץ‬

c t

.

(4.22)

Then, for example, since ‫ץ‬tc¯Sϭ0

ϭ

O

(

10Ϫ

8)

‫ץ‬

c t

S

c

,

‫ץ‬

t

S

ϭ

͑

‫¯ץ‬t

ϩ

‫ץ‬

c t

͒͑

¯S

ϩ

S

c

͒

Ϸ

‫ץ‬

tc S

cϩ

Sថ

,

and ‫¯ץ‬t Sc (4.23)

where the overdot is used (exclusively) for ¯t derivatives.
The two terms on the right-hand side of Eq. (4.23) are similar in magnitude, because Sc /¯SϭO(10Ϫ8).
The governing equations (4.1)–(4.6) become

“ • ͑ ¯␳ V͒ ϭ 0, ‫ץ‬t͑¯␳V͒ϩ“•͑¯␳VVϪ␲J ␯Ϫ␲J B͒

(4.24)

ϭ Ϫ¯␳ “ ⌸ ϩ C¯␳ ¯gϪ 2¯␳ ⍀ÃV,

(4.25)

‫␳¯͑ץ‬Sc͒ϩ“•͑¯␳ScVϩIcS͒ϭϪ¯␳Sថ ϩ¯␴ϪS ϩ␴D,

(4.26)

‫␰␳¯͑ץ‬c͒ϩ“•͑¯␳␰cVϩIc␰͒ϭϪ¯␳␰ថ ,

(4.27)

“ • Bϭ 0,

(4.28)

‫ ץ‬t Bϭ “ Ã͑ VÃBϪ ¯␩ “ ÃB͒ ,

(4.29)

where ⌸ϭPc /¯␳ϩUc is the reduced pressure, Uc is the perturbation in gravitational potential created by ␳c ,
and

C ϭ Ϫ ¯␣ S S c Ϫ ¯␣ ␰ ␰ c

(4.30)

is the CO density; see Eq. (2.14). It is of order 10Ϫ8; see Sec. V.B. The Lorentz force ␳FB has been expressed as

the divergence of the magnetic stress tensor:

¯␳ FB ϭ “ • ␲I B ,

where

␲

B ij

ϭ

␮

Ϫ1 0

͑

B

i

B

j

Ϫ

1 2

B

2

␦

i

j

͒

.

(4.31)

[The steps from Eq. (4.2) to Eq. (4.25) are not immedi-

ate; see Braginsky and Roberts (1995).] In addition to

¯␴

S Ϫ

deﬁned

by

Eq.

(4.20),

entropy

is

produced

convec-

tively at the rate

22A velocity V¯ arises through gravitational settling; it is very tiny ͓V¯ ϭO(10Ϫ8Vc)͔, and (like ¯I␰) we ignore it.

␴DϭqD/T¯ , where qDϭqJϩq␯.

(4.32)

The approximation (4.24) for mass conservation is known as the anelastic equation. It ﬁlters out sound (seismic) waves, which are uninteresting since they cross the core in minutes rather than in the decades to millenia of geomagnetic phenomena.

D. Core turbulence

Four molecular diffusivities appear in core MHD: the
kinematic viscosity ␯M , the magnetic diffusivity ␩, the thermal diffusivity ␬M , and the compositional diffusivity DM between the two components of the alloy. (For a recent estimate of ␯M , see de Wijs et al., 1998.) Three analogs of the magnetic Reynolds number (3.13) arise,
namely, the Reynolds number, the Pe´ clet number, and
the mass Pe´ clet number:

U¯ L¯

U¯ L¯

U¯ L¯

Reϭ , Peϭ , Mpϭ .

¯␯ M

¯␬ M

D¯ M

(4.33)

We have discussed the signiﬁcance of Rm in Sec. III.B. Using estimates (1.24) and (1.25), and values listed in Table I, we ﬁnd that

ReϷ2.5ϫ108, PeϷ5ϫ107, MpϷ2.5ϫ1011. (4.34)

Given such giant numbers, it cannot be doubted that small-scale turbulent eddies exist in the ﬂuid outer core, as they do in the Earth’s atmosphere and oceans. Although these small eddies have little energy compared with that in the larger scales, they are, by many orders of magnitude, more effective in transporting heat, composition, and momentum than molecular diffusion. The computing resources that would be necessary to resolve numerically the full spectrum of turbulent length and time scales does not exist, and progress can only be made by being less ambitious. We shall therefore aim to model the macroscales by averaging Eqs. (4.24)–(4.29) over the microscale, a process that introduces turbulent ﬂuxes of macroscale momentum, heat, and composition which have to be parametrized.
We divide each variable into a macroscale and a microscale part, as in Eq. (3.28). The average of the term ¯␳VV in Eq. (4.25) is

͗¯␳ VV͘ ϭ¯␳ ͗ ͑ ͗ V͘ ϩ Vt ͒ ͑ ͗ V͘ ϩ Vt ͒ ͘

ϭ¯␳ ͗ V͘ ͗ V͘ ϩ¯␳ ͗ Vt Vt ͘ .

(4.35)

Here ¯␳͗VtVt͘ is known as the Reynolds stress tensor. As in Sec. III.C, we shall develop a local turbulence
theory, based on the ‘‘Reynolds analogy,’’ the basic idea of which is that the transport of macroscale ﬁelds by chaotic, subgrid-scale ‘‘eddies’’ is similar to their transport by chaotic molecular motions and can therefore be represented mathematically in a similar way though with turbulent diffusion coefﬁcients that are much greater than their molecular counterparts. In this way ‘‘eddy’’ transport coefﬁcients are used to account crudely for the mixing done by the unresolved scales of motion.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1101

Momentum transport provides the most famous appli-

cation of local turbulence theory: since the stresses due

to

molecular

viscosity

¯␯ M

are,

by

Eq.

(4.8),

␲

␯ ij

ϭ¯␳

␯

M

(

ٌ

i

V

j

ϩ

ٌ

jV

i

Ϫ

2 3

ٌ

k

V

k

␦

i

j

)

,

the

last

term

in

Eq.

(4.35) can be analogously represented by

␲

t i

jϵ

Ϫ

͗¯␳

V

t i

V

t j

͘

ϭ¯␳

␯

t͑

ٌ

i

͗

V

j

͘

ϩ

ٌ

j

͗

V

i

͘

Ϫ

2 3

ٌ

k͗

V

k

͘

␦

ij

͒

,

(4.36)

where ␯t is the ‘‘turbulent viscosity.’’ The Reynolds
stresses then depend only on the local gradient of the
macroscale ﬂow. The magnetic diffusivity ␩ greatly exceeds ␯M , ␬M ,
and DM , and the microscale magnetic Reynolds number RmϭU tLt/␩ is small. This means (see Sec. III.D) that molecular diffusion of large-scale magnetic ﬁelds is more signiﬁcant than turbulent diffusion, i.e., when we average VÃB in Eq. (4.29) and obtain

͗ VÃB͘ ϭ ͗ V͘ Ã͗ B͘ ϩ ͗ Vt ÃBt ͘ ,

(4.37)

the ﬁnal term, ͗F͘ϭ͗VtÃBt͘, is small in comparison with ¯␩“Ã͗B͘. This does not mean that the magnetic Reynolds stress, arising from the average of ␲I B in Eq. (4.25), may also be neglected; Bt acts as an additional
brake on the macroscale ﬂow and, when we combine the
turbulent kinetic and magnetic Reynolds stresses and represent them by the ␲I t deﬁned in Eq. (4.36), we should increase ␯t appropriately.
The Reynolds analogy is also applied to two terms
that appear when Eqs. (4.26) and (4.27) are averaged:

IcS t

ϵ

͗¯␳

S

t c

Vt

͘

ϭ

Ϫ¯␳

␬

t

“

͗

S

c͘

,

(4.38)

Ic␰tϵ͗¯␳␰ct Vt͘ϭϪ¯␳␬t“͗␰c͘.

(4.39)

Since turbulence should transport all passive scalars in the same way, the same turbulent diffusivity (␬t) ap-
pears in both expressions. The energy dissipation qD can also be divided into
macroscale and microscale parts:

qDϭ͗qJ͘ϩ͗q␯͘ϩqt,

(4.40)

where ͗qJ͘ϭ␮0¯␩͗J͘2, ͗q␯͘ϭ͗eij͗͘␲i␯j͘,
are the macroscale dissipations and

(4.41)

q

t

ϭ

͗

␮

0

¯␩

͑

Jt

͒

2ϩ

e

it j ͑

␲

␯ ij

͒

t

͘

(4.42)

is the turbulent dissipation. Braginsky and Roberts (1995) showed that, since ¯␯MӶ¯␩, the viscous turbulent dissipation in Eq. (4.42) is negligible compared with the
ohmic.
Braginsky and Roberts (1995; Appendix C) deduced an alternative expression for qt based on Eq. (4.13)
which demonstrated that the energy loss (4.42) can be
replaced directly by the rate of working of the buoyancy
force and not by the turbulent cascade of classical shear
ﬂow turbulence. They showed that

q

t

ϭ¯␳

␬

t¯g•

͑

¯␣

S
“

͗

S

c͘

ϩ

¯␣

␰

“

͗

␰

c

͘

͒

.

(4.43)

If ¯␣S and ¯␣␰ are constants, this can be written as

q t ϭ Ϫ¯␳ ␬ t¯g• “ ͗ C ͘ .

(4.44)

Although ␬t and ␯t depend on position and time

through the local strength of the turbulence, we expect that in the main body of the core they will be O(U tLt). If we take Ltϭ104 m and U tϭ10Ϫ4 m sϪ1 as typical, we obtain ␬tϷ␯tϷ1 m2 sϪ1, which may be compared with

the tiny molecular diffusivities listed in Table I. We may

now, with one exception, remove the molecular ﬂuxes of

momentum, entropy, and composition in favor of the corresponding turbulent ﬂuxes.23 The one exception is

the conduction of heat down the adiabat, which cannot

be ignored because T¯ /͗Tc͘ϳ108.

E. Working equations and boundary conditions

We summarize our ﬁnal, working equations that gov-
ern the large-scale ﬁelds in the Earth’s core, the equa-
tions used in the simulations reported in Sec. VI. The
notation will be simpliﬁed: the angle brackets denoting the mean ﬁelds will be omitted but implied, and ␯ and ␬ will be turbulent diffusivities, no longer carrying a bar or a superscript t; the overbar, and the subscript c , will also be omitted except where confusion might arise. We
have

“ • ͑ ␳ V͒ ϭ 0, d t Vϭ Ϫ “ ⌸ ϩ C gϪ 2 ⍀ÃVϩ FB ϩ F␯ ,

(4.45) (4.46)

␳dt␰cϩ“•I␰ϭϪ␳␰ថ , I␰ϭϪ␳␬“␰c ,

(4.47) (4.48)

␳

d

t

S

cϩ

“

•

IS

ϭ

Ϫ

␳

Sថ

ϩ

¯␴

S Ϫ

ϩ

␴

D

,

(4.49)

ISϭϪ␳␬“Sc ,

(4.50)

“ • Bϭ 0,

(4.51)

‫ ץ‬t Bϭ “ Ã͑ VÃB͒ Ϫ “ Ã͑ ␩ “ ÃB͒ ,

CϭϪ␣SScϪ␣␰␰c ,

␳

F

␯ i

ϭ

ٌ

j␲

␯ ij

,

␲

␯ ij

ϭ

2

␳

␯

͑

e

i

j

Ϫ

1 3

e

k

k

␦

i

j

͒

,

␴DϭqD/T,

(4.52) (4.53) (4.54) (4.55) (4.56)

qDϭqJϩq␯ϩqt,

(4.57)

q J ϭ ␮ 0 ␩ J2 ,

q

␯

ϭ

e

ij␲

␯ ij

,

qtϭ␳␬g•͑ ␣S“Scϩ␣␰“␰c͒.

(4.58)

23The molecular diffusivities are conceptually important in

boundary layers. It is clear from Eq. (4.38), for instance, that

␬tϭ0 on the core-mantle boundary because Vtϭ0 there. The

total diffusivity ␬totalϭ␬tϩ¯␬M is nonzero, however. As r increases within a thin thermal layer on the CMB, ␬total de-

creases

and

‫ץ‬r͗Sc͘

increases

in

compensation,

so

that

I

S r

ϭϪ␬total‫ץ‬r͗Sc͘ does not change. Since we are not particularly

interested in the boundary layer, we ignore it and apply Eq.

(4.38) throughout the core, right up to its boundaries, and we

similarly ignore the boundary layers of ͗␰c͘.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1102

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

As

usual,

␮ 0 Jϭ “ ÃB,

␳ FB ϭ JÃB,

and

e

ij

ϭ

1 2

(

ٌ

i

V

j

ϩٌjVi);

Eq.

(4.20)

gives

¯␴

S Ϫ

.

Also

⌸ϭPc /␳ϩUc .

Solutions to Eqs. (4.45)–(4.58) are subject to bound-

ary conditions; those governing B were dealt with in Sec.

III. The velocity should obey the no-slip conditions

VCMBϭ 0,

(4.59)

VSICϭ ⍀SICÃr,

(4.60)

where ⍀SIC is the angular velocity24 of the solid inner core. This is determined by solving the equation of mo-

tion for the inner core (considered as a rigid body) as it

moves under ⌫, the sum of the magnetic, viscous, topo-

graphic, and gravitational torques to which it is sub-

jected. The topographic torque arises because, in reality,

the inner core boundary does not have an axisymmetric

ﬁgure; it has ␾-dependent bumps, on which the hydro-

dynamic pressure creates a torque.

We observed in Sec. IV.B that the ICB is more generally a surface of constant U¯ rather than of constant r.

The most obvious departure from spherical symmetry is

the equatorial bulge of about 3 km created by the cen-

trifugal force. As for a spinning top, the bulge provides a

restoring force; the axis of rotation of the SIC does not
tip over progressively under the action of ⌫x and ⌫y , but precesses about Oz. Here, as in the simulations de-

scribed in Sec. VI, let us consider only the z components of ⌫ and ⍀SIC , together with the r␾ (and ␾r) stress components. The magnetic stresses on the ICB are so

large compared with the viscous and topographic

stresses that (in the absence of a gravitational torque)

they must, paradoxically, integrate almost to zero in the
computation of ⌫, and ⍀SIC must adjust itself to bring that about. This consequence of Lenz’s law led to a pre-

diction (see Sec. VI) that the SIC is rotating in a pro-

grade direction relative to the mantle at a few degrees

per year. This conclusion was later questioned by Buf-

fett (1997), who pointed out that mass in the mantle is

not distributed with perfect spherical symmetry, and the mantle inhomogeneities make U¯ in the core (weakly) ␾

dependent. Since the ICB is an equipotential,

24Equation (4.59) is a consequence of our choice of reference frame. The core exerts magnetic, topographic, viscous, and gravitational torques on the mantle that cause ⍀ to change, resulting in variations in the length of day (Sec. I.C). The magnetic torques arise because the mantle is not, as we have supposed, electrically insulating, although its conductivity is small. The main torque on the mantle by the ﬂuid outer core may arise indirectly, through the magnetic stresses exerted by the FOC on the inner core boundary and transmitted to the mantle by a gravitational torque acting as a catalyst; see Buffett (1996). Whatever its cause, it must be correctly incorporated into simulations such as those of Sec. VI. This is done by replacing the reference frame ﬁxed in the mantle by the reference frame in which the total angular momentum of the Earth is zero. The relative motion between these frames is very small. Instead of Eqs. (4.59) and (4.60), Kuang and Bloxham (1997) require that the FOC create no viscous stress on CMB and ICB; see Sec. VII.B.

␾-dependent bumps are created on it that are directly linked to the inhomogeneities in the mantle that create them. The gravitational interaction between the bumps and the imperfections tends to lock the rotation of the SIC to that of the mantle.
Bumps on the SIC are not avoided even when we return to our simplifying assumption of a spherically symmetric reference state. Although ¯rICB is then a function only of ¯t in the reference state, the ICB advances more rapidly wherever cold descending convection currents impinge on it, and more slowly where hot rising currents leave it, and bumps are created on the ICB (Glatzmaier and Roberts, 1998). Correspondingly, the radial components of I␰ and IS vary with ␪ and ␾ on the ICB. It is convenient to generalize the deﬁnition of the overbar to mean the horizontal average of even convective quantities and to write I␰ϭ¯I␰ϩI␰h, ISϭ¯ISϩISh, and

rICBϭ¯rICB͑¯t ͒ϩrIhCB͑ ␪,␾,tc ,¯t ͒,

(4.61)

where the superscript h stands for the horizontally varying part. We then ﬁnd that

¯I

␰ r

ϭ

␳

Ϫ

⌬¯␰

rថ ICB ,

at

rϭ¯rICB ,

(4.62)

¯I

S r

ϭ

␳

Ϫ

⌬¯S

rថ

ICB

,

at

rϭ¯rICB .

(4.63)

Suppose that there is no mass exchange across the coremantle boundary; then rCMBϭconst, if the slow but inevitable contraction of the Earth as it cools is ignored.
We then have

¯I

␰ r

ϭ

0,

at

rϭrCMB ,

(4.64)

¯I

S r

ϭ

͑¯I

q r

Ϫ

I

rad͒ /

T

,

at

rϭrCMB .

(4.65)

Conditions (4.62) and (4.64) ensure conservation of the light constituents, the ﬂuxes of ␰ at the SIC being precisely what is required to account for the volumetric source on the right-hand side of Eq. (4.47). Similarly, the difference between the total ﬂow of S from the CMB implied by Eq. (4.65) and the total ﬂow of S into the FOC across the ICB, as given by Eq. (4.63), is the integral over the FOC of the right-hand side of Eq. (4.49).
To ﬁnd rថ ICB , we must apply the condition T(¯rICB) ϭTl(P¯ ICB ,¯␰) of phase equilibrium, where Tl(P,␰) is the liquidus, assumed known. The proper implementation of this condition is lengthy. The detailed argument is given by Braginsky and Roberts (1995). Sufﬁce it to say here that, to a good geophysical approximation, it implies

Sថ ϭϪcP⌬*rថ ICB /¯rICB ,

(4.66)

wwhayeroen⌬t*hies

dimensionless and depends in a complicated latent heat (which is rather uncertain) and

the depression of the freezing point due to alloying

(which is even more uncertain). The error in the result-

lianrggev.alIuf eq,R⌬a*nϭd

0.05, is hard to estimate but is probably

QCMB

are

known,

we

can

obtain

¯I

S r

from

Eq. (4.65) and can then determine the mean rate rថ ICB of

advance of the inner core boundary from Eqs. (4.49),

(4.63), and (4.66).

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1103

The boundary conditions on the horizontally varying ﬁelds are similar:

I

r␰ h ϭ

␳

Ϫ

⌬¯␰

‫ץ‬

t

r

h ICB

,

at

rϭ¯rICB ,

(4.67)

I

S r

hϭ

␳

Ϫ⌬

¯S

‫ץ‬

t

r

h ICB

,

at

rϭ¯rICB ,

(4.68)

and, in analogy with Eq. (4.66),

‫ץ‬tShϭϪcP⌬*‫ץ‬trIhCB/¯rICB , at rϭ¯rICB . We also have

(4.69)

Ir␰hϭ0, at rϭrCMB ,

(4.70)

I

S r

hϭ

I

q r

,h

/

T

,

at

rϭrCMB ,

(4.71)

where

I

q r

,h

is

the

horizontally

varying

part

of

the

heat

ﬂux from core to mantle.

These arguments show that the sources of ␰c and Sc on the inner core boundary are proportional to one an-

other and to the local rate of advance of this boundary.

To determine that rate self-consistently, we need to

know qR and ICq MB(␪,␾), and we need to solve the dynamo problem.

F. The Boussinesq approximation

The theory of laboratory convection commonly assumes that the reference state is uniform. This is known as the Boussinesq approximation. It is adopted in most studies of core MHD and the geodynamo because it simpliﬁes the mathematics slightly, and that is its main purpose. In the context of core MHD, it introduces inaccuracies that are typically of order 20%. These are often viewed as tolerable in view of greater uncertainties in some of the other model parameters. More seriously, the approximation totally disregards the adiabatic gradient and adiabatic heat ﬂux. Also, the analysis of the energy and entropy balances requires extra care; see Sec. 8 of Braginsky and Roberts (1995).
Equations (4.45) and (4.46) become

“ • Vϭ 0,

(4.72)

‫ ץ‬t Vϩ V• “ Vϭ Ϫ “ ⌸ ϩ C gϪ 2 ⍀ÃVϩ JÃB/ ␳ 0

ϩ ␯ ٌ 2 V.

(4.73)

Usually thermodynamic linearization (Sec. IV.C) is carried out as a perturbation from T¯ rather than from ¯S, so that, in place of Eq. (4.53),

C

ϭ

Ϫ

␣

T

cϪ

␣

␰ T

␰

c

,

(4.74)

where

␣

␰ T

ϭ

Ϫ

␳

Ϫ

1(

‫ץ‬

␳

/

‫ץ‬

␰

)

P

,

T

is

the

isothermal

coefﬁcient

of volume expansion.

The Boussinesq approximation may be thought of as a

double limit in which g→ϱ and C→0 with gC ﬁnite and

nonzero. This requires that ␣S→0 and ␣␰→0, so that it

is

unnecessary

to

distinguish

between

␣

␰ T

and

␣␰

in

Eq.

(4.74). Then

Scϭ͑ cP /T¯ ͒Tc ,

(4.75)

and expressions (4.53) and (4.74) for C coincide by Eq.
(1.8). The equation governing Tc is obtained by multiplying Eq. (4.49) by the constant T¯ /cP .

V. RMHD

A. Orders of magnitude

RMHD stands for rotating magnetohydrodynamics, a subject so different from MHD that it deserves its own acronym.25 If ⍀ϭ0, we shall refer to RMHD as ‘‘classical MHD’’ and, if Bϭ0, we shall call it ‘‘classical rotating ﬂuids.’’ MHD is a subject that weds electrodynamics to hydrodynamics, the progeny resembling neither parent closely. Similarly, in RMHD, classical rotating ﬂuids and classical MHD are married, but the offspring are again surprisingly different. RMHD, as the subject stood in 1973, was reviewed by Acheson and Hide (1973).
Two dimensionless numbers arise in the classical theory of rotating ﬂuids, the Ekman number E and the Rossby number Ro:

Eϭ␯/⍀L¯ 2,

(5.1)

RoϭU¯ /⍀L¯ .

(5.2)

These measure the viscous force ¯␳␯ٌ2V and the inertial
force ¯␳V•“V against the Coriolis force 2¯␳⍀ÃV. In es-
timating E and Ro we shall abandon estimate (1.24) and, following common practice, use instead LϭrCMBϪ¯rICB Ϸ2260 km. We ﬁnd that RoϷ10Ϫ5 and EϷ5ϫ10Ϫ14 (for ␯ϭ␯M) or 10Ϫ7 (for ␯ϭ␯tϷ␩). The core is therefore a rapidly rotating ﬂuid, deﬁned as one in which

EӶ1,

(5.3)

R o Ӷ 1,

(5.4)

and in which therefore viscous and inertial forces are generally small compared with the Coriolis force. The ratio Ro/E is the kinetic Reynolds number Re and is large; see Eqs. (4.33) and (4.34). The geodynamo is selfexcited, and therefore (Sec. III.B)

RmտO͑1͒,

(5.5)

where RmϭU¯ L¯ /␩.

(5.6)

Thus, by inequality (5.4), the magnetic Ekman number ␩/⍀L¯ 2 is small, so that inequality (5.3) is automatically satisﬁed. The smallness of Ro means that from now on we shall discard V•“V in Eq. (4.73). Looking ahead to Sec. V.C, we shall recognize that the nonlinearity equilibrating the solutions is the Lorentz force and not the inertial force.

25It does not get one, however, except in this review. In plasma physics RMHD is an abbreviation for ‘‘reduced MHD.’’

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1104

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

There are several possible levels at which B could saturate, but in all of them the magnetic Rossby number is small:

R o B Ӷ 1,

(5.7)

where RoBϭVA /⍀L¯ ,

(5.8)

where VAϭB/(␮␳0)1/2 is a typical Alfve´ n velocity. Taking Bϭ0.002 T and the characteristic density ␳0 from Table I, we ﬁnd that VAϷ1.7 cm sϪ1, so that RoBϷ5 ϫ 10Ϫ 4 .

The saturation level of main interest here is that of

the strong-ﬁeld dynamo in which the Lorentz and Coriolis forces have similar magnitudes.26 Their ratio is ap-

proximately

V

2 A

/2⍀

U¯ L¯ ϭ

⌳

/

R

m

,

where

⌳ϭV A2 /2⍀␩ϭB2/2⍀␮0␳0␩

(5.9)

is the Elsasser number,27 which is independent of the

length scale L¯ . The strong-ﬁeld regime is therefore one

in which

⌳ϭO͑ 1 ͒, VAϷ͑ 2⍀␩ ͒1/2, Bϭ͑ 2⍀␮0␳0␩ ͒1/2, (5.10)

which gives Bϭ0.002 T for the core; the characteristic

current density is JϭB/␮0L¯ ϭ0.004 A mϪ2. Also, since

V

2 A

Ϸ

2

⍀

U¯ L¯ ,

the

Alfve´ n

number

U¯ /V¯ A

is

approximately

ͱRo, which is small. The magnetic energy density EB

therefore greatly exceeds EV, the kinetic energy density

(in the rotating reference frame): EB/EVϷRoϪ1. In some

of the simulations of Sec. VI, EB/EKϳ103. Energy is not

equipartioned.

Our plan now is to give little more than a thumbnail

sketch of some of the concepts and phenomena that are

signiﬁcant in Sec. VI. We start in Sec. V.B with the clas-

sical theory of rotating ﬂuids. This subject has an im-

mense literature, and it is obviously impossible to do

more than extract some of its ﬂavor through a few

simple examples. More details may be found in the clas-

sic text of Greenspan (1968), in Roberts and Soward

(1978), and in several books on the ﬂuid dynamics of

atmosphere and oceans, e.g., Gill (1982), Monin (1990),

and Pedlosky (1979). In Sec. V.C, we shall consider how

magnetic ﬁelds change the results of Sec. V.B. In Sec.

V.D, we discuss some matters of more direct relevance

to dynamo theory.

Simplicity is sought in this section, and we shall ini-

tially exploit the Boussinesq approximation (Sec. IV.F);

we shall usually exclude compositional buoyancy and

shall replace C by Ϫ␣Tc . We shall also generally ignore the solid inner core, so that the ﬂuid outer core ﬁlls r

ϽrCMB .

26A weak-ﬁeld dynamo is one in which the Hartmann number MϭVAL /ͱ(␯␩) is O(1), meaning that Lorentz and viscous forces are similar in magnitude. This would give BϷ0.3 nT for
the core, which is even less than the ﬁeld seen at the Earth’s
surface (Sec. I.B). 27As far as we are aware, Elsasser never wrote down his num-
ber (5.9), although he did identify the scale (5.10) of B; see
Elsasser (1946).

B. Classical rotating ﬂows

In this section we suppose that Bϵ0. Inequality (5.3)
suggests that solutions should be developed asymptotically, in the limit E→0. The plan therefore is to generate a mainstream28 expansion for solutions within the core, excluding boundary layers on the core-mantle boundary. The viscous force does not appear in the equation governing the leading-order mainstream solution. This lowers the differential order of the system, and we can require only that

Vrϭ0, on rϭrCMB .

(5.11)

The resulting solution will not obey the conditions (4.59)
on 1rÃV. It is necessary to develop a boundary layer solution that matches to the mainstream value of 1rÃV at its inner edge, while satisfying condition (4.59) in full
on the CMB. Initially we shall consider the mainstream
only and shall return to the boundary layer later.
Let us temporarily disregard buoyancy. The surviving inertial term ‫ץ‬tV in Eq. (4.73) is signiﬁcant only at high frequencies, as for inertial waves, which are determined
by condition (5.11) and

“ • Vϭ 0,

(5.12)

‫ ץ‬t Vϭ Ϫ “ ⌸ Ϫ 2 ⍀ÃV.

(5.13)

This is an eigenvalue problem yielding an inﬁnity of frequencies of order ⍀, though all less than 2⍀.
Consider next low frequencies. For steady motions (‫ץ‬tVϭ0) we ﬁnd, by operating on Eq. (5.13) by “Ã, that

2 ⍀• “ Vϭ 0.

(5.14)

This embodies the Proudman-Taylor theorem: The slow steady motion of a rotating inviscid ﬂuid is twodimensional with respect to the rotation axis. Taking ⍀ϭ⍀1z , we have

VϭV͑x,y ͒.

(5.15)

Because the boundary is axisymmetric, the only solution (5.15) that obeys condition (5.11) is

VϭV˜ G͑ s,t ͒1␾ ,

(5.16)

where (s,␾,z) are cylindrical coordinates; t has been included in Eq. (5.16) in view of later developments, but clearly the low-frequency inertial modes, though time dependent, have approximately the form (5.16). As in Sec. III.C, a tilde is used to denote axisymmetry. A ﬂow of type (5.16) is important in the angular momentum balance of the Earth (Sec. I.C) and is termed geostrophic. It is axisymmetric, zonal, and constant on geostrophic cylinders C(s) of constant radius s. Figure 8 shows a typical geostrophic cylinder in the case when a solid inner core is present. It also shows a very signiﬁcant imagi-

28
In asymptotic theory, the mainstream would usually be called ‘‘the outer solution’’ and the boundary layer ‘‘the inner solution’’ (see, for example, van Dyke, 1964), but such descriptions are confusing in the context of the core-mantle boundary.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1105

convective instability, but, if the heat sources are gradu-
ally increased, weak convection occurs as soon as the Rayleigh number Ra exceeds a critical value Rac . The Rayleigh number is a dimensionless measure of the thermal forcing, deﬁned here by

Raϭg␣␤L¯ 4/␯␬,

(5.21)

FIG. 8. Sketch of a geostrophic cylinder C(s) of radius s together with its northern and southern spherical end caps, N(s) and S(s). Also shown is the tangent cylinder C(rICB).

nary cylinder C(¯rICB), which touches the SIC on its equator and which is therefore called the tangent cylinder.
We next consider the action of the buoyancy force
through two examples, in the ﬁrst of which T˜ c is axisymmetric and is given. The steady inviscid ﬂow driven by
T˜ c obeys

“•V˜ ϭ0,

(5.17)

0ϭϪ“⌸˜ Ϫg␣T˜ cϪ2⍀ÃV˜ . Operating on Eq. (5.18) by “Ã, we see that

(5.18)

2⍀•“V˜ ϭϪ␣gÃ“T˜ c .

(5.19)

Since g(ϭϪg1r) is radial, we recover solution (5.16) unless T˜ c depends on latitude, in which case

V˜ ϭ͓V˜ T͑ s,z,t ͒ϩV˜ G͑ s,t ͔͒1␾ ,

(5.20)

where we have again included a t dependence for later convenience. The ﬂow V˜ T(s,z,t)1␾ is called the thermal wind. Its magnitude is O(gC˜ /⍀) and, if this is to be comparable with our assumed characteristic velocity U¯ ϭ5ϫ10Ϫ4 m sϪ1, the pole-equator difference in C˜ must be O(10Ϫ8). This is the origin of our estimate of C in Sec. IV.C. The corresponding T˜ c is of order 10Ϫ4 K, as in our estimate of ␦T in Sec. II.C. The state (5.20) may be subject to asymmetric baroclinic instabilities. These are studied in, for example, Pedlosky (1979). We shall touch on their magnetic analogs in Sec. V.C.
The second example is one of convective stability and requires us to restore viscosity in the mainstream, although we continue to assume that EӶ1. (For simplicity, we suppose that the Prandtl number ␯/␬ is not small.) Imagine that heat sources are distributed uniformly in the core. If these are weak, they create only a spherically symmetric temperature distribution Tc(r) that carries heat out of the core by thermal conduction. Although the associated density distribution is top heavy, the diffusion of heat and momentum prevents

where ␤ is a typical gradient of Tc(r). For RaϭRac , and also for Ra modestly in excess of Rac , convection takes the form of a ‘‘cartridge belt’’ of two-dimensional cells, often called Taylor cells, regularly spaced round the axis of rotation and drifting in longitude about that axis (Roberts, 1968; Jones et al., 2000).
Taylor cells are seen clearly in Fig. 9, which is taken from the dynamo model of Kageyama and Sato (1997). Adjacent cells rotate in opposite directions (about their axes) in a sequence of cyclonic and anticyclonic vortices, with vorticity respectively parallel and antiparallel to ⍀. The effect of these motions on the magnetic ﬁeld will be discussed in Sec. VI.B. The name ‘‘Taylor cell’’ is a useful reminder of the Proudman-Taylor theorem, which the ﬂow is trying to obey by being as 2D as possible. We see from our ﬁrst example that, if ␯ϭ0, small-amplitude motion must be geostrophic, since the thermal forcing Tc(r) is independent of latitude. But geostrophic motions have no radial components and cannot carry heat outwards. Convection can occur only if the viscous forces are large enough to break the rotational constraint of the theorem. Thus, although Rac would be O(1) if E were O(1), the critical Rayleigh number is large when EӶ1; in fact RacϭO(EϪ4/3), and the number of cells in the cartridge belt is of order EϪ1/3, i.e., the scale LЌ of the motions perpendicular to ⍀ is O(E1/3). Convective heat transport is mainly in the s direction, i.e., away from the rotation axis (Busse, 1970); this is signiﬁcant for the simulations of Sec. VI.
We now consider the boundary layer on the CMB, which is known as an Ekman layer. Dimensional analysis of Eq. (4.73) correctly indicates that its thickness ␦␯ is of order29 ͱ (␯/͉⍀r͉)ϰE1/2L¯ , where ⍀rϭ1r•⍀ϭ⍀ cos ␪. The Ekman layer is not passive; it controls the mainstream in the sense that the geostrophic ﬂow can only be determined through the Ekman layers. It does this through a process called Ekman pumping. To match the leading-order mainstream solution [which we temporarily denote by Vm(0s)] to condition (4.59), the Ekman layer has to pump ﬂuid radially, with a velocity of order E1/2Vm(0s) , into the mainstream. This provides a boundary
29There is obviously a complication at the equator of the CMB, where ⍀rϭ0 and ␦␯ϭϱ. The Ekman layer has a passive singularity at the equator; ␦␯ is not inﬁnite but is of order E2/5 (which is much greater than E1/2 for E→0); see Stewartson (1966). For an analysis of the Ekman layer, see Greenspan (1968).

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1106

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

FIG. 9. The cartridge belt of Taylor cells in the rotating, convective dynamo of Kageyama and Sato (1997). Cyclonic cells are colored blue and anticylonic pink; the equatorial plane is also shown. This ﬁgure illustrates how a zonal magnetic-ﬁeld line is distorted by the motion. The line is red where ϪV•(JÃB) is positive. This is where kinetic energy is transformed into magnetic energy; see Eq. (3.12). The white arrows denote V on the ﬁeld line. From Kageyama and Sato, 1997, with permission [Color].

condition that must be obeyed by the second term, Vm(1s) , in the expansion of the mainstream. Therefore Vm(1s) is O(E1/2Vm(0s)), rather than the smaller O(EVm(0s)) that Eq. (4.73) might have superﬁcially led one to expect a priori. A signiﬁcant application of this result is considered in Sec. V.D and in the next paragraph.
The Ekman layers adjust the rotation of a ﬂuid to that of its boundary, by a process called spin up which is surprisingly rapid. Suppose that mantle and core are spinning together and that the angular velocity ⍀ of the mantle is then increased slightly. Eventually the core must adjust to the new angular velocity, and it does this by creating an Ekman layer that sucks ﬂuid out of the mainstream, so drawing the preexisting z-directed vortex lines in the mainstream (as seen in the inertial frame) towards the rotation axis and increasing their

density to the new requirement set by the boundary.30 Since the rate of Ekman pumping is proportional to E1/2, the process of spin-up is essentially complete in a time ␶su of order EϪ1/2(2␲/⍀), i.e., only 3200 daysϷ 9 yr for ␯ϭ␯t.
Finally, when an SIC is present, the tangent cylinder divides the FOC into three regions, the northern interior
30We are here appealing to the Kelvin-Helmoltz theorem, an analog of the frozen-ﬂux theorem, according to which vortex lines in an incompressible inviscid ﬂuid are material curves that move with the ﬂuid. In the present case, the vortex lines are all parallel to Oz both initially and ﬁnally, and the Ekman suction essentially crowds them inwards as a whole. (The inertial waves described earlier may be visualized in the inertial frame as waves that travel along these vortex lines.)

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1107

of the tangent cylinder, the southern interior, and the exterior. The dynamics in these regions tend to be rather different, and they adjust to one another in complicated transitional regions surrounding the tangent cylinder. This was already apparent in Stewartson’s (1966) solution of the ‘‘Proudman problem,’’ which is to determine the ﬂow that is set up when the SIC spins at a slightly different rate from the mantle; complicated ‘‘Stewartson layers’’ surround the tangent cylinder.

C. Magnetic effects

We now consider how the results of Sec. V.B are
changed when a magnetic ﬁeld B is present. As before,
we ignore the inertial term V•“V and consider the limit E→0, focusing ﬁrst on the mainstream.
As is too well known to describe here, disturbances in classical MHD are transmitted by Alfve´ n waves. These
are dispersionless and travel, in a uniform ﬁeld B0 , with velocity VAϭϮB0 /ͱ(␮0␳0). Rotation has a profound effect when inequality (5.7) holds. The waves are then of
two distinct types: inertial waves (Sec. V.B) with frequencies of order ⍀, and MAC waves.31 This acronym
highlights the forces that (together with the pressure
gradient) are signiﬁcant: magnetic, Archimedean, and
Coriolis, though we have not yet included the Archimedean (buoyancy) force; the inertial force ‫ץ‬tV is conspicuous by its absence. MAC waves are dispersive; their fre-
quencies are of order V A2 /2⍀L¯ 2, which, by inequality (5.7), are much smaller than ⍀. When the inertial term ‫ץ‬tV is deleted from the equation of motion, a quasihydrostatic force balance remains: at leading order,

0ϭϪ“⌸Ϫg␣TcϪ2⍀ÃVϩJÃB/␳0 .

(5.22)

The inertial forces so essential for the Alfve´ n wave are insigniﬁcant for the MAC wave. It is this fact that makes RMHD so different from classical MHD. Time dependence enters only through the induction equation

‫ץ‬tBϭ“Ã͑ VÃB͒.

(5.23)

[The diffusive terms of Eqs. (4.52) and (4.73) have been ignored.] For small BϪB0 , the linearized forms of Eqs. (5.22) and (5.23) provide the dispersion relation for
MAC waves. Their time scales are of order

␶MACϭ2⍀L¯ 2/V A2 ϭ2⍀␮0␳0L¯ 2/B2,

(5.24)

or about 4000 years for Bϭ0.002 T. This is roughly the

time scale of the westward drift (Sec. I.C). Braginsky

(1964d) suggested that the observed geomagnetic secu-

lar variation is a manifestation of MAC waves, an idea

later developed by Hide (1966).

We now reconsider the same two examples of convec-

tive instability that were discussed in Sec. V.B, but now

in the presence of an axisymmetric zonal applied ﬁeld

31To call these fast and slow waves (as is often done) runs some risk of confusion with magnetoacoustic waves that are known by the same names.

B˜ 0ϭB˜ 0(s,z)1␾ . In the ﬁrst example, the ﬂow is driven by an axisymmetric temperature ﬁeld T˜ c(r,␪), which depends on the colatitude ␪. As in Sec. V.B, a zonal ﬂow is created, but it now has three parts: a thermal wind
obeying Eq. (5.19), a geostrophic ﬂow, and a magnetic
wind

V˜ BϭB˜ 02/2⍀␮0␳0s, so that Eq. (5.20) is generalized to

(5.25)

V˜ ϭ͓V˜ T͑ s,z,t ͒ϩV˜ B͑ s,z,t ͒ϩV˜ G͑ s,t ͔͒1␾ .

(5.26)

Braginsky (1967) presented a linear stability analysis

of state (5.26) that assumed all forms of diffusion to be

absent but that is otherwise general. If the thermal forc-

ing T˜ c is weak, the perturbations create MAC waves that travel longitudinally around the rotation axis. If it is

strong, the perturbations grow without limit, according

to the linear theory, but Braginsky (1964d) suggested

that nonlinearities would equilibrate these perturbations

at ﬁnite amplitude and that, since they are preferentially

asymmetric and therefore evade Cowling’s theorem (Sec. III.B), they would provide the nonzero ¯E␾ needed to maintain the geodynamo. Unstable waves tend to

travel westward (Acheson, 1972; see also Roberts and

Stewartson, 1975). The shortness of the time scale (5.24)

has encouraged speculations that the brevity of polarity

transitions (Sec. I.C) is a manifestation of more drastic

MAC instabilities. Many further studies of RMHD in-

stabilities have been completed since 1967, several of

which concern resistive instabilities such as tearing

modes. The subject has been recently reviewed by Fearn

(1998).

In RMHD, MAC waves take over the role of Alfve´ n

waves almost completely. We say ‘‘almost’’ because

there is one important class of motions for which this is

untrue: torsional waves. It is easy to see that the Coriolis

force associated with the geostrophic mode (5.16) is con-

servative and is therefore totally ineffective, since it can

be absorbed into the reduced pressure ⌸. Weaker forces, in particular ‫ץ‬tV˜ G , that would otherwise be ne-
glected become signiﬁcant. The torsional wave is essen-

tially an Alfve´ n wave in which the geostrophic cylinders C(s) turn about Oz and are linked to each other by the

component Bs(s,␾,z) of B that threads them together. The restoring force on cylinder C(s) depends on the in-

tegral

T( s )

of

B

2 s

/

␮

0

over

its

surface;

the

inertia

of

C( s )

is proportional to the integral m(s) of ␳. The torsional

wave therefore travels with the local wave speed Vtors(s)ϭͱ(T/m)ϷVA , where VA is a mean of the Al-
fve´ n velocity computed from B˜ s . This provides one of
the more rapid time scales of the macroscale ﬁelds:

␶torsϭrCMB /Vtors .

(5.27)

If we take Bsϭ2ϫ10Ϫ4 T, we obtain ␶torsϷ32 yr. Torsional waves are responsible for carrying the z compo-

nent of angular momentum across the core. As men-

tioned in Sec. I.C, several analyses of geophysical data

claim to have detected them (e.g., that of Zatman and

Bloxham, 1997).

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1108

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

Braginsky (1970) provided the complete theory of dis-
sipationless torsional waves. Roberts and Soward (1972)
showed that the waves are damped by the Ekman layers and, unless maintained, decay during a spin-up time ␶su . Although strictly speaking Eq. (5.22) applies only to the
ageostrophic ﬂow VϪV˜ G , it may also be applied to V˜ G if attention is focused on time scales long compared with ␶su . Otherwise part of the inertial acceleration, namely, (‫ץ‬V˜ G /‫ץ‬t)1␾ , must be restored to Eq. (5.22). As a beneﬁcial byproduct, the ‘‘stiffness’’ of numerical simula-
tions is reduced by including (‫ץ‬V˜ G /‫ץ‬t)1␾ (Jault, 1995). Consider now the second example, the convective in-
stability of a rotating sphere containing heat sources, but
now permeated by a magnetic ﬁeld. This ﬁeld can be
more effective than viscosity in breaking the rotational
constraint. Considering the critical Rayleigh number Rac as a function of the ﬁeld strength B, it is found that Rac is smallest in the strong-ﬁeld range, ⌳ϭO(1). The minimum Rac is O(EϪ1), which is much smaller for E Ӷ1 than the O(EϪ4/3) critical value for Bϭ0. We see that the magnetic ﬁeld facilitates convection.32 This is not totally surprising. In the absence of ﬁeld, the scale LЌ of the convective motions in directions perpendicular to ⍀
has to be small ͓LЌϭO(E1/3L¯ )͔ to break the rotational constraint, and such motions are energetically expensive. When ⌳ϭO(1), the scale of the cells is much
larger ͓LЌϭO(L¯ )͔, and these are less costly. (The increased scale of convection can be inferred from Fig. 9.
Taylor cells, though evident, are not as numerous as
they would be if B were zero.) The optimal case, Rac ϭO(EϪ1), may be restated as Rac*ϭO(1), where

Ra*ϭg␣␤L¯ 2/2⍀␬

(5.28)

is a Rayleigh number that is independent of ␯, as is the optimizing ﬁeld strength B given by Eq. (5.10). Viscosity plays no role in the mainstream solution. Solutions for RaϾRac have been derived by Walker and Barenghi (1999).
Although it might seem that marginal stability calculations (Rmϭ0) of the type just described have no direct bearing on the geodynamo ͓RmտO(1)͔, they have proved to be a reliable guide in predicting the existence of strong-ﬁeld convective dynamos such as that of Jones and Roberts (2000). A few words of caution are appropriate here. It is often said that the North-seeking property of the magnetic compass needle establishes that the Coriolis force dominates the dynamics of the core. This is an overstatement. More precisely, the Coriolis force is the only force acting that has a preferred direction. The Lorentz force is equally dominant, and, in acting to counter Coriolis forces and break the rotational constraint, it too acquires the same preferred direction. This

32It is sometimes speculated that convective dynamos have a general thermodynamic reason for their existence: a body of convecting ﬂuid generates a dynamo ﬁeld so that it can cool as rapidly as possible.

FIG. 10. Illustration of how Lorentz forces act to oppose Coriolis forces. Fluid velocity, Coriolis force, and Lorentz force are plotted in the equatorial plane. From Sakuraba and Kono, 1999, with permission.
is well illustrated by Fig. 10, which is taken from the dynamo simulation of Sakuraba and Kono (1999).
Consider next the boundary layers. In classical MHD at high Hartmann number MϭVAL¯ /(␯␩)1/2, asymptotic methods apply in which both ␯ and ␩ are set to zero in determining the leading-order mainstream, and this is matched to thin ‘‘Hartmann layers’’ on the boundaries. In dynamo theory, however, we are almost always concerned with a mainstream in which RmϭO(1) and in which therefore ohmic diffusion is not conﬁned to boundary layers. (An exception arose in Sec. III.D.) The Ekman-Hartmann layers relevant to core dynamics therefore have the character of a magnetically modiﬁed Ekman layer rather than of a Hartmann layer; e.g., see Loper (1970).
Finally we temporarily reintroduce the SIC. The effect of a magnetic ﬁeld on Stewartson (1957) layers surrounding the tangent cylinder has been analyzed by Hollerbach (1994, 1996b), who found that Bs tends to thicken them. Nevertheless, regions of large shear near

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1109

the tangent cylinder appear in the simulations of Sec. VI and seem to be signiﬁcant in the dynamo mechanism. The MHD of the three regions of the ﬂuid outer core separated by the tangent cylinder can be very different (Sec. V.B). The SIC also has another signiﬁcant effect on core MHD: a large-scale magnetic ﬁeld threading the SIC cannot change drastically in a time shorter than the electromagnetic time constant of the SIC, which is about 2500 years. This tends to prevent large-scale ﬁelds in the FOC from changing more rapidly than this (Hollerbach and Jones, 1993).

D. The Taylor state and model z

We have argued in Sec. V.C that, over time scales long compared with ␶su , inertial forces are negligible and that, at least to leading order, viscosity can also be ignored. This was the basis of our second convection example, which we emphasized was, in its optimal marginal state, independent of ␯ in all respects. The question arises therefore whether anelastic RMHD dynamos can be constructed which, except in boundary layers, satisfy

“ • ͑ ␳ V͒ ϭ 0,

(5.29)

0ϭ Ϫ “ ⌸ ϩ C gϪ 2 ⍀ÃVϩ JÃB/ ␳ ,

(5.30)

or the equivalent for the uniform ␳ (Boussinesq) models. By integrating the ␾ component of Eq. (5.30) over the curved surface of the geostrophic cylinder C(s) and ap-
plying Eq. (5.29), we obtain

͵ ͵ C͑s͒͑JÃB͒␾dSϭ2⍀ C͑s͒␳V•dS

͵ ϭϪ2⍀

Nϩ

␳
S

V•

dS,

(5.31)

where N(s) and S(s) are the northern and southern spherical caps at the ends of the cylinder C(s); see Fig. 8. By applying the boundary condition (5.11), Taylor (1963) obtained

͵C͑ s ͒ ͑ JÃB͒ ␾ d S ϭ 0,

(5.32)

which is known as Taylor’s condition or Taylor’s constraint. Unless this is satisﬁed, Eq. (5.30) has no solution obeying condition (5.11); if it is satisﬁed, there is an inﬁnity of solutions.
An alternative approach applies “Ã to Eq. (5.30) to obtain

2⍀‫ץ‬zVϭϪ“Ã͑ CgϩJÃB/␳ ͒.

(5.33)

Regarding CgϩJÃB/␳ as known, we can determine Vs and Vr uniquely by two simple integrations along every line segment parallel to Oz; condition (5.11), applied at
the ends of the segment where it meets the CMB, deter-
mines the two arbitrary ‘‘constants’’ of integration.
Equation (5.29) is, in component form,

␳‫␾ץ‬V␾ϭϪ‫ץ‬s͑s␳Vs͒Ϫs‫ץ‬z͑␳Vz͒.

(5.34)

This gives V␾Ј uniquely from VsЈ and VzЈ , so that the asymmetric ﬂow is completely determined. For the axi-
symmetric ﬂow, there is a difﬁculty: when the V˜ s and V˜ z determined above are substituted into the right-hand
side of Eq. (5.34), it is generally found to be nonzero. It
is only when condition (5.32) is obeyed that the right-
hand side of Eq. (5.34) is zero. When this happens, we
can determine V˜ ␾ by integrating the ␾ component of Eq. (5.33) along the line segments, but, since there are
no restrictions on V˜ ␾ at the ends of the segment, there is an arbitrary ‘‘constant’’ of integration, corresponding to
an arbitrary geostrophic ﬂow V˜ G(s,t)1␾ . This ﬂow is dynamically innocuous, since its Coriolis force can be absorbed into “⌸ as in Sec. V.C.
Condition (5.32) is related to the angular momentum balance about Oz. If the ‫ץ‬tV were restored to Eq. (5.30), one might hope that the system would respond to a failure of condition (5.32) by generating torsional waves that, if sufﬁciently damped, would adjust V˜ G until, in a comparatively short time of order ␶su , the condition would be satisﬁed. The magnetic ﬁeld subsequently would evolve slowly in obedience to

“ • Bϭ 0,

(5.35)

‫ ץ‬t Bϭ “ Ã͑ VÃB͒ Ϫ “ Ã͑ ␩ “ ÃB͒ ,

(5.36)

and (JÃB)␾ would evolve slowly also. Condition (5.32) will be gradually violated unless some action is taken. From Eq. (5.36), Taylor (1963) developed an evolution
equation for V˜ G which ensured that condition (5.32), if satisﬁed initially, would be satisﬁed for all t. This deter-
mined the evolution of V˜ G uniquely, so completing the determination of V.
Although Taylor’s prescription provides a clear-cut program for the construction of a strong-ﬁeld dynamo, no successful three-dimensional dynamo has yet been found by following it; see Fearn and Proctor (1987). It was shown by Roberts and Stewartson (1975) that one should not take it for granted that the Taylor prescription will work; see also Jones and Roberts (2000). The Taylor method may run into difﬁculties through treating the Ekman layers as passive layers that can always be constructed (after the mainstream has been found) in such a way that conditions (4.59) are satisﬁed in full.
One of the tacit assumptions behind the Taylor idea is that there exists a solution of the (viscous) RMHD equations that, as E→0, has an E-independent mainstream form, as in the second convection example of Sec. V.B. Braginsky (1975, 1978) questioned this assumption and proposed an alternative, in which the mainstream solution has a geostrophic part proportional to an inverse power of E as E→0. He called this model z because the geostrophic ﬂow increasingly makes the lines of force of the meridional ﬁeld Bp almost parallel to the z axis. The idea behind model z may be clariﬁed by regarding the large geostrophic ﬂow, of order EϪ1/2(JÃB/2⍀␳)ϭEϪ1/2⌳(␩/L¯ ), as the leading-order mainstream ﬂow Vm(0s) . Ekman pumping from N(s) and

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1110

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

S(s) created by this ﬂow is of order E1/2 times smaller (Sec. V.C), i.e., it is O(⌳␩/L¯ ) and is a part of Vm(1s) . Condition (5.32) no longer follows; the integral on the right-hand side of Eq. (5.31) is now in general nonzero and is of the same order as the left-hand side, to which it can be made equal, so depriving Taylor’s condition of any signiﬁcance. The remainder of Vm(1s) generally contains asymmetric parts, such as the ﬁeld-generating MAC waves envisaged in Sec. V.C, and these evade Cowling’s theorem.
It should by now be clear that the difﬁculty in solving Eqs. (5.29) and (5.30) arises from the omission of the viscous term F␯ and is conﬁned to the axisymmetric part V˜ of V. This has led to a number of investigations of axisymmetric MHD dynamos in which either C in Eq. (5.30) or V˜ T in Eq. (5.26) is speciﬁed as the energy source. The aim of these models is to understand better the role of the Taylor state, model z, and indeed other types of force balance. Interest in the electrodynamics is subsidiary, but nevertheless Cowling’s theorem must be evaded, and this is most often done in the simplest possible way, by including an ␣ effect, thus leading to Eqs. (3.22) and (3.23) with E˜ ϭ␣B˜ . See Hollerbach (1996a) for a recent review of this area.
VI. MHD DYNAMO SIMULATIONS
A. The development of models
Here we focus only on dynamo simulations that simultaneously solve for the thermodynamic variables, the ﬂuid ﬂow, and the magnetic ﬁeld in three dimensions (3D) with full time dependence and feedbacks. The resulting magnetic ﬁeld is maintained by rotating convection for several magnetic diffusion times ␶␩ . This excludes many studies of mean-ﬁeld and kinematic dynamos, convection in the presence of an externally applied ﬁeld, and solutions in a drifting frame of reference, though all of these have materially aided the development of dynamo theory. The comparatively simple model of Childress and Soward (1972) should, however, be mentioned, since very considerable progress can be made analytically using this model; see also Soward (1974).
The Childress-Soward model is an MHD dynamo driven by convective motions in a plane layer rotating about the vertical, heated from below and cooled from above (the rotating Be´ nard layer). Its success rests on the fact (Sec. V.B) that, provided ⌳ϽO(1), the scale LЌ of the motions in the (horizontal) direction perpendicular to ⍀ is small compared with the depth L¯ of the layer.33 Such microscale motions induce a macroscale ﬁeld in a way similar to that described in Sec. III.C for
33In spherical convection, the cartridge belt of Taylor cells (Sec. V.B) also deﬁnes a two-scale motion, suggesting that the Childress-Soward idea can again be usefully employed to ﬁnd an MHD dynamo by mainly analytic means; see Busse (1975).

turbulent ﬂow, though their ␣ effect is 2D rather than 3D, and ␣ is a function of z. As the Rayleigh number (5.21) is increased beyond Rac , a bifurcation is reached at which kinematic dynamo action occurs; a further increase in Ra results in ﬁnite B and a concomitant increase in the horizontal scale LЌ of the motions (Sec. V.C), though ⌳ remains small and the dynamo is of the weak ﬁeld variety. Eventually, as Ra is increased further, an asymptote is reached where B→ϱ and LЌ→L¯ . The result B→ϱ is a symptom of the violation of the assumption LЌӶL¯ on which the analysis rests; in reality, B does not become unbounded, but the asymptote shows very convincingly that a strong-ﬁeld regime (⌳տ1) exists in which LЌϭO(L¯ ); see Fautrelle and Childress (1982). The analytic method cannot follow the solution into this regime for which, as for all strong-ﬁeld dynamos, numerical computation is required.
Numerical integrations of planar models were carried out by Meneguzzi and Pouquet (1989), Brandenburg et al. (1990), St. Pierre (1993), and Jones and Roberts (2000). That of St. Pierre (1993) was the ﬁrst to reach a parameter regime in which the magnetic energy EB exceeds the kinetic energy EV of the convection that maintains it: EB/EVϷ10 (Sec. V.A). It also demonstrated how the two major forces in the problem, the Coriolis and Lorentz forces, nearly balanced each other locally (see Fig. 10 above). Planar models avoid the complication34 of spherical geometry and are convenient for testing the effects of various physical parametrizations (Vel´ımsky´ and Matyska, 2000), but obviously cannot represent global modes.
Gilman and Miller pioneered spherical MHD simulations in the early 1980s with the development of a global solar dynamo model (Gilman and Miller, 1981; Gilman, 1983). This model employed the Boussinesq approximation (Sec. IV.F), which was singularly inappropriate for the interior of the Sun. Glatzmaier (1984, 1985a, 1985b) developed an alternative anelastic model (Sec. IV.C), which allowed for large variations of density with depth and which, by ﬁltering out sound waves, allowed much larger numerical time steps than could have been taken had compressibility been fully included. These early solar dynamo models produced cycles of magnetic reversals in some ways similar to the migration of large-scale ﬁeld in the solar cycle.
The ﬁrst Earth-like magnetic ﬁeld was generated by Glatzmaier and Roberts (1995a, 1995b) using a 3D global model designed to simulate the core MHD. The original version included only thermal buoyancy and used the Boussinesq approximation. The current version (Glatzmaier and Roberts, 1996a, 1996b, 1997, 1998;
34Spherical models make use of a spectral transform in ␪ that is not fast, but planar models can make use of the fast Fourier transform in all three coordinate directions. It is therefore practical to compute planar MHD models on a workstation or a PC, truncating the Fourier modes after, for example, 64 terms in all three coordinate directions (Jones and Roberts, 2000).

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1111

Glatzmaier et al., 1999) accounts for both thermal and
compositional buoyancy and uses the anelastic approxi-
mation, with a reference state ﬁtted to the Earth and
convection driven by a prescribed, Earth-like, heat ﬂux
through the CMB, precisely as explained in Sec. IV above, the momentum ﬂux ␳V playing a role similar to the velocity V in the Boussinesq approximation; cf. Eqs.
(4.45) and (4.72). The magnetic ﬁeld outside the core in
this model had an intensity, structure, and time depen-
dence similar to the geomagnetic ﬁeld (Sec. I.C). A
Boussinesq model was later developed by Kuang and
Bloxham (1997, 1999) that used different velocity and thermal boundary conditions and operated in a somewhat different parameter regime but one that also produced an Earth-like magnetic ﬁeld outside the core. Both of these models were able to operate with relatively small diffusion coefﬁcients, which increased the effects of buoyancy and Coriolis forces; moreover, in rough agreement with estimates for the Earth, EB/EV Ϸ103; see Sec. V.A.
The Glatzmaier-Roberts model is an improved version of Glatzmaier’s original solar model; it computes the magnetic ﬁeld within a ﬁnitely conducting solid inner core, the importance of which was demonstrated with a 2D mean-ﬁeld model by Hollerbach and Jones (1993), and it treats the Coriolis force implicitly in the time integration, which makes it possible to operate in more extreme (less diffusive) parameter regimes. It uses the spectral transform method, the spatial resolution chosen (i.e., the number of modes retained in the spectral expansions) being dictated by the time span to be simulated. High resolution is affordable only for time spans of a few thousand years (Roberts and Glatzmaier, 2000), but more than a million years have been simulated at low resolution, using a numerical time step of about 15 days (Glatzmaier et al., 1999). In this review we present results from a medium-resolution case: 49 radial (Chebyshev) levels in the ﬂuid outer core (plus 17 in the solid inner core), 144 latitudinal levels, and 144 longitudinal levels. This corresponds to a rhomboidal truncation of spherical harmonics up to order 47 and degree 95. We set the core size, rotation rate, density proﬁle, CMB heat ﬂux, and magnetic diffusivity (␩ϭ2 m2 sϪ1) to Earth-like values and use viscous, thermal, and compositional eddy diffusivities to account for mixing by the unresolved turbulence (see Sec. IV.D). The diffusivity ␬t of heat and composition is taken as ␩ ; the viscous eddy diffusivity ␯t is 750 times greater. In addition, we use a ‘‘hyperdiffusivity’’ that increases these values (according to the square of the spherical harmonic degree) such that the diffusivity experienced by the degree 95 modes is ten times greater (ﬁve times greater for the magnetic diffusivity) than that experienced by the low-degree (large-scale) modes. This means that the Ekman number E is 5.4ϫ10Ϫ6 for the large scales and 5.4ϫ10Ϫ5 for the small scales. Here, and for the remainder of this review, we follow the usual practice of taking L¯ to be the radial depth of the FOC:

L¯ ϭrCMBϪ¯rICBϷ2.26ϫ106 m.

(6.1)

In order to include torsional waves, the axisymmetric

inertial terms are retained, but the asymmetric accelera-

tions are discarded.

Simulations such as these are speciﬁc to the Earth.

They are expensive and do not provide a complete un-

derstanding of the fundamentals of convective dynamos

in rotating spherical shells. For example, they are

strongly driven and create ﬁelds on many time and

length scales; they do not provide information about the

bifurcation structure of the solutions as the sources of

buoyancy are gradually strengthened from the marginal

state. They also sample parameter space sparsely. Many other models have been studied recently,35 especially

ones in which diffusion plays a more prominent role and

in which, therefore, the ﬂow and ﬁeld structures are

more dominated by larger scales. Several assume that

the SIC is electrically insulating. Apart from the

Kageyama-Sato (1995, 1997) model, which is fully com-

pressible, all these simulations employ the Boussinesq

approximation. The Kageyama-Sato model also differs

from the others by employing a perfect gas equation of

state. It is therefore similar to the early solar dynamo

models, e.g., those of Gilman and Miller (1981), Gilman

(1983), and Glatzmaier (1984, 1985a, 1985b).

In order to reach the more geophysically realistic pa-

rameter regimes of the Glatzmaier-Roberts and Kuang-

Bloxham models without their huge computational cost,

a

‘‘2

1 2

D’’

Boussinesq

model

was

developed

by

Jones

et al. (1995), Sarson et al. (1997, 1998), Sarson and Jones

(1999), Morrison and Fearn (2000), and Sarson (2000).

These do what in the solar physics community used to

be called modal simulations, in which the spatial resolution is adequate in the r and ␪ directions, but which

retain only the longitudinal wave numbers 0 and 2 (or

wave numbers 0, 2, 4, 6, and 8 in the case of Sarson,

2000). Although the longitudinal structure is coarse, the

axisymmetric parts of the ﬁelds appear to be quite simi-

lar to those of the fully 3D simulations of Glatzmaier

and Roberts.

B. Some results
Magnetic ﬁeld is generated by twisting and shearing of existing magnetic ﬁeld into new magnetic ﬁeld, and by magnetic diffusion that causes ﬁeld lines to break and reconnect into new topologies, as in the turbulent ␣ effect described in Sec. III.C. Taylor cells create a largescale ␣ effect of the type discussed in Sec. III.D. The dynamo simulation of Kageyama and Sato (1997) illustrates this nicely; see Fig. 9. Taylor cells of rotating ﬂuid outside the tangent cylinder have local angular velocity (vorticity) either parallel (cyclonic) or antiparallel (anticyclonic) to ⍀. In addition, due to the spherical bound-

35See Kageyama and Sato (1995, 1997), Kida et al. (1997), Busse et al. (1998), Christensen et al. (1998, 1999), Kitauchi and Kida (1998), Grote et al. (1999, 2000), Kageyama et al. (1999), Olson et al. (1999), Sakuraba and Kono (1999), and Kutzner and Christensen (2000).

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1112

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

aries, ﬂuid drifts along these cells, so that the motions have a left-handed helical sense in the northern hemisphere, and a right-handed sense in the southern hemisphere. They twist the zonal magnetic-ﬁeld lines, whether directed eastward or westward, into righthanded helices in the northern hemisphere and lefthanded helices in the southern hemisphere. The average over longitude ␾ is a meridional ﬁeld Bp (see also Glatzmaier, 1985a, and Olson et al., 1999). This completes
step 1 of the two-stage process maintaining B˜ described in Sec. III.C.
In addition, zonal ﬁeld has to be generated from meridional ﬁeld by step 2 of the two-stage process of Sec. III.C. Kageyama and Sato (1997) ﬁnd that in their simulation the zonal ﬁeld is mainly created from the meridional ﬁeld by the zonal shear outside the tangent cylinder. In the terminology of Sec. III.C, B˜ is maintained by an ␣␻-dynamo process. Olson et al. (1999), however, demonstrate that weakly driven convective dynamos, which are dominated by Taylor cells outside the tangent cylinder, usually sustain their B˜ by an ␣2-dynamo mechanism. Their more strongly driven (less diffusive) dynamos maintain B˜ by an ␣␻ process, step 2 being dominated by a thermal wind shear inside the tangent cylinder. Indications of this are seen in Fig. 11 from Olson et al. (1999), where (in the northern hemisphere) a left-handed helical upﬂow in the polar region horizontally diverges below the core-mantle boundary. This correlates with enhanced convective heat ﬂux and weak radial magnetic ﬁeld at the poles. Reducing viscous, thermal, compositional, and magnetic diffusion also strengthens the dynamo process and increases the kinetic and magnetic energies. The ratio of these energies depends on the relative values of these diffusivities (Busse et al., 1998).
Both the Glatzmaier-Roberts and the KuangBloxham models operate with relatively small diffusion coefﬁcients. The convective velocity is very time dependent but typically has a maximum of a few mm sϪ1. Convection is so effective that the maximum variation in T (on a sphere of constant radius) is only about 10Ϫ3 K. The maximum magnetic-ﬁeld intensity (in the deep interior of the core) can be as large as 50 mT and the magnetic energy (integrated throughout the core) is typically more than 2000 times greater than the kinetic energy of the convection that maintains it: EB/EVϭ2ϫ103; see Sec. V.A.
Large magnetic ﬁelds with large spatial gradients produce strong Lorentz forces, which need to be approximately balanced by the other major force in this problem, the Coriolis force. This was illustrated in Fig. 10 above, where the Coriolis and Lorentz forces were seen to nearly balance in one of the downwellings.
We now show results from the Glatzmaier-Roberts integrations. The simulated ﬂuid ﬂow outside the tangent cylinder is a combination of axisymmetric zonal and meridional circulations and nonaxisymmetric Taylor cells; see the snapshot of Fig. 12, where the equatorial projection of the ﬂow is shown. The Taylor columns are

smaller in scale, less well deﬁned, and more time dependent than they are in more diffusive solutions. Although the patterns of the thermodynamic perturbations, ﬂuid ﬂow, and magnetic ﬁeld are continually changing in all directions, one can identify a westward phase propagation. Even on the inner core boundary, the patterns of heat ﬂux, composition ﬂux, and inner core growth rate tend to propagate westward relative to the rotating frame of reference (Glatzmaier and Roberts, 1998). The phase velocity depends on location and time but is typically about 0.1° yrϪ1, similar to the observed westward drift of the Earth’s ﬁeld (Sec. I.C).
The axisymmetric zonal ﬂow outside the tangent cylinder is rather constant (Fig. 13), although there is a weak minimum in the equatorial region. Much greater variations exist inside the tangent cylinder where, relative to the mean rotation rate of the model Earth, the zonal ﬂow is westward near the mantle and eastward near the solid inner core. This is principally a thermal wind, whose existence may be understood in the following way. The efﬁcient Taylor column convection outside the tangent cylinder maintains that region at a fairly uniform temperature and composition compared with the less efﬁcient convective ﬂow inside the tangent cylinder (Fig. 13), where thermal and compositional buoyancy causes ﬂuid near the polar axis to ﬂow outward along that axis. Mass conservation requires the ﬂuid to return near the tangent cylinder, to ﬂow toward the rotation axis near the ICB and away from it near the CMB. Since the angular momentum is approximately conserved, ﬂuid inside the tangent cylinder spins up (moves eastward) near the ICB, and spins down (moves westward) near the CMB (Fig. 13). Although the ﬂuid near the ICB ﬂows eastward, the patterns of thermodynamic perturbations, ﬂuid velocity, and magnetic ﬁeld there all drift westward; see above.
The radial magnetic ﬁeld threads the ﬂuid outer core and solid inner core together and provides a potentially strong magnetic torque between them, as described in Sec. IV.E. To balance (on average) the small viscous torque on the ICB, the SIC must rotate at roughly the same average rate of the ﬂuid just above it; see Sec. IV.E. This mechanism, similar to an induction motor, is what maintains the super-rotation of the SIC predicted by Glatzmaier and Roberts (1995a, 1996b). This prediction motivated seismologists to look for, and ﬁnd, evidence of inner core super-rotation (Song and Richards, 1996; Su et al., 1996). The initial prediction and those original seismic estimates placed the rate at roughly 2° yrϪ1. More recent seismic analyses (Creager, 1997; Souriau, 1998) suggest a rate of order 0.1° yrϪ1 or less. The maximum, according to the Glatzmaier-Roberts simulation described here, is about 1° yrϪ1. Simulations (Buffett and Glatzmaier, 2000) that include a parametrized gravitational coupling (Sec. IV.E) between density variations in the mantle and an estimated topography of the inner core boundary give average superrotation rates as small as 0.02° yrϪ1.
Since the SIC rotates as a solid body while the zonal ﬂow just above it does not, a local shear exists in the

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1113

FIG. 11. Heat ﬂux from the core, radial components of ﬂuid velocity and magnetic ﬁeld, and helicity plotted on a surface below the core-mantle boundary. Blue means outwardly directed heat ﬂux and red inwardly directed; reds are positive, blues are negative for Vr and Br ; red helicity is right handed, blue is left handed. From Olson et al., 1999, with permission [Color].

zonal ﬂow near the ICB. Shear ﬂow also exists across the tangent cylinder due to the transition between the different styles of convection inside and outside this cylinder. These strong shear ﬂows, not seen in the more diffusive dynamos, are responsible for much of the zonal ﬁeld generation in the Glatzmaier-Roberts simulations, which appears to maintain its B˜ by an ␣␻ mechanism.
The Glatzmaier-Roberts geodynamo simulations were the ﬁrst to demonstrate spontaneous magnetic dipole re-

versals. The initial reversal (Glatzmaier and Roberts, 1995b) occurred about 3.5ϫ104 yr into the original simulation and took a little more than 103 yr to complete.
Before and after the reversal the dipole polarities were
opposite inside and outside (roughly) the tangent cylin-
der. (A somewhat similar dual-polarity conﬁguration is
seen in the modal calculations of Sarson, 2000.) We con-
tinued our simulation at this point using the anelastic model. After about 104 yr, the ﬁeld outside the tangent

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1114

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

FIG. 12. Fluid velocity plotted in the equatorial plane for a snapshot from the Glatzmaier-Roberts simulation. The lengths of the arrows are proportional to the ﬂow speed.

cylinder reversed again, leaving a single dipole polarity throughout the core (Glatzmaier and Roberts, 1996a). This new polarity conﬁguration is apparently more stable, since, although the ﬁeld was time dependent, it maintained this new polarity state for the next 2.3 ϫ105 yr. It then quickly reversed in less than 103 yr (although it took a further few thousand years for the axial

dipole to become dominant again). It remained in this polarity state for the next 1.7ϫ105 yr, quickly reversed
back, and has been in this latest state for the remaining 105 yr of the simulation (Glatzmaier et al., 1999). This 5.5ϫ105 yr simulation (over 27 magnetic dipole diffusion times) took 16 million numerical time steps; 3ϫ105 yr of
this are illustrated in Fig. 14, where it can be seen that,

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

FIG. 13. Longitudinally averaged density perturbation, meridional circulation, and angular velocity (relative to the rotating frame) for a snapshot from the Glatzmaier-Roberts simulation. Solid contours are positive; broken are negative so that, for the angular velocity, solid contours represent superrotation. In the center plot, the lengths of the arrows are proportional to the meridional ﬂow speed.

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1115

FIG. 14. 300 000 years in the middle of a 550 000-year Glatzmaier-Roberts simulation showing the evolution of the magnetic dipole (outside the core) in terms of its South-pole trajectory (in an equal-area projection with the North geographic pole at the top), its pole latitude, and its dipole moment. The South magnetic pole of the dipole is plotted once per 100 years and the dipole moment is in units of 1022 A m2. From Glatzmaier et al., 1999, with permission.
as in the palaeomagnetic reversal records, the dipole moment decreases signiﬁcantly during reversals. The frequency of reversals and their durations in this simulation also compare well with those of the Earth. It is tempting to regard a ﬁeld reversal as just another ﬂuctuation (albeit a large one) of a system driven so strongly that it varies stochastically on all scales. The second panel of Fig. 14 makes this view hard to defend. The level of secular variation both before and after the reversal is extremely low, but is quite typical of the model at all times. A reversal occurs like a ‘‘bolt from the blue.’’
During the ﬁrst reversal, the ﬁeld reversed inside the tangent cylinder about a thousand years before it did so outside, but during the next two reversals the opposite sequence occurred. That is, when viewed from the Earth’s surface, the ﬁeld would appear to have completed the reversal a thousand years or so before it actually reverses deep within the core (Glatzmaier et al.,

1999). The four snapshots in Fig. 15 illustrate how the

radial component of the ﬁeld at the surface and the axi-

symmetric part of the ﬁeld throughout the core change

during a typical reversal.

Occasional, spontaneous dipole reversals have also

occurred in recent modal calculations (Sarson and

Jones, 1999; Sarson, 2000) and in other 3D simulations

(Kagayama et al., 1999); they differ in details, as they

surely have in all past geomagnetic reversals. The rever-

sals seen in the modal calculations are associated with

ﬂuctuations in the (axisymmetric) meridional circula-

tion, possibly due to buoyancy surges originating near

the inner core boundary. It is not clear, however, what

triggers what, since buoyancy surges continually occur in

these simulations, and the ﬂuctuations in the meridional

circulation sometimes appear after the magnetic dipole

begins to decrease. We also see changes in the structure

of the meridional circulation Vp during reversals, but these usually occur after the dipole begins to lose inten-

sity.

Some of the ﬂuctuations continually occurring in the

Glatzmaier-Roberts simulations are strong enough to

produce a ﬁeld that is locally of the opposite polarity.

Only once in many attempts does one survive long

enough to challenge the original polarity seriously. In

some of these cases the ﬁeld reverses outside the tangent

cylinder but not inside it. Then, instead of a complete

reversal, the ﬁeld outside the tangent cylinder quickly

reverses back to its original polarity (Glatzmaier et al.,

1999). This may provide an explanation of cryptochrons

(Sec. I.C). The modal calculations of Sarson (2000) also

show a tendency for the surface ﬁeld to reverse more

frequently than the entire ﬁeld inside the core.

The geomagnetic ﬁeld has apparently been of con-

stant dipole polarity during long ‘‘superchrons’’ (tens of

millions of years). Since this is roughly the time scale for

mantle convection (which is a million times slower than

core convection), it has been suggested that changes in

the thermal structure of the lowermost mantle (due to

the accumulation of subducted lithospheric slabs) may

inﬂuence the geodynamo; see Sec. II.C. The Glatzmaier-

Roberts model has been used to test this hypothesis

(Glatzmaier et al., 1999). Eight cases were simulated, in each of which the total heat ﬂow QCMB from core to
mantle was 7.2 TW. They differed only in the way that

this heat ﬂow was distributed over the CMB, i.e., in the

choice of ICq MB(␪,␾). Spontaneous reversals occurred in

all cases but one. The essence of the results was

(a)

The

westward

drift

is

less

disturbed

when

I

q CMB

is

axisymmetric;

(b) When ICq MB(␪) is unsymmetric with respect to the

equatorial plane, reversals are frequent;

(c)

When

I

q CMB

is

symmetric

with

respect

to

the

equato-

rial

plane

͓ ICq MB(␪)ϭICq MB(␲Ϫ␪) ͔,

and

when

I

q CMB

is

greater in the polar regions than in the equatorial re-

gions (rather than the reverse), the periods of constant

polarity are longer, the duration of reversals (when they

occur) is shorter, and the reversals are between states in

which the ﬁeld has a greater intensity and a weaker

secular variation.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1116

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

FIG. 15. Four snapshots from the Glatzmaier-Roberts simulation spanning a magnetic reversal. The top row shows the radial component of the ﬁeld (red directed outward and blue inward) plotted at the surface of the model Earth. The bottom row shows the longitudinally averaged ﬁeld through the core at the same times. Poloidal ﬁeld is illustrated as lines of force on the left of each image; blue is clockwise directed and yellow is counter-clockwise directed. On the right of each image are contours of toroidal ﬁeld; red is eastward and blue is westward. The outer circle is the core-mantle boundary; the inner circle is the inner core boundary. From Glatzmaier et al., 1999, with permission [Color].

Conditions (c) promote the rise of buoyant ﬂuid in the tangent cylinder and therefore a strong, stable, thermal wind therein. Prescribed ICq MB(␪,␾) that are less compatible with this favored dynamics seem to spawn more successful ﬂuctuations that enhance the secular variation of the ﬁeld and its reversal frequency.
VII. THE FUTURE
This ﬁnal section is devoted to some contentious issues and to some unresolved matters that hopefully will be targets for future research.
A. Turbulence, diffusion, and hyperdiffusion
From a numerical point of view, the necessity for eddy diffusivities of heat, composition, and momentum represents a failure to resolve the solution adequately. If one could gradually improve the resolution, one would also gradually diminish the assumed eddy diffusivities until, when the microscale Lt mainly responsible for the diffusion of the macroscale ﬁelds was reached, it would no longer be necessary to invoke turbulent diffusivities at all. Braginsky and Roberts (1995) estimate that Lt

ϳ2 km and, since the currently attainable resolution is

very much coarser than this, turbulent diffusivities are a

practical necessity. It is hard to estimate the eddy diffusivities ␬t and ␯t.

In addition to the heuristic approach of Braginsky and

Meytlis (1990), two methods have been tried: direct nu-

merical simulation and closure approximation. In the

former, the microscale ﬁelds are solved in isolation in a

small volume of the core, the macroscale ﬁelds being speciﬁed; ␬t is estimated from the statistics of the solu-

tion. See, for example, Matsushima et al. (1999). In the

second approach, the equations governing the mi-

croscale

ﬁelds

␰

t c

,

S

t c

,

Vt,

etc.

are

simpliﬁed

and

solved

using a (possibly drastic) closure approximation, and the

mean

ﬂuxes

͗

␰

t c

Vt

͘

,

͗

S

t c

Vt͘

,

͗ Vt Vt ͘ ,

etc.

are

evaluated

in

terms of the macroscale ﬁelds, as described in Sec. IV.D.

The resulting estimates of ␬t and ␯t are, however, no

better than the closure approximation used to obtain

them. In the absence of reliable estimates of ␬t and ␯t, the

usual expedient is to suppose that all diffusivities,

whether for scalar or for vector ﬁelds, are the same or

are comparable. Most simulations of convective dyna-

mos are of this ‘‘equidiffusional’’ type, but their diffu-

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1117

sivities are many orders of magnitude greater than ␩ ϭ2 m2 sϪ1, so that their Ekman numbers E (viscous, thermal, and magnetic) are typically 10Ϫ3 – 10Ϫ4, rather than 5ϫ10Ϫ9. The Kuang-Bloxham model is equidiffu-
sional, all diffusivities being (for the large scales) about 1.5ϫ104 m2 sϪ1, so that EϷ4ϫ10Ϫ5. All diffusion time scales are therefore at least 7500 times too short relative to the rotation period of the Earth. An alternative way of describing this choice of parameters (J. Bloxham, private communication) is to say that all diffusivities are 2 m2 sϪ1 but that the accelerations experienced by the ﬂuid outer core (and solid inner core) are artiﬁcially damped by 7500 and that the viscous forces are artiﬁcially enhanced by the same factor. In effect, Kuang and Bloxham have increased all their Ekman numbers from 5 ϫ10Ϫ9 to 4ϫ10Ϫ5. However, they have shown (W. Kuang, private communication and Kuang, 1999) that their viscous and inertial terms are at least an order of magnitude smaller than the Coriolis, Lorentz, and buoyancy forces (outside the boundary layers). They therefore claim that the momentum balance has not been seriously disturbed, especially since only the axisymmetric parts of the inertial terms are retained in their model. In the Glatzmaier-Roberts model, ␬t and ␩ are both set to 2 m2 sϪ1, so the thermal and magnetic Ekman numbers are 5ϫ10Ϫ9 (for the large scales). However, to resolve the Ekman boundary layers (Sec. V.B), a much larger value for ␯t, 1500 m2 sϪ1, is used.
The phrase ‘‘for the large scales’’ was used twice in the last paragraph. This was because the choices of diffusivities made there do not overcome the numerical difﬁculties faced by Earth-like simulations. With one exception, these have all required additional ‘‘hyperdiffusion,’’ i.e., diffusivities that increase with decreasing length scale. Hyperdiffusivity allows the large (global) length scales to experience appropriately little diffusion (i.e., to be more Earth-like) and evades the small scales that would otherwise have to be numerically resolved. Solutions with hyperdiffusion appear to be somewhat more Earth-like than those that avoid hyperdiffusion by uniformly increasing the diffusivities on all scales while retaining the same numerical resolution. This may be because the total amount of diffusion is much less. The magnitude and wave-number dependence of hyperdiffusion has, however, so far been quite arbitrary, and the resulting energy spectra may be signiﬁcantly different from those that would be obtained at a much higher spatial resolution using a smaller constant diffusivity. One can avoid hyperdiffusivity by drastically increasing the spatial resolution (Roberts and Glatzmaier, 2000), but such simulations are too expensive to integrate over long times and they still need an eddy viscosity that is almost three orders of magnitude larger than the other diffusivities.
The need to increase ␯ and to add hyperdiffusion has unfortunate repercussions on the energy arguments of Sec. II: Q␯ need no longer be small, and hyperdiffusion is not allowed for in the energy balance (2.3). In a snapshot of the Glatzmaier-Roberts simulation, it was found that QJϷ0.3 TW, QtϷ0.1 TW, and Q␯Ϸ0.7 TW, but QD

should be approximately 2.1 TW according to Sec. II.C. The 1-TW discrepancy seems rather too large to be explained as a ﬂuctuation and is more probably the result of the hyperdiffusion.
Two further difﬁculties should be mentioned. Suppose that, as argued by Braginsky (1999), there is a stable layer (an inverted ‘‘ocean’’) at the top of the core. How should it be simulated numerically? Recall that the energy dissipated by the turbulence is directly supplied by buoyancy forces; see Eq. (4.44). Clearly, by the second law, qt must be positive, but, when the density distribution is bottom heavy, i.e., when the CO density C increases downwards, qt is negative, according to Eq. (4.44). This represents a breakdown in local theory and of the Ansa¨tze (4.48) and (4.50). This breakdown is not surprising: if turbulence is driven only by a locally topheavy density distribution, it cannot arise when the distribution is locally bottom heavy. If there is turbulence in a stable layer, its cause must be different, such as baroclinic instability, or the overshoot of turbulence from an adjacent unstable layer. So far, no serious effort has been made to ﬁnd a substitute for the Ansa¨tze (4.48) and (4.50) inside the stable layer, and we can only hope that the ocean proves to be either nonexistent or unimportant in core MHD.
Second, we took the Reynolds analogy too far in Sec. IV.D when we assumed that, as for molecular motions, the microscale velocities are completely chaotic, and in particular have no preferred direction. Because of the effects of the Coriolis and Lorentz forces, core turbulence is expected to be highly anisotropic, the convective eddies being elongated in the directions of ⍀ and B; see Braginsky (1964d), Braginsky and Meytlis (1990), Braginsky and Roberts (1995), and St. Pierre (1996). This means that turbulent transport should also be enhanced in these directions, so that the isotropic laws (4.48) and (4.50) should be replaced by

I␰ϭϪ␳␬I•“␰c ,

(7.1)

ISϭϪ␳␬I•“Sc ,

(7.2)

where the diffusivity, now the tensor ␬I, is as before the same for ␰c as for Sc .
Anisotropy considerably complicates momentum transfer. In principle a fourth-rank tensor I␯ substitutes for the scalar ␯, and Eq. (4.55) is replaced by

␲i␯jϭ␳␯ijklekl .

(7.3)

Since eijϭeji , the symmetries ␯jiklϭ␯ijlkϭ␯ijkl can be assumed, and the requirement q␯ϵeij␯ijkleklу0 of positive entropy production demands a further reduction to

‘‘only’’ 21 independent components that have to be es-

timated. There is an added complication: turbulent

transport of vector ﬁelds, such as magnetic ﬁeld or mo-

mentum, introduces more than simple diffusion. For ex-

ample

(Sec.

III.C),

in

addition

to

a

contribution

␩

t
“

ÃB

to the turbulent emf F, the ␣ effect introduces a term

␣B, proportional to the macroscale B and not to its gra-

dient. In an analogous way, studies of turbulence in the

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1118

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

Sun have led to the recognition of the ‘‘⌳ effect’’ (aka the ‘‘AKA effect,’’ standing for anisotropic kinetic alpha effect), in which

␲i␯jϭ␳⌳ijk␻k ,

(7.4)

where ␻ϭ“ÃV is the local macroscale vorticity; see Chap. 4 of Ru¨ diger (1989). When this operates, as it does in solar convection, the macroscale ﬂow cannot be solid-body rotation.
The importance of these tensor diffusivities in computer simulations of the geodynamo again depends on the spatial resolution one can afford. Presumably the more modes that are numerically resolved, the less important are the details of the transport by the remaining modes until, when the resolved part of the spectrum contains most of the energy, all anisotropic turbulent transport will be adequately represented. Until that day arrives, one can hope that, since RoӶ1, turbulent momentum transport is not signiﬁcant in the core, except in boundary layers and possibly near the tangent cylinder. If so, it may not be necessary to apply the Ansa¨tze (7.3) and (7.4). It is doubtful whether the Ansa¨tze (7.1) and (7.2) can be similarly ignored, but so far no numerical simulations have taken them into account.

B. Boundary conditions

Away from boundaries, the eddy viscosity plays a mi-

nor role in both the Glatzmaier-Roberts and the Kuang-

Bloxham models, but it is still orders of magnitude

greater than what one wants. The Ekman layers on the

core-mantle and inner core boundaries present the dy-

namo theorist with something of a dilemma. On the one

hand, if he takes ␯ to be the molecular viscosity ␯M , he ﬁnds that ␦␯ is only about 10 cm and that the pumping velocities E1/2U¯ are only about 10Ϫ10 m sϪ1; see Sec. V.B.

He may ﬁnd it hard to believe that these can be signiﬁ-

cant, and he may agree with Kuang and Bloxham (1999)

that they should be ignored, and the no-slip conditions

(4.59) and (4.60) replaced by the conditions that the tan-

gential

components

␲

␯ r␪

and

␲

␯ r␾

of

viscous

stress

vanish—the so-called viscous stress-free conditions.

On the other hand, such a step would remove the Ek-

man layers that are responsible for spin-up (Sec. V.B)

and for helping to create an ␣ effect. Moreover, the use

of ␯M in computing the Ekman pumping is questionable. The Ekman layers observed in the upper layers of the

ocean are described by a turbulent viscosity about three

orders of magnitude greater than the molecular viscosity

of seawater (Hunkins, 1966). This suggests that we

should similarly describe the Ekman layers in the core with the turbulent viscosity ␯tϷ2 m2 sϪ1, leading to ␦␯ Ϸ160 m and a spin-up time ␶su of only about 9 yr (Sec.
V.B); the viscous coupling between the ﬂuid and its

boundaries is especially potent for the geostrophic part

of V. In contrast, the viscous stress-free conditions give

␶suϭϱ, and only magnetic, topographic, and gravitational torques can couple the ﬂuid to its boundaries.

This suggests that use of the viscous stress-free condi-

tions in place of the no-slip conditions is a more serious step than one might intuitively have suspected.
No-slip conditions were used in the GlatzmaierRoberts simulations mainly for consistency: since the transport of momentum by the unresolved eddies is included within the ﬂuid, why ignore the stronger coupling between the ﬂuid and the solid boundaries? Nevertheless, the eddy viscous torque on the boundaries of the Glatzmaier-Roberts simulations is too large, so the other extreme of eliminating it totally, as in the KuangBloxham model, is certainly worth investigating. Test runs with the Glatzmaier-Roberts model with zero viscous torque on the boundaries produced, as in the Kuang-Bloxham solution, a smaller inner core superrotation rate and less magnetic-ﬁeld generation near the inner core boundary.
A comparison of the magnetic ﬁelds generated in the Kuang-Bloxham simulation and those in the current Glatzmaier-Roberts solution (with no-slip boundaries) is made in Fig. 16. Snapshots of the radial component of the generated ﬁelds, plotted at the CMB, and the axisymmetric parts of the ﬁelds, plotted throughout the core, are displayed. Both solutions produce dipoledominated magnetic ﬁelds at the CMB. The lower diffusion and higher resolution of the current GlatzmaierRoberts model, however, produces a solution with a more prominent small-scale magnetic structure and a more distinct tangent cylinder effect. The zonal ﬁeld is mainly conﬁned to the interior of the tangent cylinder, whereas the meridional ﬁeld is mainly outside it. The Kuang-Bloxham solution and the original (lowresolution) Glatzmaier-Roberts solution are dominated by larger-scale structures.
C. The road ahead
It will be apparent to all readers who have had the stamina to read this review through that geodynamo theory already explains, at least qualitatively, virtually every known facet of the geomagnetic ﬁeld described in Sec. I.C. They will also recognize that this is remarkable in view of the approximations and uncertainties currently inherent in the theory, some more of which are adumbrated below.
First, most of the physical parameters on which simulations rely are uncertain to a greater or lesser extent, especially the strength of the radioactive heat sources in the core, the heat ﬂow from the core to the mantle, the latent heat of crystallization of core ﬂuid, and the depression of its freezing point by alloying. This is hardly surprising when we recall that the principal alloying element in the core is still unknown and is likely to remain so for some time to come. Interesting issues are raised by the apparent youth of the solid inner core. According to the simulation of Glatzmaier and Roberts, ¯␶Ϸ1.3 Gyr, and this is not very different from other recent

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1119

FIG. 16. Snapshots of the radial component of the magnetic ﬁeld at the core-mantle boundary (in equal-area projections with red outward and blue inward) and the longitudinally averaged ﬁeld throughout the core (with toroidal ﬁeld contours on the left and poloidal ﬁeld lines of force on the right of each image). The top row is from the Kuang-Bloxham simulation and the bottom row is from the Glatzmaier-Roberts simulation [Color].

estimates.36 The age of the geomagnetic ﬁeld exceeds 3 Gyr, however (Sec. I.C). A dynamo can certainly operate in a completely ﬂuid core (Sakuraba and Kono, 1999), but one might have expected it to produce a qualitatively different, and perhaps weaker, magnetic ﬁeld than the one observed today. Even if we took QD ϭ0, we would, if we insisted that QRϭ0, ﬁnd from Sec. III.C that ¯␶ is only 3.7 Gyr. Moreover, QDϭ0 is unreasonable since the convection that is needed to maintain the adiabatic state of the core must dissipate energy
36Stevenson et al. (1983) found that ¯␶Ϸ2.3 Gyr, Buffett et al. (1996) that ¯␶Ϸ2.8 Gyr, and Labrosse et al. (1997) that ¯␶ Ϸ1.7 Gyr. These authors modeled the thermal history of the Earth globally and did not set up a theory on which a detailed model could be constructed. Consequently, they omitted the turbulent dissipation, which, following Braginsky and Roberts (1995), is included by Glatzmaier and Roberts.

(and, if there is no magnetic ﬁeld, this involves a very large viscous dissipation; see Sec. V.C). Is it really true that the SIC is much younger than the Earth? Or is QR signiﬁcantly nonzero? Or are some of our other parameter estimates wide of the mark?
Second, spectral transform methods are currently used by virtually every modeler. Despite their many virtues, it is hard to attain the numerical resolution necessary to represent solutions in such a lightly diffused system as the core; large eddy diffusivities have to be used. Presumably, as computer technology advances and becomes cheaper, resolution will improve and smaller diffusivities will become practical. Eventually, at some time in the remote future, it will be possible to construct models using only molecular diffusivities. In the meantime, how is progress to be maintained? One possibility is to continue to parametrize the effects of the turbulent microscales on the global macroscales and to hope that

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1120

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

advances in turbulence theory will lead to improved parametrizations. A second method, currently favored by several new groups, is to use very different numerical methods that are more suitable for lightly damped systems and are more readily applied in parallel computation than spectral transform methods. Some exploratory work has already been done in this direction, but only for 2D mean-ﬁeld models (Braginsky and Roberts, 1994; Jault, 1995).
Improved numerical techniques should help answer other puzzles, such as the roles of Taylor’s constraint and model z in geodynamo theory (Sec. V.D). Current Earth-like models cannot provide an unequivocal answer because their Ekman numbers E are too large, but the question has been given increasing attention over the past decade through the integration of 2D meanﬁeld models. It appears that some models are Taylorlike, some have a model-z character, while others are of neither type; see Hollerbach (1996a).
The effect of adding new physical ingredients should be investigated. For example (Sec. II.A), the effects of the luni-solar precession of the Earth’s mantle needs to be studied with 3D geodynamo simulations.
Finally, models similar to the geodynamo models reviewed here form a basis for the study of dynamos in other planets and satellites in our solar system and in planets now being discovered around other stars, and we may expect new research in this direction also.
Geodynamo simulations have made giant strides during the past decade and, with more players joining this exciting activity every year and with ever increasing computer capabilities, we are conﬁdent that many of the outstanding issues will be resolved during the ﬁrst decade of the new millennium.
ACKNOWLEDGMENTS
We thank Robert Coe and David Loper for their criticisms of a draft of this review, and Weijia Kuang for supplying half of Fig. 16, and for his comments on parts of Sec. VII. We thank Maureen Roberts for producing Figs. 1, 6, 7, and 8. Our work on the geodynamo is supported by the Institute of Geophysics and Planetary Physics, the Los Alamos LDRD Program, by the University of California Research Partnership Initiatives Program, by the NSF Geophysics Program, and by the NASA HPCC/ESS Grand Challenge Program. Computing resources were provided by the NPACI, NCSA, GSFC, and the Los Alamos ACL supercomputing centers.
REFERENCES
Acheson, D. J., 1972, ‘‘On the hydromagnetic stability of a rotating ﬂuid annulus,’’ J. Fluid Mech. 52, 529–541.
Acheson, D. J., and R. Hide, 1973, ‘‘Hydromagnetics of rotating ﬂuids,’’ Rep. Prog. Phys. 36, 159–221.
Bender, C. M., and S. A. Orszag, 1978, Advanced Mathematical Methods for Scientists and Engineers (McGraw-Hill, New York).

Birch, F., 1952, ‘‘Elasticity and constitution of the Earth’s interior,’’ J. Geophys. Res. 57, 227–286.
Blackett, P. M. S., 1947, ‘‘The magnetic ﬁeld of massive rotating bodies,’’ Nature (London) 159, 658–666.
Braginsky, S. I., 1963, ‘‘Structure of the F layer and reasons for convection in earth’s core,’’ Dokl. Akad. Nauk SSSR 149, 1311–1314 [Sov. Phys. Dokl. 149, 8–10 (1963)].
Braginsky, S. I., 1964a, ‘‘Self-excitation of a magnetic ﬁeld during the motion of a highly conducting ﬂuid,’’ Zh. Eksp. Teor. Fiz. SSSR 47, 1084–1098 [Sov. Phys. JETP 20, 726–735].
Braginsky, S. I., 1964b, ‘‘Theory of the hydromagnetic dynamo,’’ Zh Eksp. Teor. Fiz. SSSR 47, 2178–2193 [Sov. Phys. JETP 20, 1462–1471 (1965)].
Braginsky, S. I., 1964c, ‘‘Kinematic models of the Earth’s hydromagnetic dynamo,’’ Geomagn. Aeron. (Russia) 4, 732– 747 [Geomagn. Aeron. (USA) 4, 572–583 (1964)].
Braginsky, S. I., 1964d, ‘‘Magnetohydrodynamics of the Earth’s core,’’ Geomagn. Aeron. (Russia) 4, 898–916 [Geomagn. Aeron. (USA) 4, 698–712 (1964)].
Braginsky, S. I., 1967, ‘‘Magnetic waves in the Earth’s core,’’ Geomagn. Aeron. (Russia) 7, 1050–1060 [Geomagn. Aeron. (USA) 7, 851–859 (1967)].
Braginsky, S. I., 1970, ‘‘Torsional magnetohydrodynamic vibrations in the Earth’s core and variations in day length,’’ Geomagn. Aeron. (Russia) 10, 3–12 [Geomagn. Aeron. (USA) 10, 1–8 (1970)].
Braginsky, S. I., 1975, ‘‘Nearly axially symmetric model of the hydromagnetic dynamo of the Earth, I,’’ Geomagn. Aeron. (Russia) 15, 149–156 [Geomagn. Aeron. (USA) 15, 122–128 (1975)].
Braginsky, S. I., 1978, ‘‘Nearly axially symmetric model of the hydromagnetic dynamo of the Earth,’’ Geomagn. Aeron. (Russia) 18, 340–351 [Geomagn. Aeron. (USA) 18, 225–231 (1978)].
Braginsky, S. I., 1999, ‘‘Dynamics of the stably stratiﬁed ocean at the top of the core,’’ Phys. Earth Planet. Inter. 111, 21–34.
Braginsky, S. I., and V. P. Meytlis, 1990, ‘‘Local turbulence in the Earth’s core,’’ Geophys. Astrophys. Fluid Dyn. 55, 71–87.
Braginsky, S. I., and P. H. Roberts, 1994, ‘‘From Taylor state to model-z,’’ Geophys. Astrophys. Fluid Dyn. 77, 3–13.
Braginsky, S. I., and P. H. Roberts, 1995, ‘‘Equations governing Earth’s core and the geodynamo,’’ Geophys. Astrophys. Fluid Dyn. 79, 1–97.
Braginsky, S. I., and P. H. Roberts, 2000, in Advances in Nonlinear Dynamos, edited by A. Ferriz Mas and M. Nu˜nez (Gordon and Breach, London).
Brandenburg, A., A. Nordland, P. Pulkkinen, R. F. Stein, and I. Tuominen, 1990, ‘‘Simulation of turbulent cyclonic magneto-convection,’’ Astron. Astrophys. 232, 277–291.
Buffett, B. A., 1996, ‘‘A mechanism for decade ﬂuctuations in the length of day,’’ Geophys. Res. Lett. 23, 3803–3806.
Buffett, B. A., 1997, ‘‘Geodynamic estimates of the viscosity of the Earth’s inner core,’’ Nature (London) 388, 571–573.
Buffett, B. A., and G. A. Glatzmaier, 2000, ‘‘Gravitational breaking of inner-core rotation in geodynamo simulations,’’ Geophys. Res. Lett., in press.
Buffett, B. A., H. E. Huppert, J. R. Lister, and A. W. Woods, 1996, ‘‘On the thermal evolution of the Earth’s core,’’ J. Geophys. Res. 101, 7989–8006.
Bullard, E. C., 1949, ‘‘The magnetic ﬁeld within the earth,’’ Proc. R. Soc. London, Ser. A 197, 433–453.
Busse, F. H., 1970, ‘‘Thermal instability in rapidly rotating systems,’’ J. Fluid Mech. 44, 441–460.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1121

Busse, F. H., 1975, ‘‘A model of the geodynamo,’’ Geophys. J. R. Astron. Soc. 42, 437–459.
Busse, F. H., 2000, ‘‘Homogeneous dynamos in planetary cores and in the laboratory,’’ Annu. Rev. Fluid Mech. 32, 383–408.
Busse, F. H., E. Grote, and A. Tilgner, 1998, ‘‘On convection driven dynamos in rotating spherical shells,’’ Stud. Geophys. Geod. 42, 1–6.
Chabot, N. L., and M. J. Drake, 1999, ‘‘Potassium solubility in metal: The effects of composition at 15 kbar and 1900° on partitioning between iron alloys and silicate melts,’’ Earth Planet. Sci. Lett. 172, 323–335.
Childress, S., and A. D. Gilbert, 1995, Stretch, Twist, Fold: The Fast Dynamo (Springer, Berlin).
Childress, S., and A. M. Soward, 1972, ‘‘Convection-driven hydromagnetic dynamo,’’ Phys. Rev. Lett. 29, 837–839.
Christensen, U., P. Olson, and G. A. Glatzmaier, 1998, ‘‘A dynamo model interpretation of geomagnetic ﬁeld structures,’’ Geophys. Res. Lett. 25, 1565–1568.
Christensen, U., P. Olson, and G. A. Glatzmaier, 1999, ‘‘Numerical modeling of the geodynamo: A systematic parameter study,’’ Geophys. J. Int. 138, 393–409.
Cowling, T. G., 1933, ‘‘The magnetic ﬁeld of sunspots,’’ Mon. Not. R. Astron. Soc. 140, 39–48.
Creager, K., 1997, ‘‘Inner core rotation rate from small scale heterogeneity and time-varying travel times,’’ Science 278, 1284–1288.
de Wijs, G. A., G. Kresse, L. Voc˘adlo, D. Dobson, A. Alfe` , M. J. Gillan, and G. D. Price, 1998, ‘‘The viscosity of liquid iron at the physical conditions of the Earth’s core,’’ Nature (London) 392, 805–807.
Dudley, M. L., and R. W. James, 1989, ‘‘Time-dependent kinematic dynamos with stationary ﬂows,’’ Proc. R. Soc. London, Ser. A 425, 407–429.
Dziewonski, A. M., and D. L. Anderson, 1981, ‘‘Preliminary reference Earth model,’’ Phys. Earth Planet. Inter. 25, 297– 356.
Elsasser, W. M., 1946, ‘‘Induction effects in terrestrial magnetism. Part II. The secular variation,’’ Phys. Rev. 70, 202–212.
Fautrelle, Y., and S. Childress, 1982, ‘‘Convective dynamos with intermediate and strong ﬁelds,’’ Geophys. Astrophys. Fluid Dyn. 22, 235–279.
Fearn, D. R., 1998, ‘‘Hydromagnetic ﬂow in planetary cores,’’ Rep. Prog. Phys. 61, 175–235.
Fearn, D. R., and M. R. E. Proctor, 1987, ‘‘On the computation of steady, self-consistent spherical dynamos,’’ Geophys. Astrophys. Fluid Dyn. 38, 293–325.
Fowler, T. K., 1999, ‘‘Nuclear power—fusion,’’ Rev. Mod. Phys. 71, S456–459.
Gailitis, A., 1993, ‘‘Magnetic ﬁeld generation by the axisymmetric conducting ﬂuid ﬂow in a spherical cavity of a stationary conductor,’’ Magn. Gidrodin. 29, 3–14 [Magnetohydrodynamics 29, 107–115 (1993)].
Gailitis, A., O. Lielausis, S. Dementev, E. Platacis, A. Ciferons, G. Gerbeth, T. Gundrum, F. Stefani, M. Christen, H. Ha¨ nel, and G. Will, 2000, ‘‘Detection of a ﬂow induced magnetic ﬁeld eigenmode in the Riga dynamo facility,’’ Phys. Rev. Lett. 84, 4365–4368.
Gill, A. E., 1982, Atmosphere-Ocean Dynamics (Academic, New York).
Gilman, P. A., 1983, ‘‘Dynamically consistent nonlinear dynamos driven by convection in a rotating spherical shell. II. Dynamos with cycles and strong feedbacks,’’ Astrophys. J., Suppl. Ser. 53, 243–268.

Gilman, P. A., and J. Miller, 1981, ‘‘Dynamically consistent nonlinear dynamos driven by convection in a rotating spherical shell,’’ Astrophys. J., Suppl. Ser. 46, 211–238.
Gire, C., J.-L. Le Moue¨ l, and J. Ducruix, 1983, ‘‘Evolution of the secular variation ﬁeld from the beginning of the century,’’ Nature (London) 307, 349–352.
Glatzmaier, G. A., 1984, ‘‘Numerical simulations of stellar convective dynamos. I. The model and the method,’’ J. Comput. Phys. 55, 461–484.
Glatzmaier, G. A., 1985a, ‘‘Numerical simulations of stellar convective dynamos. II. Field propagation in the convection zone,’’ Astrophys. J. 291, 300–307.
Glatzmaier, G. A., 1985b, ‘‘Numerical simulations of stellar convective dynamos. III. At the bottom of the convection zone,’’ Geophys. Astrophys. Fluid Dyn. 31, 137–150.
Glatzmaier, G. A., R. S. Coe, L. Hongre, and P. H. Roberts, 1999, ‘‘The role of the Earth’s mantle in controlling the frequency of geomagnetic reversals,’’ Nature (London) 401, 885–890.
Glatzmaier, G. A., and P. H. Roberts, 1995a, ‘‘A threedimensional convective dynamo solution with rotating and ﬁnitely conducting inner core and mantle,’’ Phys. Earth Planet. Inter. 91, 63–75.
Glatzmaier, G. A., and P. H. Roberts, 1995b, ‘‘A threedimensional self-consistent computer simulation of a geomagnetic ﬁeld reversal,’’ Nature (London) 377, 203–209.
Glatzmaier, G. A., and P. H. Roberts, 1996a, ‘‘An anelastic evolutionary geodynamo simulation driven by compositional and thermal convection,’’ Physica D 97, 81–94.
Glatzmaier, G. A., and P. H. Roberts, 1996b, ‘‘Rotation and magnetism of Earth’s inner core,’’ Science 274, 1887–1891.
Glatzmaier, G. A., and P. H. Roberts, 1997, ‘‘Simulating the geodynamo,’’ Contemp. Phys. 38, 269–288.
Glatzmaier, G. A., and P. H. Roberts, 1998, ‘‘Dynamo theory then and now,’’ Int. J. Eng. Sci. 36, 1325–1338.
Greenspan, H. P., 1968, The Theory of Rotating Fluids (Cambridge University, Cambridge, England).
Grote, E., F. H. Busse, and A. Tilgner, 1999, ‘‘Convectiondriven quadrupolar dynamos in rotating spherical shells,’’ Phys. Rev. E 60, R5025–R5028.
Grote, E., F. H. Busse, and A. Tilgner, 2000, ‘‘Regular and chaotic spherical dynamos,’’ Phys. Earth Planet. Inter. 117, 259–272.
Gubbins, D., and P. H. Roberts, 1987, in Geomagnetism, 2, edited by J. A. Jacobs (Academic, London, England), pp. 1–183.
Hide, R., 1966, ‘‘Free oscillations of the Earth’s core and the theory of the geomagnetic secular variation,’’ Philos. Trans. R. Soc. London, Ser. A 259, 615–647.
Hide, R., 1978, ‘‘How to locate the electrically conducting ﬂuid core of a planet from external magnetic observation,’’ Nature (London) 271, 640–641.
Hide, R., and S. R. C. Malin, 1979, ‘‘The size of Jupiter’s electrically conducting ﬂuid core,’’ Nature (London) 280, 42–43.
Hide, R., and S. R. C. Malin, 1981, ‘‘On the determination of the size of the Earth’s core from observations of the geomagnetic secular variation,’’ Proc. R. Soc. London, Ser. A 374, 15–33.
Hollerbach, R., 1994, ‘‘Magnetohydrodynamic Ekman and Stewartson layers in a rotating spherical shell,’’ Proc. R. Soc. London, Ser. A 444, 333–346.
Hollerbach, R., 1996a, ‘‘On the theory of the geodynamo,’’ Phys. Earth Planet. Inter. 98, 163–185.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

1122

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

Hollerbach, R., 1996b, ‘‘Magnetohydrodynamic shear layers in a rapidly rotating plane layer,’’ Geophys. Astrophys. Fluid Dyn. 82, 237–253.
Hollerbach, R., and C. A. Jones, 1993, ‘‘Inﬂuence of the Earth’s core on geomagnetic ﬂuctuations and reversals,’’ Nature (London) 365, 541–543.
Hughes, D. W., 1993, in Solar and Planetary Dynamos, edited by M. R. E. Proctor, P. C. Matthews, and A. M. Rucklidge (Cambridge University, Cambridge, England), pp. 153–159.
Hunkins, K., 1966, ‘‘Ekman drift currents in the Arctic Ocean,’’ Deep-Sea Res. 13, 607–620.
Ivers, D. J., and R. W. James, 1984, ‘‘Axisymmetric antidynamo theorems in compressible non-uniform conducting ﬂuids,’’ Philos. Trans. R. Soc. London, Ser. A 312, 179–218.
Jacobs, J. A., 1953, ‘‘The Earth’s inner core,’’ Nature (London) 172, 297–300.
Jackson, A., J. Bloxham, and D. Gubbins, 1993, in Dynamics of Earth’s Deep Interior and Earth Rotation, edited by J.-L. Le Moue¨ l, D. E. Smylie, and T. Herring, Geophys. Mongr. Sec. No. 72, IUGG Vol. 12 (Am. Geophys. Union, Washington, D.C.), p. 97.
Jault, D., 1995, ‘‘Model z by computation and Taylor’s condition,’’ Geophys. Astrophys. Fluid Dyn. 79, 99–124.
Jault, D., G. Gire, and J.-L. Le Moue¨ l, 1988, ‘‘Westward drift, core motions and exchanges of angular momentum between core and mantle,’’ Nature (London) 333, 353–356.
Jones, C. A., A. W. Longbottom, and R. Hollerbach, 1995, ‘‘A self-consistent convection-driven geodynamo model, using a mean-ﬁeld approximation,’’ Phys. Earth Planet. Inter. 92, 119–141.
Jones, C. A., and P. H. Roberts, 2000, ‘‘Convection-driven dynamos in a rotating plane layer,’’ J. Fluid Mech. 404, 311–343.
Jones, C. A., A. M. Soward, and A. I. Mussa, 2000, ‘‘The onset of thermal convection in a rapidly rotating sphere,’’ J. Fluid Mech. 405, 157–179.
Kabantsev, A. A., V. B. Reva, and V. G. Sokolov, 1999, ‘‘Experimental veriﬁcation of the MHD dynamo in the axisymmetric linear machine,’’ Fusion Technol. 35, 185–189.
Kageyama, A., M. M. Ochi, and T. Sato, 1999, ‘‘Flip-ﬂop transitions of the magnetic intensity and polarity reversals in the magnetohydrodynamic dynamo,’’ Phys. Rev. Lett. 82, 5409– 5412.
Kageyama, A., and T. Sato, 1995, ‘‘Computer simulation of a magnetohydrodynamic dynamo II,’’ Phys. Plasmas 2, 1421– 1431.
Kageyama, A., and T. Sato, 1997, ‘‘Generation mechanism of a dipole ﬁeld by a magneto-hydrodynamic dynamo,’’ Phys. Rev. E 55, 4617–4626.
Kida, S., K. Araki, and H. Kitauchi, 1997, ‘‘Periodic reversals of magnetic ﬁeld generated by thermal convection in a rotating spherical shell,’’ J. Phys. Soc. Jpn. 66, 2194–2201.
Kitauchi, H., and S. Kida, 1998, ‘‘Intensiﬁcation of magnetic ﬁeld by concentrate-and-stretch of magnetic ﬂux lines,’’ Phys. Fluids 10, 457–468.
Kono, M., and H. Tanaka, 1995, in The Earth’s Central Part: Its Structure and Dynamics, edited by T. Yukutake (Terrapub, Tokyo, Japan), pp. 75–94.
Krause, F., and K.-H. Ra¨ dler, 1980, Mean-ﬁeld Magnetohydrodynamics and Dynamo Theory (Pergamon, Oxford, England).
Kuang, W., 1999, ‘‘Force balances and convective state in the Earth’s core,’’ Phys. Earth Planet. Inter. 116, 65–79.

Kuang, W., and J. Bloxham, 1997, ‘‘An Earth-like numerical dynamo model,’’ Nature (London) 389, 371–374.
Kuang, W., and J. Bloxham, 1999, ‘‘Numerical modeling of magnetohydrodynamic convection in a rapidly rotating spherical shell: Weak and strong ﬁeld dynamo action.’’ J. Comput. Phys. 153, 51–81.
Kutzner, C., and U. Christensen, 2000, ‘‘Effects of driving mechanisms in geodynamo models,’’ Geophys. Res. Lett. 27, 29–32.
Labrosse, S., J.-P. Poirier, and J.-L. Moue¨ l, 1997, ‘‘On the cooling of Earth’s core,’’ Phys. Earth Planet. Inter. 99, 1–17.
Landau, L. D., and E. M. Lifshitz, 1987, Fluid Mechanics, 2nd Edition (Pergamon, Oxford, England).
Langel, R. A., 1987, in Geomagnetism, Vol. I, edited by J. A. Jacobs (Academic, London), pp. 249–512.
Langel, R. A., and R. H. Estes, 1982, ‘‘A geomagnetic ﬁeld spectrum,’’ Geophys. Res. Lett. 9, 250–253.
Langel, R. A., and R. H. Estes, 1985, ‘‘The near-Earth magnetic ﬁeld at 1980 determined from MAGSAT data,’’ J. Geophys. Res. 90, 2495–2509.
Larmor, J., 1919, ‘‘How could a rotating body such as the Sun become a magnet?’’ Rep. Brit. Assoc. Adv. Sci. 159–160.
Le Moue¨ l, J.-L., and V. Courtillot, 1981, ‘‘Core motions, electromagnetic core-mantle coupling and variations in earth’s rotation: new constraints from geomagnetic secular variation impulses,’’ Phys. Earth Planet. Inter. 24, 236–241.
Le Moue¨ l, J.-L., J. Ducruix, and C. H. Duyen, 1982, ‘‘The worldwide character of the 1969–1970 impulse of the secular variation rate,’’ Phys. Earth Planet. Inter. 28, 337–350.
Le Moue¨ l, J.-L., T. R. Madden, J. Ducruix, and V. Courtillot, 1981, ‘‘Decade ﬂuctuations in geomagnetic westward drift and earth rotation,’’ Nature (London) 290, 763–765.
Loper, D. E., 1970, ‘‘On the unsteady hydromagnetic perturbations at the Earth’s core-mantle interface,’’ Phys. Earth Planet. Inter. 4, 129–137.
Loper, D. E., 1978, ‘‘The gravitationally powered dynamo,’’ Geophys. J. R. Astron. Soc. 54, 389–404.
Loper, D. E., and P. H. Roberts, 1981, ‘‘A study of conditions at the inner core boundary of the Earth,’’ Phys. Earth Planet. Inter. 24, 302–307.
Love, J. J., and D. Gubbins, 1996, ‘‘Optimized kinematic dynamos,’’ Geophys. J. Int. 124, 787–800.
Lund, S. P., G. Acton, B. Clement, M. Hastedt, M. Okada, and T. Williams, 1998, ‘‘Geomagnetic ﬁeld excursions occurred often during the last million years,’’ EOS Trans. Am. Geophys. Union 79, 178–179.
Malkus, W. V. R., 1963, ‘‘Precessional torques as the cause of geomagnetism,’’ J. Geophys. Res. 68, 2871–2886.
Malkus, W. V. R., 1968, ‘‘Precession of the Earth as the cause of geomagnetism,’’ Science 160, 259–264.
Matsushima, M., T. Nakajima, and P. H. Roberts, 1999, ‘‘The anisotropy of local turbulence in the Earth’s core,’’ Earth Planets Space 51, 277–286.
McElhinny, M. W., 1973, Paleomagnetism and Plate Tectonics (Cambridge University, Cambridge, England).
Meneguzzi, M., and A. Pouquet, 1989, ‘‘Turbulent dynamos driven by convection,’’ J. Fluid Mech. 205, 297–318.
Merrill, R. T., M. W. McElhinny, and P. L. McFadden, 1996, The Magnetic Field of the Earth, Paleomagnetism, the Core and the Deep Mantle (Academic, London).
Moffatt, H. K., 1969, ‘‘The degree of knottedness of tangled vortex lines,’’ J. Fluid Mech. 35, 117–129.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

P. H. Roberts and G. A. Glatzmaier: Geodynamo theory and simulations

1123

Moffatt, H. K., 1978, Magnetic Field Generation in Electrically

Conducting Fluids (Cambridge University, Cambridge, En-

gland).

Monin, A. S., 1990, Theoretical Geophysical Fluid Dynamics

(Kluwer, Dordrecht).

Morrison, G., and D. R. Fearn, 2000, ‘‘The inﬂuence of Ray-

leigh number, azimuthal wave number and inner core radius

on

2

1 2

D

hydromagnetic

dynamos,’’

Phys.

Earth

Planet.

Inter.

117, 237–258.

Olson, P., U. Christensen, and G. A. Glatzmaier, 1999, ‘‘Nu-

merical modeling of the geodynamo: Mechanisms of ﬁeld

generation and equilibration,’’ J. Geophys. Res. 104, 10383–

10404.

Parker, E. N., 1955, ‘‘Hydromagnetic dynamo models,’’ Astro-

phys. J. 122, 293–314.

Pedlosky, J., 1979, Geophysical Fluid Dynamics (Springer,

New York).

Ponomarenko, Y. B., 1973, ‘‘On the theory of hydromagnetic

dynamos,’’ Zh. Prikl. Mekh. Tekh. Fiz. #6, 47–51 [J. Appl.

Mech. Tech. Phys. 14, 775–778 (1973)].

Ricard, Y., M. Richards, C. Lithgow-Bertolloni, and Y. LeS-

tunff, 1993, ‘‘A geodynamic model of mantle density hetero-

geneity,’’ J. Geophys. Res. 98, 21859–21909.

Roberts, P. H., 1968, ‘‘On the thermal instability of a rotating-

ﬂuid sphere containing heat sources,’’ Philos. Trans. R. Soc.

London, Ser. A 263, 93–117.

Roberts, P. H., 1972, ‘‘Kinematic dynamo models,’’ Philos.

Trans. R. Soc. London, Ser. A 272, 663–703.

Roberts, P. H., and G. A. Glatzmaier, 2000, ‘‘A test of the

frozen ﬂux approximation using geodynamo simulations,’’

Philos. Trans. R. Soc. London, Ser. A 358, 1109–1121.

Roberts, P. H., and T. H. Jensen, 1993, ‘‘Homogeneous dyna-

mos: Theory and practice,’’ Phys. Fluids B 5, 2657–2662.

Roberts, P. H., and A. M. Soward, 1972, ‘‘Magnetohydrody-

namics of the Earth’s core,’’ Annu. Rev. Fluid Mech. 4, 117–

152.

Roberts, P. H., and A. M. Soward, 1978, eds., Rotating Fluids

in Geophysics (Academic, New York).

Roberts, P. H., and K. Stewartson, 1975, ‘‘On double-roll con-

vection in a rotating magnetic system,’’ J. Fluid Mech. 68,

447–466.

Roberts, P. H., and M. Stix, 1971, The Turbulent Dynamo. A Translation of a Series of Papers by F. Krause, K.-H. Ra¨dler

and M. Steenbeck (Report TN/IA-60, National Center for At-

mospheric Research, Boulder, Colorado). Ru¨ diger, G., 1989, Differential Rotation and Stellar Convection:

Sun and Solar-type Stars (Gordon and Breach, New York).

Sakuraba, A., and M. Kono, 1999, ‘‘Effect of the inner core on

the numerical solution of the magnetohydrodynamic dy-

namo,’’ Phys. Earth Planet. Inter. 111, 105–121.

Sarson, G. R., 2000, ‘‘Reversal models from dynamo calcula-

tions,’’ Philos. Trans. R. Soc. London, Ser. A 358, 921–942.

Sarson, G. R., and C. A. Jones, 1999, ‘‘A convection driven

geodynamo reversal model,’’ Phys. Earth Planet. Inter. 111,

3–20.

Sarson, G. R., C. A. Jones, and A. W. Longbottom, 1997, ‘‘The

inﬂuence of boundary region heterogeneities on the geody-

namo,’’ Phys. Earth Planet. Inter. 101, 13–32.

Sarson, G. R., C. A. Jones, and A. W. Longbottom, 1998,

‘‘Convection-driven geodynamo models of varying Ekman

number,’’ Geophys. Astrophys. Fluid Dyn. 88, 225–259.

Song, X., and P. Richards, 1996, ‘‘Seismological evidence for differential rotation of the Earth’s inner core,’’ Nature (London) 382, 221–224.
Souriau, A., 1998, ‘‘New seismological constraints on differential rotation of the inner core from Novaya Zemlya events recorded at DRV, Antarctica,’’ Geophys. J. Int. 134, F1–F5.
Soward, A. M., 1972, ‘‘A kinematic theory of large magnetic Reynolds number dynamos,’’ Philos. Trans. R. Soc. London, Ser. A 272, 431–462.
Soward, A. M., 1974, ‘‘A convection-driven dynamo I. Weakﬁeld case,’’ Philos. Trans. R. Soc. London, Ser. A 275, 611– 651.
Soward, A. M., 1990, ‘‘A uniﬁed approach to a class of slow dynamos,’’ Geophys. Astrophys. Fluid Dyn. 53, 81–107.
Stacey, F. D., 1992, Physics of the Earth, 3rd ed. (Brookﬁeld, Brisbane, Australia).
Steenbeck, M., and F. Krause, 1966, ‘‘The generation of stellar and planetary magnetic ﬁelds by turbulent dynamo action,’’ Z. Naturforsch. A 21A, 1285–1296. (English translation in Roberts and Stix, 1971).
Stevenson, D. J., 1984, ‘‘The energy ﬂux number and three types of planetary dynamo,’’ Astron. Nachr. 305, 257–264.
Stevenson, D. J., T. Spohn, and G. Schubert, 1983, ‘‘Magnetism and thermal evolution of the terrestrial planets,’’ Icarus 54, 466–489.
Stewartson, K., 1957, ‘‘On almost rigid rotations,’’ J. Fluid Mech. 3, 17–26.
Stewartson, K., 1966, ‘‘On almost rigid rotations, Part 2,’’ J. Fluid Mech. 26, 131–144.
Stieglitz, R., and U. Mu¨ ller, 2000, ‘‘Experimental demonstration of a homogeneous two-scale dynamo,’’ Phys. Fluids (to be published).
St. Pierre, M. G., 1993, in Theory of Solar and Planetary Dynamos, edited by M. R. E. Proctor, P. C. Matthews, and A. M. Rucklidge (Cambridge University, Cambridge, England), pp. 295–302.
St. Pierre, M. G., 1996, ‘‘On the local nature of turbulence in Earth’s outer core,’’ Geophys. Astrophys. Fluid Dyn. 83, 293–306.
Su, W., A. M. Dziewonski, and R. Jeanloz, 1996, ‘‘Planet within a planet: Rotation of the inner core of the Earth,’’ Science 274, 1883–1887.
Taylor, J. B., 1963, ‘‘The magnetohydrodynamics of a rotating ﬂuid and the Earth’s dynamo problem,’’ Proc. R. Soc. London, Ser. A 274, 274–283.
van Dyke, M., 1964, Perturbation Methods in Fluid Mechanics (Academic, New York).
Vanyo, J. P., 1991, ‘‘A geodynamo powered by luni-solar precession,’’ Geophys. Astrophys. Fluid Dyn. 59, 209–234.
Velı´msky´ , J., and C. Matyska, 2000, ‘‘The inﬂuence of adiabatic heating/cooling on magnetohydrodynamic systems,’’ Phys. Earth Planet. Inter. 117, 197–207.
Verhoogen, J., 1961, ‘‘Heat balance of the Earth’s core,’’ Geophys. J. R. Astron. Soc. 4, 276–281.
Vidale, J. E., and P. S. Earle, 2000, ‘‘Fine-scale heterogeneity in the Earth’s inner core,’’ Nature (London) 404, 273–275.
Walker, M. R., and C. F. Barenghi, 1999, ‘‘Nonlinear magnetoconvection and the geostrophic ﬂow,’’ Phys. Earth Planet. Inter. 111, 35–46.
Zatman, S., and J. Bloxham, 1997, ‘‘Torsional oscillations and the magnetic ﬁeld within the Earth’s core,’’ Nature (London) 388, 760–763.

Rev. Mod. Phys., Vol. 72, No. 4, October 2000

