zotero/storage/UCJNFSY8/.zotero-ft-cache

This article was downloaded by: [New York University] On: 07 January 2015, At: 16:58 Publisher: Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK
Geophysical & Astrophysical Fluid Dynamics
Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/ggaf20
The geodynamo, past, present and future
Paul H. Roberts a & Gary A. Glatzmaier b a Institute of Geophysics and Planetary Physics, University of California , Los Angeles, CA, 90095, USA b Institute of Geophysics and Planetary Physics, University of California , Santa Cruz, CA, 95064, USA Published online: 19 Aug 2006.
To cite this article: Paul H. Roberts & Gary A. Glatzmaier (2001) The geodynamo, past, present and future, Geophysical & Astrophysical Fluid Dynamics, 94:1-2, 47-84, DOI: 10.1080/03091920108204131
To link to this article: http://dx.doi.org/10.1080/03091920108204131
PLEASE SCROLL DOWN FOR ARTICLE
Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content.

This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http://www.tandfonline.com/page/terms-and-conditions

Downloaded by [New York University] at 16:58 07 January 2015

Geophys. Astrophys. Huid D y m i c s . Vol. 94, pp. 47-84 Reprints availabledirectly from the publisher Photocopying permitted by License only

02001 OPA (Overseas Publishers Association)N.V.
Published by License under the Gordon and Breach Science
Publishers imprint. Printed in Malaysia.

Downloaded by [New York University] at 16:58 07 January 2015

THE GEODYNAMO, PAST, PRESENT AND FUTURE
PAUL H. ROBERTSa7*and GARY A. GLATZMAIERb
aZnstitute of Geophysics and Planetary Physics, University of California, LQSAngela, CA 90095, USA; bInstitute of Geophysics and Planetary
Physics, University of California, Santa Cruz, CA 95064, USA
(Received 23 February 2000;In final form 6 November 2000)
The geodynamo simulation of Glatzmaier and Roberts (1996, Physica D97,81) is driven by the cooling of the model Earth, which releases latent heat and light components of core fluid at the freezing surface of the inner core as it advances outwards. At some time in the past, the inner core was only a quarter of its present size and at some time in the future it will be twice its present size. The geodynamo operating during those epochs are studied here, the three models (past, present and future) being tied together in an evolutionary sense. The time taken for the models to evolve from past to future depends on the cooling rate, which is controlled by the dynamics of the mantle and is not studied here. All three models generate external fields of comparable strength and all three appear to be close to Taylor states. Unexpectedly, the future model showed considerable variability in time, while the past model does not. Deviations from axisymmetry in the external field increase with inner core radius and the relative predominance of the centered dipole over other multipole components declines.
Keywords: Dynamo theory; Geodynamo; Taylor states
1. INTRODUCTION
Jacobs (1953) proposed that the solid inner core (SIC) was created by the freezing of the fluid outer core (FOC) during the secular cooling of the Earth since its creation. Verhoogen (1961) realized that the latent heat released at the inner core boundary (ICB) during freezing would provide an energy source to power core convection. Braginsky (1963, 1964) pointed out that the freezing would also release light
*Corresponding author. e-mail: roberts@math.ucla.edu
47

48

P. H. ROBERTS AND G . A. GLATZMAIER

Downloaded by [New York University] at 16:58 07 January 2015

components of core fluid, and that these would enhance the buoyancy. Since then there have been many further contributions to the theory, some concerned with global aspects of convection and the overall evolution of the core (e.g., Gubbins et al., 1979; Buffett et al., 1996; Labrosse et al., 1997) and others that developed detailed equations for combined compositional and thermal convection in the core (e.g., Braginsky, 1964; see also Braginsky and Roberts, 1994, 1995, 2000, which are collectively referred to here as ‘BR’).
Glatnaier and Roberts (1996a, here ‘GRl’) used BR as the basis of an evolutionary dynamo model, which was taken further in Glatnaier and Roberts (1996b, 1997, here ‘GR2’ and ‘GR3’), in Roberts and Glatzmaier (2000), in Glatzmaier et al. (1999) and in the present paper. The model creates thermal and compositional buoyancy at the ICB through specifiedcooling at the core-mantleboundary (CMB). In this it differs from other three-dimensional geodynamo simulations, such as those of Glatzmaier and Roberts (1995a,b, 1996c), Kuang and Bloxham (1997), Christensen et al. (1999), Sakuraba and Kono (1999), and Kutmer and Christensen (2000) that are Boussinesq models, steadily driven by specified thermal sources at the ICB and/or in the bulk of the fluid, and also differs in a similar way from the 2;-dimensional model of Sarson and Jones (1999), i.e., a model in which only one of the asymmetric Fourier modes is included in the expansion of the solution. The GR models allow for the compressibility of the Earth, so that the core density increases with depth in close agreement with the Preliminary Reference Earth Model (PREM) of Dziewonski and Anderson (1981).
The present paper is concerned with the geodynamo at a time when the radius ~ I C Bof the SIC was a quarter of what it is today, and also with the geodynamo at an epoch in the remote future when rIcB will be 70% of the radius, rCMB, of the CMB. It will be convenient to specify
these models by the abbreviations ‘S’, ‘G’and ‘L‘,standing for ‘small’,
‘generic’and ‘large’, and also to refer to them as ‘the ancient Earth’, the ‘present day Earth’ (PDE) and the ‘future Earth’ models. The geometry of the S and L models is very different from the G models of G R l - 3, but the governing differential equations and boundary conditions are identical. There is therefore no point in repeating the mathematical formulation, which the reader can find in GRl. Technical details basic to case G have been described in GRl and

THE GEODYNAMO

49

Downloaded by [New York University] at 16:58 07 January 2015

GR2; this case was also one of the 8 cases discussed in Glatzmaier er al. (1999). By now, this case has been integrated for a total of more than a half a million years of simulated time. Cases S and L were started from case G and were integrated for about 30,000 years in case S, and for more than 90,000 years in case L. Appendix A describes how these models differ from case G. We shall now summarize the ten principal assumptions on which the models rest.
Basis of the models:
1. The Earth is cooling on a geologically long timescale (ts, say). As it cools, the inner core grows, releasing latent heat and light component of core fluid which, for simplicity, we assume is a binary iron-rich alloy. These sources, together with internal energy loss, drive the core into convection so weakly that the basic state of hydrostatic equilibrium is scarcely affected, but so strongly that the basic state is well-mixed. The FOC is therefore isentropic and
chemically homogeneous. Its specific entropy, 3, and the mass
fraction of light component, t,depend only on Is. The density p ,
pressurep , gravitational field g, sound speed iis, etc., are by design
close to those of PREM and, like that model, the small changes
created by centrifugal forces are ignored, so that p , p, g, iis, etc., depend only on distance, r, from the geocenter (and on rs). (Basic
state variables carry an overbar as above; contributions to the variables arising from convection are barless, e.g., the total density
+ is p p.)
2. The inner core is passive. By this we mean that it contains no radioactive heat sources (see also point 5 below), and that the deposition of Joule heat, through the electric currents flowing through it, is insufficient to create sub-solidus convection. Its
t ( r ) composition is therefore “frozen in”, being (at a particular
distance r from the geocenter) the composition that it had when that particular layer was earlier added to the surface of the ICB. The diffusivity of heat is much greater than that of composition, and it would be unreasonable to suppose that S ( r ) is, like t ( r ) , frozen into the SIC. As in GRl-3, we suppose for simplicity that
the adiabatic heat flux is continuous at the ICB; see Eq.(12) of
GRl. This approximation is discussed in more detail in Appendix A but, in view of the uncertainties in the parameters mentioned in

50

P.H. ROBERTS AND G . A. GLATZMAIER

Downloaded by [New York University] at 16:58 07 January 2015

point 3 below, it would be hard to justify a more elaborate model. Because of the mobility of the FOC, the temperature T ~ CoBf the ICB is almost uniform over its surface. The SIC is free to turn about the polar axis in response to the magnetic and viscous torques to which the FOC subjects it. 3. There is no exchange of material between mantle and core; ~ c M Bis a constant independent of rs, as is the pressurepCMBand the mean gravitational field gCMBat the top of the core. The total mass of the core is then time-independent, and the total individual masses of the light and heavy components of the core are unchanging, no matter whether they are in solid or liquid phase. The redistribution of mass leads to increasing central condensation
and to secular changes in 3, t,p ( r ) , p ( r ) , g ( r ) and UJr) that are
easily computed in the way described in the Appendix. It may be recalled that ~ ( ~ I c Bto) gether with determines in principle the
mean temperature T(TIcB)of the phase boundary, i.e., the ICB.
Modeling the mean temperature T ( r )is problematical. Because so little is known about the physical properties of core fluid at these great pressures, even the principal alloying element being unknown, several of the parameters involving freezing of the alloy are very uncertain. In the absence of sure knowledge to the contrary, we have supposed that some significantparameters have the same values for all three models, such as the specific heats, the coefficient of volume expansion, the latent heat of crystallization, the rejection factor (which is the fraction of light component released on freezing), the depression of the freezing point by the
light constituent, the fraction, f, of the density jump at the JCB boundary that can be attributed to the compositional discontinuity there, etc. (we may note that the remaining fraction, 1-f, representing contraction of the core fluid as it freezes at the ICB,
together with the contraction associated with the slow cooling of SIC and FOC, results in a secular decrease in rCMB, which we are ignoring; see above.)
Having made these assumptions, the secular decrease of in the FOC is determinate. From our knowledge of PDE, we can uniquely specify the basic state from ~ I C B I. n this sense, our S, G and L models are uniquely linked together; see Appendix A. How

THE GEODYNAMO

51

Downloaded by [New York University] at 16:58 07 January 2015

FICB changes with ts cannot be determined at this stage, so that the
S and L models cannot be dated.
4. The heat flux, F, from the core to the mantle is known. This is a short cut through which we avoid facing the much more complicated problem of solving for the thermal state of mantle and core as a coupled system, obtaining F as part of the resulting
solution. The larger part ( F ) of F is the heat flux down the
adiabatic gradient defined by the basic state of the core; it is
independent of colatitude, 8,and longitude, 4, on the CMB. The
remainder, which may be positive or negative, is the convective heat flux. Because the core fluid is so mobile, TCMBis almost uniform over the entire CMB but, because of convection in the
mantle, F will in general depend on 8, and 4; see GR3. Because tS
is the timescale of mantle convection, it is also the timescale on which F vanes, in the reference frame locked to the mantle. For simplicity it is assumed here that F is uniform over the CMB so that F = Q c ~ ~ / 4 7 r R & w~ ,here QcMB is the net heat flux to the mantle; models where this assumption was not made are reported in GR3 and in Glatzmaier et al. (1999). 5. There are no radioactive sources in the FOC; F is entirely the result of the cooling of the Earth. The mantle acts as a valve controlling the core. The greater the F, the more rapidly the core cools, the faster ICB advances, and the greater are the thermal and compositional sources on the ICB driving convection. It should be
emphasized that these sources are determined by F and the
convective motions themselves; they are not specified separately.
Moreover at the ICB they inevitably depend on f3 and 4, since the
ICB advances faster where the cold convection currents descend onto it than where the hot convection currents leave it; the thermal and compositional sources on the ICB are enhanced or suppressed in proportion to ~ I C B .This creates a modest topography on the ICB; see Glatzmaier and Roberts (1998). Even in the absence of the SIC, cooling at the CMB creates a destabilizing entropy gradient that can power a convective dynamo. Because the buoyancy sources are then distributed less favorably, throughout the FOC rather than at the bottom of the FOC, the mechanism is inherently less efficient than the convection mechanism operating

52

P. H. ROBERTS A N D G. A. GLATZMAIER

Downloaded by [New York University] at 16:58 07 January 2015

in the PDE and the future Earth, and we found it necessary to

assume a larger Q c M B for the S model than for the G and L

models. The specified F determines i ‘ l c ~w, hich in principle allows

rICB to be determined as a function of rs, and the S and L models

to be dated. The assumption that there is no radioactivity in

the core is consistent with the current belief of many geochemists

(e.g., Stacey, 1992), but tends to require that the inner core be

only 1 or 2 billion years old. As BR and Braginsky (1997) point

out, this difficulty would be less severe if small amounts of

radioactivity were present. Chabot and Drake (1999) recently esti-

mated that 40Kdissolved in the core currently generates 10” W of

heat.

6. The angular velocity of the Earth depends on tS. Assuming that

the SIC grows at a volumetrically constant rate and that the age of

the inner core, T~~~(say), is approximately 1.3 x lo9yr (GR3), the

angular velocity, R, of the Earth at a time when the inner core

radius was rICB/4was 9.2 x 10-’rad/s, where we have assumed

that the tidal deceleration of the Earth is = -5.2 x

rad/s2

(Munk and Macdonald, 1960; Lambeck, 1980). We take

R=9.2 x lO-’rad/s for model S,7.29 x lO-’rad/s for model G

and 3.64 x lO-’rad/s for model L. The distribution of shallow

seas between now and about 8.4 x lo9 years in the future when

(again assuming constant volumetric growth) ~ I C B=O . ~ ~ C M Bi,s

completely unknown, as is therefore the tidal deceleration. The

reduction of 0 to half its present value is therefore completely

arbitrary. The character of the solutions would not be affected

drastically until R was much smaller than this.

7. Convection alone powers the geodynamo. The dynamical effects

of the luni-solar precession are ignored.

8. Core convection operates on a timescale (I,-, say) short compared

with tS. In simulating core convection and the geodynamo, the

basic states of mantle and core are “frozen in time”. (In math-

ematical parlance, we employ two-timescale methods.) If it were

economically possible to simulate the core over periods compar-

able with ts, it would be necessary to update the basic state

continuously, treating ts as a passive parameter. If, as in the

present paper, we wish to understand the geodynamo at a remote

epoch in the past or future, we need only specify the basic state,

THE GEODYNAMO

53

Downloaded by [New York University] at 16:58 07 January 2015

together with Fat that epoch, and to regard these as constants, as

indeed they are on the tc timescale.

9. It is impossible to resolve numerically the small-scale fields and

flows expected to be significant in the core. Following Braginsky

and Meytlis (1990) we represent their effect in transporting large-

scale entropy, chemical composition and momentum by turbulent

fluxes, as is commonly done in so-called 'local' models of

turbulence. Turbulence also adds significantly to the energy

dissipation, and this is recognized in the energy budget of our

models. Although, as Braginsky and Meytlis point out, core

turbulence is probably highly anisotropic, we suppose that the

turbulent thermal diffusivity, the turbulent compositional diffu-

sion coefficient, and the turbulent viscosity are scalars and not

tensors. Braginsky and Meytlis suggest that they should be of the

same order of magnitude as the magnetic diffusivity, 7j = 2 m2s-',

and we adopt this value for the turbulent diffusivity of heat and

composition. For numerical reasons (since we specify the Earth's

dimensions and rotation rate), the kinematic viscosity had to be

further increased, to F = 1450m2s-', and hyper-diffusivity

(Glatzmaier and Roberts, 1995a) had to be added to all diffusion

coefficients to damp the smaller scales. Our Prandtl number,

P = F/E,where K is thermal diffusivity, is therefore unrealistically

large. Although our angular velocities are Earthlike (see point 6),

our Ekman number E = F/2n(?CMB - FICB)*for model G is at

least 2 x

which is 3 orders of magnitude greater than for the

real Earth, even if we adopt F = 2 m2s-' as a turbulent viscosity

rather than its probable molecular value (10-6m2s-1).

In Section 3 we shall make some comparisons between our

results and those of Kutzner and Christensen (2000), which studies

Boussinesq dynamos driven by thermal buoyancy alone in the

absence of electrical conduction in the SIC. Their models, which

are based on that of Glatzmaier (1984), also are fully three

dimensional (3D). By restricting themselves to much larger

viscosities ( E - 3 x lo-"), they were able to avoid using hyperdif-

fusion and to study convection at P of unit order. The Rayleigh

number they use (a nondimensional measure of the buoyant

driving in the core) is also smaller than ours and is very probably

also much smaller than that of the FOC.

Downloaded by [New York University] at 16:58 07 January 2015

54

P. H. ROBERTS A N D G. A. GLATZMAIER

10. Seismicwaves are unimportant to the dynamo and are filtered out by replacing the mass conservation equation by the anelastic approximation, which improves on the Boussinesq models by recognizing that fluid density depends on p and T.
Paleomagnetism indicates that the Earth has been magnetic for at least 3.6 x 109yr (e.g., see Merrill et al., 1996). Although it has reversed polarity often and irregularly, it seems to have changed remarkably little in other respects. For example, Kono and Tanaka (1995) estimate that its dipole moment has been within a factor of 3 of its current value over most of geological time. This is surprising considering how greatly the internal structure of the Earth has changed, especially in respect of the SIC whose age, TSIC, is less, and possibly much less, than the age of the Earth. Stevenson et al. (1983) estimated that ~ 1 ~ 5 2x .1309yr, Labrosse et al. (1997) that ~ ~ ~ ~x 51019yr., B7uffett et al. (1996) that ~ ~ ~x 1~09yr~and2 . 8 GR3 that 7sIcc 1.3 x lo9yr.
Hollerbach and Jones (1993) have argued that, despite its small size, the SIC is at present influential in core magnetohydrodynamics (MHD). They argued that, since its magnetic diffusion time T& (x2400yr) is long compared with tc, the SIC of PDE suppresses fluctuations in the magnetic field B in the FOC and reduces the frequency of polarity reversals.' Our early simulations (Glatzmaier and Roberts, 1995a, b, 1996c)supported their speculation, in the sense that the magnetic field generated in the outer core was more stable and dipole-dominated when the inner core was conducting, as compared with an electrically insulating inner core of the same size. For our
.iIc much smaller ancient inner core x 150 yr, making it seem unlikely
that the MHD state of the FOC is significantly smoothed by such a
.ilC small SIC. For our future core, x lo4 yr, making it seem plausible
that SIC smoothing action is great. It has long been speculated that core flows and fields were very different in the past when the SIC was small or absent (e.g., Braginsky, 1964,Stevenson et al., 1983). Buoyant driving from the bottom of the FOC would then be weak or absent,

'We here took T& = &-B/7?r), just as though the SIC were surrounded by insulator, and assumed that, as in our models, r) is the same in both inner and outer cores; in reality, the magnetic diffusivity in the inner core is probably only about 75% of that in the outer core; see BR.

THE GEODYNAMO

55

and we expect (see point 5 above) that the convective efficiency would be less than for PDE. We may also expect that the MHD dynamo will function differently in the future.2 One objective of the present paper is to examine the remaining questions by simulating the geodynamo at a time when the radius of the SIC was only a quarter of its present size, and also when it will have doubled. We hope that, by contrasting the
S, G and L solutions reported below, an increased understanding will
emerge of how the structure and temporal behaviors of the geomagnetic field depend on the gross geometry of the core, an objective similar to that of Sakuraba and Kono (1999). One of the more interesting results of the present investigation is that the magnetic field is more stable for our model of the ancient Earth than for our present Earth model, which in turn is more stable than the future Earth model.

Downloaded by [New York University] at 16:58 07 January 2015

2. RESULTS
We suppose that the total heat flux, QcMB,from the future core to the
mantle is the same as we assumed in GR1 and GR2 for the PDE,
namely 7.2TW; we suppose that QCMB=1O.OTW for the ancient
core.3 The total heat flux is divided between the adiabatic heat flux
p,",, (associated with molecular thermal diffusion down the
temperature gradient) and the superadiabatic heat flux @gB
(associated with turbulent thermal diffusion down the entropy gradient)B;M:@ is estimated in Appendix A. These are listed in Table I along with other quantities that are of interest later.
We now show a number of typical snapshots of the S, G and L
solutions starting (Fig. 1) with the horizontal averages 3 and of S
and E. (Recall here - see Section 1 - that these refer to the convective
r, contributions alone; the volume averages of entropy and composition,
3 and representing the slowly evolving reference state, have been
'One possibility that we do not examine is that, as the SIC grows beyond PDE and the mass fraction of light constituent in the FOC increases, the composition of the FOC may become eutectic. If so, compositional buoyancy will be lost, and only thermal sources will be available to drive core convection. We shall assume that compositional driving is present at all times up to, and including, model L.
'When we assumed Q c ~ e = 7 . 2 T Wfor the S core, it appeared, possibly for the reasons given at the end of Section 1, that the magnetic field B was dying out. To be sure of obtaining a working dynamo, we therefore increased Q ~ MtBo IOTW.

56 Quantity

P. H. ROBERTS A N D G . A. GLATZMAIER

TABLE I ProDerties of solutions

Small

Generic

Larne

305.37 9.121 10.0 5.89 4.11 0.024 0.18 0.8 1.4 2.1 1.8 0.5
64.
6.44 0.47 60. 0.13 3.5 0.47
1 .o
1.4

1221.5 9.044 7.2 5.38 1.82 0.76 13. 0.7 1.1 5.6 9.9 2.1
225. 3.25 2.76
143. 0.12 1.5 9.2 1.7 3.9

2443.0 7.298 1.2 4.43 2.17 4.38 50. 1.4 0.6 2.1 2.5
1 .o
18. 5.09 12.3 63. 0.031 1.2 2.5 0.7 0.6

Unit

~

~~

km

kg.m2

Tw

Tw

Tw

1030J 10~~~

103kg/s 10 ~2 kg/m3

IO-~~/S

"/yr 10 'm/s
~
G lozzA m2

G A/m2 10 ' G m/s
~
10 - l o Gm/s2

Downloaded by [New York University] at 16:58 07 January 2015

subtracted.) The surface of radius r ( > rICB) is required to transmit
outwards the entropy generated within it by Joule heating (for both the resolved and unresolved scales), conduction of heat down the adiabat, and viscous regeneration of heat. This internally generated entropy increases with r as must, therefore, the flux of entropy. This leads to the increasingly negative gradient o f ? seen in the right-hand
portion of each curve. The initial sharp decrease in ? is due to the
entropy source on the ICB. This is comparatively small for the ancient Earth; because of its small area, the ICB is a relatively ineffective source of thermal buoyancy for the S core. The slope of-the curve at
the CMB is proportional to @gB,but the slope of the curve is zero
at the CMB because there is no mass flux into the mantle, by assumption 3 of Section 1. The slope of the curve at the ICB is larger for the PDE than for the future Earth. Nevertheless, because its area is 4 times greater, compositional buoyancy is more significant for the future Earth. The ICB of the ancient core is 16 times smaller in area than that of the PDE, and compositional buoyancy is comparatively unimportant for the S model.
Since we employ the anelastic approximation (see point 10 of
Sect. l), V . (pV)= 0, and therefore a streamfunction 1c, exists for the

THE GEODYNAMO

51

Downloaded by [New York University] at 16:58 07 January 2015

0

CMB
Radius

Light Constituent

0

Radlus

CMB

FIGURE 1 The horizontal means of the convective contributions to the specific entropy and the mass fraction of light constituent from snapshots of simulationsof the geodynamo in the ancient, present day and future core. The overall maxima are
respectively 2.25 x 10-4Jkg-'K"-' and 5.05 x lo-*.

mean mass flux, ( p v ~in)m, eridian planes. (Asin Tab. I and below, we
denote the axisymmetricpart of a fieldfby (f).) The top three panels
of Figure 2 show equi - 1c, contours for the S , G and L snapshots. More
precisely, we may call these curves 'streaklines', since they are time dependent and show only the instantaneous direction of the flux. The fluid flow is counterclockwise round the contour when shown as full curves and clockwise when displayed as dashed curves. The bottom three panels of Figure 2 show the contours of the angular velocity
(relative to the rotating frame of reference) (V+)/s,where s=r cos 8 is distance from the rotation axis and (r, 8, 4) are polar coordinates. The
zonal flow is eastward when full curves are used and westward where dashed curves are shown; it is remarkably symmetric with respect to the equatorial plane, particularly in case S. The zonal flow, and the

58

P. H. ROBERTS AND G. A. GLATZMAIER

Meridional Mass Circulation

Downloaded by [New York University] at 16:58 07 January 2015

Differential Rotation

FIGURE 2 Streaklines of the mean mass flux in meridian planes and contours of the

' mean angular velocity for the three snapshots. The contour separations of the former are

2 x lo9kg s ~

for the S core, and 1.4x 10" kg s I for the G and L cores; for the latter,

the contour intervals are approximately 0.6"yr I in all 3 cases. The maxima of the

meridional mass flux are about 0.4,2.9 and 3.5 x

10" kgs ~

I; the maxima of I(VM)lare

2.1, 5.6 and 2.1 x 10-4ms I; the maxima of the angular velocity are 2.9, 9.9 and

2.5"yr-',andthemaximaof~(V,)~are5.1,21 and 10 x 1 0 - 4 m s ' , f o r t h e S , G a n d L

snapshots respectively.

largest gradients of the meridional flow, are confined in case G to the (TC), i.e., the cylinder parallel to the rotation axis and touching the inner core at its equator; the S core is too small to display any such effect. The zonal flow in cases G and L is eastward near and on the

Downloaded by [New York University] at 16:58 07 January 2015

THE GEODYNAMO

59

ICB. The maximum values of ( V M )(,V + )and ( V+)/s,together with the maxima of other fields considered below, are listed in Table I.
The SIC responds to the magnetic stresses (and to some extent the viscous stresses) to which the FOC subjects it. The resulting magnitude and sense of rotation of the SIC was studied in depth in GR2. Sufficeit to say here that, in obedience to Lenz's law, the SIC tends to turn in
' the same direction as the adjacent overlying fluid, which is prograde in
cases G and L. The angular velocity of the SIC is S2s1c=2.4" yr- for the G snapshot, but is less for the L snapshot (R~*c=0.55"yr-'), probably because the zonal wind shear at the ICB is less in case L. The
zonal shear is also quite small in the S model and probably does
comparatively little to generate zonal field by the w-effect. The magnitude of the angular velocity for the S snapshot is correspond-
', ingly small: Rsrc= -0.43" yr- where the minus sign indicates
westward motion. The moment of inertia of the SIC, ZSIC,scales as
GI,, but the magnetic torque on it scales as riIc (for the same magnetic
field strength). This means that Rsrc varies very little for the L core, but vanes markedly for the S core. The time dependencies of Rsrc in
the S and L cases are contrasted in Figure 3. Taken over the periods
employed in that figure, the mean and standard deviation of RsIc for
the S core are -0.053"yr-' and 0.55"yr-', and for the L core are
0.74" yr-' and 0.044"yr- I. According to GR2, the corresponding values for model G are 2.6" yr-' and 0.4"yr-I.
The buoyancy generated at the ICBs of the all three models tends to maintain meridional circulationsinside the TC that are consistent with

Small inner core

Large inner core

time (Yr)

o so0 loon 1500 2000 2500 time Qr)

' FIGURE 3 The angularvelocity,Rs~co,f the solid inner core in yr - as a functionof
time for the S and L models.

Downloaded by [New York University] at 16:58 07 January 2015

60

P. H. ROBERTS AND G. A. GLATZMAIER

the zonal flows shown. That is, for cases G and L, fluid near the z-axis is less dense than that on the TC, at the same radius (Fig. 4). Therefore, on average, fluid rises near the z-axis, causing fluid to move toward the z-axis when near the ICB and away from it when near the CMB. Coriolis forces resulting from the circulations maintain eastward zonal flow near the ICB and westward near the CMB. The gradient normal to the z-axis of the longitudinally-averaged density perturbation (Fig. 4)inside the TC for case S is essentially the opposite of those for the other two cases. As a result, for case S, fluid tends to sink along the z-axis; the additional drift toward the axis as the fluid sinks maintains the eastward zonal flow there. Because of the small inner core, the circulation along the axis for case S is also part of the circulation along the TC, which is similar to the separate circulation cells near the TC for cases G and L.
Figure 5 displays snapshots of the axisymmetric magnetic field in the three cases. The meridional field lines are shown as full curves when directed in the counterclockwise sense and as dashed curves

Density Perturbation

FIGURE 4 Contours of the convective contributions to the density perturbation for the three snapshots. Solid contours represent more dense fluid relative to the mean; broken contours represent less dense fluid. The contour separations are 1.0, 0.7 and 0.4 x 10W4kgm- for the S , G and L cases, respectively.

THE GEODYNAMO

61

Meridional Fields

Downloaded by [New York University] at 16:58 07 January 2015

Zonal Fields
FIGURE 5 Lines of force of the axisymmetric meridional field and contours of the axisymmetriczonal field for the three snapshots. The dux separations for the former are lOOTkm2 for the S and L cases, and 80Tkm’ for the G model; for the latter they are 0.7mT for the S model and 1mT for the G and L models. The maxima of I(BM)Iare 6.4, 23 and 18mT, and of I(Em)l are 6.0, 14.3 and 6.33mT. for S, G and L snapshots respectively.
when clockwise. The zonal field is eastward where the contours are full curves and westward otherwise; it is fairly antisymmetricwith respect to the equatorial plane. The zonal field, and the largest gradients of the meridional field, are confined to the TC in case G, but the S core is too small to display any such effect. It is striking how little the fields spread into the SIC in case L. This is not totally surprising. Field is not generated in the SIC and any currents that diffuse into its interior

Downloaded by [New York University] at 16:58 07 January 2015

62

P. H . ROBERTS AND G. A. GLATZMAIER

represent an ohmic loss of magnetic energy detrimental to the functioning of the dynamo.
The argument of Hollerbach and Jones (1993), that the SIC will act as a moderator of the field that discourages reversals, appears to have less force for the L core than for the G core. There are apparently two related reasons for this. First, because the depth of the FOC is less, the preferred length scales are somewhat shorter (even with our hyperdiffusion). Second, time scales, such as the convective overturning time and the effective magnetic diffusion time, are reduced. As a result the L dynamo is quite variable, reversing twice, at times separated by about 15,000yrs, during the period of approximately 90,000yrs of simulated time over which it was integrated. (A further illustration of its variability is shown below.) These reductions in length and time scales means that the magnetic field does not penetrate as far into the SIC for the L model. The average field strength at the ICB is only 34 G and the average angle made by B with the horizontal is only 12" for the L snapshot; the corresponding values for the G snapshot are 56 G and 16", indicating a greater penetration of B into the SIC.
The axisymmetric part of the current density, (4,where J = IJI, is
shown for the 3 snapshots in Figure 6, for the fluid core only. It is striking, though not surprising in view of Figure 5, how strongly the

Electric Current Density

' FIGURE 6 Contours of mean current density in the fluid outer core for the three
snapshots. The contour intervals are 1.3 x 10 A m for the S and G models but twice ~
that for the L model. The maxima of 1(1J1)1are 0.13.0.12 and 0.031A m - 2 , for S, G and L snapshots respectively.

Downloaded by [New York University] at 16:58 07 January 2015

THE GEODYNAMO

63

currents are confined to the TC in case G. The power spectra at the CMB for the snapshots are shown in Figure 7. They are nearly flat at the CMB for 1=2 -5 , but then decrease, principally because of our inclusion of hyperdiffusivity. When hyperdiffusion is not present, as in the model of Roberts and Glatzmaier (2000), the spectra remain almost flat to significantly higher harmonic degree 1, before they too decline. It is clear that the smaller the radius of the inner core the more the field is dominated by the central dipole. It is also true that that dipole is increasingly axial, its angle of tilt from the geographical axis being 0.4”, 3.9” and 17.4” respectively for the S, G and L snapshots. The mean of the dipole moment, m, and the tilt (averaged over a few thousand years at the end of the three simulations) are shown in Table I. Figure 8 illustrates the variability of the L model through plots of m and the dipole tilt as functions of time. The m =0 and rn = 1 quadrupolar components of the field can be allowed for by displacing the centered dipole to a magnetic center, to create a “displaced” or “eccentric” dipole (e.g., see Connerney, 1993). The magnetic center of the L snapshot is situated at only 0.1 rCMB from the geocenter, much as it is for the real Earth today; in our S and G snapshots however it is situated at 0.02 rCMB from the geocenter.
How do the dynamos function electromagnetically? Cowling’s theorem shows that the &component of the electromotive force, &= (V’xB’), produced by the asymmetric components V’ =V - (V) and B’ =B - (B) of velocity and field, is essential in maintaining a field outside the core that has an axisymmetric part. Lacking

the dynamo can at best maintain an asymmetric B outside the core. Contours of €+ are shown in Figure 9.
Simple axisymmetric dynamos have often been constructed in which &@ is represented by a constant or spatially-varying multiple, a,of (B,) or in which & is taken to be a(B). These are called “a-effect models” and are usually motivated by analyses of microscale
induction, i.e., the creation of (BM)by the part of & arising from the microscale parts of V and B. These suggest that there is a connection
between (Y and H , the simplest proposal being that a / H is approximately constant and negative, where H is the local helicity of the turbulence. Although there are many astrophysical contexts in

0
W Jm-?, . -2 - 0

-3.

-4-

ooooooooo
0

a

-5 -6 .

00
00 00
OO

-7 -

00-

-8

Downloaded by [New York University] at 16:58 07 January 2015

0
W Jrn-?, .

-2 .
0

-3.

-4 .

o

o

~

o~o
000

~

b

-5 -

00 00

-6

00 00

-7 -

00.

-0

0
W Jm-?l

-2-0

-3-

ooooo ooooo

-4 -

00

C

0

-5.

00

0

-6

OO

-7 -

OO-

-8-

FIGURE 7 Power spectra at the core-mantle boundary for (a) the small core, (b) the
present core, and (c) the large core in units of J m '.

THE GEODYNAMO

65

Moment (Am2)
6r

Tilt (degrees)
30 r

Downloaded by [New York University] at 16:58 07 January 2015

4.51 0 lo00 zoo0 3000 4000 5000
time (yr)

0 loo0 2000 3000 4ooo 5ow
time (yr)

FIGURE 8 (a) The dipole moment (in Am’), and (b) the angle of tilt (in degrees) between the geographical axis and the axis of the centered dipole for the L model, as functions of time t (in yr).

which microscale induction is believed to be significant, the Earth, because of the large magnetic diffusivity of its core, does not appear to be one of them (see, e.g., Braginsky and Meytlis, 1990). Nevertheless one might expect that, if instead of using the unknown microscale fields,we used the asymmetric(though macroscale)parts of V and B to evaluate E and H, we might see some vestiges of the assumed connection between a and H. We therefore define a and H by

H = (V’* V X V‘).

(2.3)

Direct investigation of the conjecture that a and H a r e closely related encounters a difficulty: a zero of &@ does not generally coincide with a zero of (B&,so that (2.2) leads to infinitiesin a. We therefore compare the &@ with H(B+)in Figure 9. The helicity tends to be large only near the SIC so that, particularly in case L, most of the (BJ shown in Figure 5 is eliminated when multiplied by H. According to Figure 9, there is no very direct connection between a and H, though both are approximately antisymmetric with respect to the equatorial plane. Previous dynamo simulations by Olson ef af. (1999) found a good correlation between H(B,) and &&. Their simulations involved, however, weaker convective driving and less rotational dominance, and so dynamo action mainly occurred outside the TC, whereas in our

66

P. H.ROBERTS AND G. A. GLATZMAIER

Downloaded by [New York University] at 16:58 07 January 2015

' FIGURE 9 Contours of E , and of the product of helicity, H, and the mean zonal field,
(B,), for the three snapshots. The contour intervals for the former are 2 x 10 'Tm sfor the S model and 6 x 10 - T m s - for the G and L models; for the latter they are 1.5 x 1 0 - ' 6 T m s C Zfor the S model and 5 x 1 0 - t 6 T m s - 2 for the G and L models.
', The maxima of lEml are about 3.4, 1.5 and 1.2 x 1 0 - 7 T m s - ' and those of IH(B,+)Iare
about 0.36, 1.3 and 1.8 x 10 I4Tm s - for S, G and L snapshots respectively.
case it arises mostly in the vicinity of the TC. We conclude that the traditional interpretation of the a-effect in terms of the local helicity may only be appropriate for dynamos that are gentler than ours. These conclusions are not much affected if we replace H by
HtOt= (V.V x V), which augments H = (V' .V X V') by the helicity

Downloaded by [New York University] at 16:58 07 January 2015

THE GEODYNAMO

61

of the mean flow R = (V) .V x (V). The structure of the contours
shown in the 3 lower panels of Figure 9 is not substantially altered, although the ratio of the maximum of HI,, to that of H is about 13, 7
and 1.4 for the S , G and L snapshots. This suggests that R has a similar structure to H , but is much larger for the S core though smaller
for the L core. The zonal field, (B&, is generated from (BM), partially by the
differential rotation w, a process often called “the w-effect”, and partially by the meridional components €M of E. The latter, when represented by a ( B M )(see above), is regarded a second (though less essential) part of the 0-effect. A question arises, “Which of these two sources of ( B 4 ) is the more important?” If the former dominates, the model is said to be of “aw-type”; if the latter of “a2-type”; if both are significantit is of “a2w-type”.(In each case, one a refers to the process in which (B,,.,) arises from (B4).)Are the dynamos described here of a2-,aw-or a2w-type?This question concerns the relative importance of

X, = [v x (V‘ x B’)Jg, and X, = [v x ((V)x (B))lg (2.4)

in generating (B4). Our integrations suggest that all 3 models are
basically of a2w-type.Although the maximum of X, exceeds that of X, in the S core, X, is mostly confined to the immediate vicinity of the
ICB, and it appears that the a-effect is slightly more effective in
creating B+ than the w-effect is. In cases G and L, both X, and X, are
confined to the vicinity and interior of the TC, and there seems to be no clear preponderance of one generation mechanism over the other.
Concerningthe dynamical functioning of the dynamos, we may ask, ‘How Taylor-like are they?’, by which we mean ‘How small is the Taylor integral ‘T(s)?’, where

(Taylor, 1963). The integral is performed over every geostrophic cylinder C(s);this is a cylinder of radius s, coaxial with the rotation axis, and of finite length, being confined to the FOC. Taylor’s result follows from the assumption that viscosity has no significant dynamical effect, even in boundary layers, a reasonable (though not

Downloaded by [New York University] at 16:58 07 January 2015

68

P. H . ROBERTS AND G. A. GLATZMAIER

necessarily correct) assumption when the Ekman number, E, is small,
as it is here. Taylor's argument applies outside the TC (s > ~ICB).For
the fluid interior of the TC, we define two analogous integrals over
cylinders, Cds) and C&), to the North and South of the SIC:

(2.6N)

(2.6s)

In Figure 10 we show

It is seen that C is relatively small outside the TC (s > rlcB), except near
s= rCMB. Inside the TC (s < r1cB)there is tendency for 7 ~ ( s )and
7&)to run parallel to each other, suggesting that each spherical cap on the N and S cylinders are exerting roughly the same torques on the SIC. For a further recent discussion of the role of the Taylor constraint, see Jones and Roberts (1999).
Because the area of the ICB is small for the ancient core, the compositional and thermal sources on the ICB are less effective in driving motions than for the PDE's core and the loss of internal energy from the bulk of the FOC plays a much greater role in driving the

I0.5

1
5
0.5

II0.5

'0

0.5

1 '0

0.5

1

0

0.5

1

SkIB

%B

dLB

FIGURE 10 The function <(s) plotted as a function of s/rCMB for the three snapshots. Inside the tangent cylinder (s < rice), the integral (2.6N)is shown by a full line and the
integral (2.65) by a dashed line.

THE GEODYNAMO

69

motions. Such matters can be quantitatively assessed through com-
parison of the individual contributions dCMe",,,B,and gMmaB de
by gravitational settling, internal energy and specific heat to Qt-MB. This matter is taken up in Appendix A where it is argued that

Downloaded by [New York University] at 16:58 07 January 2015

and q is evaluated;here Q$, is the adiabatic flux from the ICB into the FOC. In Table I above, q is listed, as are the implied values of TICB/&B and MSIC,which is the net mass flux from the ICB. These give some idea about the rate at which the system is evolving at the three epochs considered.
3. CONCLUSIONS
Although the expense of MHD integrations in three dimensions has precluded a thorough survey of parameter space, the integrations presented here have revealed a number of interesting effects that may prove to be typical. The salient differences in the solutions for the small (S), generic (G) and large (L) core models may be summarized in the following way:
(a) Dipole dominance of the external magnetic field depends strongly on rICB;see Figure 7. This statement means more than the usual observation that the dipole field decreases with distance from a planet more slowly than that of the other multipoles, so that the greater the dominance of the dipole at the surface r =rp of the planet, the deeper the source of field, i.e., the smaller rCMB/rP.Nor are we referring to the fact that, apart from the dipolar term (f= l), the power spectrum is approximately linear in all three cases, a fact that can be used to provide an estimate of rCMB/rP,as suggested by Hide (1978); see also Glatzmaier and Roberts (1996~).The basis of Hide's idea is that, when extrapolated to the core surface, the energy density contained in the harmonic of
degree f should be the same, independent of 1. Except for the dipole
term, which is anomalously large, this seems to be approximately true in the case of the Earth, though a theoretical explanation is

70

P.H.ROBERTS A N D G.A. GLATZMAIER

Downloaded by [New York University] at 16:58 07 January 2015

lacking, and it is also unclear why the dipole should be exceptional. The simulation of Glatmaier and Roberts (1996~)tended to confirm both the flatness of the spectrum for I > 1 and the prominence of the I = 1 term at the core-mantleboundary of the present day Earth. The results reported here for the small and large cores also show the same pattern but, when we state above that the dipole dominance of the external magnetic field depends strongly on rICB, we are referring to the relative prominance of 1= 1 over the other harmonics at the CMB, and this appears to decrease with increasing inner core ratio rfCB/rCMB. As the fields of other planets and satellites in the solar system become more fully explored, this may provide an indication of the presence (and perhaps radii) of the solid inner cores that they possess. Other factors may, however, dash this hope. For example,Kutmer and Christensen (2000)find that dipole dominance is more evident when the buoyancy sources driving convection are mainly at the ICB than when they are distributed throughout the FOC. Dipole dominance may also depend on the
much smaller Prandtl number P that they assume. Grote et al. (2000) suggest that dipole dominance may be a property of large Prandtl number systems such as ours, for which P x 720. Clearly
more work needs to be done at higher numerical resolution and lower viscosity to clarify this matter. (b) The axisymmetry of the external field increases with decreasing IICB/TCMB. For example, we have seen in Section 3 that the mean dipole tilts are OS", 2.8" and 12.3" respectively for the S, G and L models. Unlike (a), this conclusion is unaffected by uncertainty in rCMB/rP. The tilt of Saturn's (present) dipole axis is less than 1"; see Connerney (1993). The giant planets have, in comparison with the Earth, large internal density stratifications and have electrical conductivities that depend strongly on r; moreover, their MHD states satisfy different boundary conditions. The relevance of our models is therefore doubtful. If nevertheless we speculate, we might suggest that the small tilt in Saturn's field indicates that its dynamo operates throughout its deep fluid core. If this is the case, (a) above suggests that the dipole field is also anomalously strong, and that the assumption that the quadrupole and dipole energies are comparable when extrapolated downwards to the surface where effective dynamo action starts will lead to an underestimate

THE GEODYNAMO

71

Downloaded by [New York University] at 16:58 07 January 2015

of Saturn’s core radius. The high tilt of the magnetic axes of Uranus and Neptune and the large offsets of their magnetic centers (see Connerney, 1993) suggests that both of these planets have relatively large inner cores where dynamo action is weak or absent. (c) The ancient SIC has a much smaller effect on core MHD than that of the G and L models. There is much greater kinetic and magnetic
activity outside the tangent cylinder of the S model.
(d) In view of the stabilizing effect of the inner core, postulated by Hollerbach and Jones (1993) and previously confirmed by ourselves for G models, it may be surprising that the magnetic
field was less irregular for the S model than for the G model and
less for the G model than for the L model. Our earlier results confirmed that an electrically conducting inner core (having the dimensions of the present day Earth) stabilizes the solutions as compared with a model where the inner core is electrically insulating but in other respects is the same. In contrast, we are here comparing models all of which have the same inner core conductivity, and we have found that the smaller the inner core radius, the more stable these models are. This confirms the conclusion of Morrison and Fearn (2000) from studies of 24 D geodynamo models, i.e., models in which only one of the asymmetric Fourier modes are included in the expansion of the solution. It is however at variance with the conclusions of Sakuraba and Kono (1999), who studied two models, identical to each other except that one had an inner core of present day size and the other had no inner core whatever. They found that the latter model is less stable than the former. It is not clear why this finding is opposite to our own, but it may be pertinent to remark that their models are Boussinesq, employ a kinematic viscosity and thermal diffusivity that are 20 times greater than the magnetic diffusivity (whereas, for the Earth, 7 is by far the largest of the three diffusivities),and ignores, for their model possessing a SIC, buoyancy sources on the surface of the ICB. (e) The future core model is rather unsteady and prone to reversals. This may be attributed to lack of communication between the northern and southern hemispheres of the FOC, which tends to make these regenerate as independent systems. Other symptoms of

Downloaded by [New York University] at 16:58 07 January 2015

12

P. H . ROBERTS A N D G. A. GLATZMAIER

this communication difficulty is the greater displacement of the magnetic center from the geocenter, and the greater inclination of the centered dipole noted in (b) above.
(f) The moments of inertia of the SIC of the S,G, and L models are
GCB, proportional to and are therefore approximately in the ratios
:1 :30. The products of surface area and torque arm are
proportional to &,, and are in the more moderate ratios
0.016 :1 :8. The angular velocity, Rsrc, of the S inner core is therefore more responsive to the fluctuating torque arising from the varying state of the surrounding FOC; see GR2. Consequently, it is not surprising that Rsrc is quite variable for the S model despite the relative quiescence of its MHD state. The magnitude of RsIc is also small, which is consistent with the rather small zonal shears that arise in the FOC for this case, and there is no marked tendency for eastward rotation to predominate over westward; see Figure 3. In contrast, the solid L core rotates fairly steadily in the prograde direction. Its greater inertia damps fluctuations in Rsrc despite the greater intrinsic variation of the MHD state of the FOC; see (e). The somewhat smaller (Rsrcl, as compared with the G model, seems to be due to a smaller zonal wind at the ICB and the concomitant reduction in the torque exerted by the FOC across the ICB; see Figure 2. (g) Jault (1996) found from studies of intermediate dynamo models (i.e., axisymmetric dynamic models) that dynamo action is impeded by a large inner core. This effect is not as evident in our integrations. When Figure 7(c) is compared with Figures 7(a) and 7(b), it appears that the field energy produced by the L core is slightly greater than for the G and S cores. For a small inner core, thermal buoyancy is mainly provided volumetrically, through the cooling of the bulk fluid; for a large core, thermal and compositional buoyancy sources are more effectively placed to drive motions, namely at the bottom of the fluid core where heat and light constituents are released; see Figure 1.

Acknowledgments
We are grateful to U. Christensen and to a second but anonymous referee for their constructive criticisms. GAG was partially supported

THE GEODYNAMO

73

by NSF EAR99-02969, PHR by NSF EAR97-25627, and both of us by the NASA HPCC/ESS Grand Challenge Program and the IGPP at the Los Alamos National Laboratory. GAG was also supported by the University of California Research Partnership Initiative. Computing resources were provided by the San Diego Supercomputing Center, the National Center for Supercomupting Applications, the Texas Advanced Computing Center, and the Goddard Space Flight Center. We thank Maureen Roberts for computational assistance.

Downloaded by [New York University] at 16:58 07 January 2015

References
Braginsky, S. I., “Structure of the F layer and reasons for convection in the Earth’s core,” Soviet Phys. Dokl. 149, 8- 10 (1963).
Braginsky, S. I., “Magnetohydrodynamics of Earth’s core,” Geomagn. Aeron. 4, 698-712 (1964).
Braginsky, S. I., “On a realistic geodynamo model,” J. Ceomagn. Geoelectr. 49, 1035- 1048 (1997).
Braginsky, S . 1. and Meytlis, V. P., “Local turbulence in the Earth’s core,” Geophys. Astrophys. Fluid Dynam. 55, 71 -87 (1990).
Braginsky, S. I. and Roberts, P. H., “Equations governing convection in Earth’s core and the geodynamo,” Rep?., IGPP, UCLA (1994).
Braginsky, S. I. and Roberts, P. H., “Equations governing convection in Earth’s core and the geodynamo,” Geophys. Asrrophys. Fluid D y m . 79, 1-97 (1995).
Braginsky, S. I. and Roberts, P. H., “On the theory of convection in the Earth’s core,” In: Advances in Non-linear Dynamos. (A. Femz-Mas and M. Niinez Eds.), London: Gordon and Breach, to appear (2000).
Buffett, B. A., Huppert, H. E., Lister, J. R. and Woods, A. W., “On the thermal evolution of the Earth’s core,” J . Geophys. Res. 101, 7989-8006 (1996).
Chabot, N. L. and Drake, M. J., “Potassium solubility in metal: The effects of composition at 15kbar and 1900” on partitioning between iron alloys and silicate melts,” Earth Planet. Sci. L e t f s . 172, 323-335 (1999).
Christensen, U., Olson, P. and Glatzmaier, G. A., “Numerical modelling of the geodynamo: systematic parameter study,” Geophys. J. Inr. 138, 393-409 (1999).
Connerney, J. E. P., “Magnetic fields of the outer planets,” J. Geophys. Res. 98, 18, 659- 18,679 (1993).
Diewonski, A. M. and Anderson, D. L., “Preliminary reference Earth model,” Phys. Earth PIanet. Infer.25, 297-356 (1981).
Glatzmaier, G .A., “Numerical simulations of stellar convective dynamos. I. The model and the method,” J. Compuf.Phys. 55,461 -484 (1984).
Glatzmaier, G. A., Coe, R. S., Hongre, L.and Roberts, P. H., “The role of the Earth’s
mantle in controlling the frequency of geomagnetic reversals,” Nature 401,885 -890 (1999). Glatzmaier, G. A. and Roberts, P. H., “A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle,” Phys. Earth Planet. Inter. 91, 63-75 (1995a). Glatzmaier, G. A. and Roberts, P. H., “A three-dimensional self-consistent computer simulation of the geomagnetic field,” Nature 377, 203-209 (1995b). Glatzmaier, G . A. and Roberts, P. H., “An anelastic evolutionary geodynamo simulation driven by compositional and thermal convection,” Physica D 97, 81-94 (1996a). This paper is referred to as ‘GRI’.

74

P. H. ROBERTS AND G. A. GLATZMAIER

Downloaded by [New York University] at 16:58 07 January 2015

Glatzmaier, G. A. and Roberts, P. H., “Rotation and magnetism of Earth’s inner core,” Science 274, 1887- 1891 (1996b).This paper is referred to as ‘GR2’.
Glatzmaier, G. A. and Roberts, P. H., “On the magnetic sounding of planetary interiors,” Phys. Earth Planet. Inter. 98, 207 -220 (1996~).
Glatzmaier, G. A. and Roberts, P. H., “Simulating the geodynamo,” Conremp.Phys. 38, 269-288 (1997). This paper is referred to as ‘GR3’.
Glatzmaier, G . A. and Roberts, P. H., “Dynamo theory then and now,” Znt. J. Engng. Sci. 36, 1325- 1338 (1998).
Grote, E., Busse, F. H. and Tilgner, A., “Regular and chaotic spherical dynamos,” Phys. Earth Planet. Inter. 117, 259-272 (2000).
Gubbins, D., Masters, T. G. and Jacobs, J. A., “Thermal evolution of the Earth’s core,” Geophys. J . R. Astr. SOC.59, 57-99 (1979).
Hide, R., “How to locate the electrically conducting fluid core of a planet from external magnetic observation,” Nature 271, 640-641 (1978).
Hollerbach, R. and Jones, C. A., “Influence of the Earth’s outer core on geomagnetic fluctuations and reversals,” Nature 365,541 - 543 (1993).
Jacobs, J. A,, “The Earth’s inner core,” Nature 172, 297 (1953). Jault, D., “Sur I’inhibition de la kgedration du champ magnktique dans certains
modtles de dynamo planktaire en pn5senced’une graine solide,” Comp. Rend. Acad. Sci. Paris 323, sir. IIa, 451 -458 (1996). Jones, C. A. and Roberts, P. H., “Convectiondriven dynamos in a rotating plane layer,” J. Fluid Mech. 404,31 1-343 (1999). Kono, M. and Tanaka, H., “Intensity of the geomagnetic field in geological time: a statistical study”. In: The Earth‘s Central Part: Its Structure and Dynamics. (T. Yukutake, Ed.), pp. 75-94, Tokyo: Terrapub (1996).
Kuang, W.and Bloxham, J., “An Earth-like numerical dynamo model,” Nature Js9,
371 -374 (1997).
Kutner, C. and Christensen, U., “Effects of driving mechanisms in geodynamo models,” Geophys. Res. Lett. 27, 29-32 (2000).
Labrosse, S.,Poirier, J.-P. and Le Mouel, J.-L., “On the cooling of Earth’s core,” Phys.
Earth Planet. Inter. 99, 1-17 (1997). Lambeck, K., The Earth’s Variable Rotation. Cambridge: University Press (1980).
Memll, R. T., McElhinny, M. W.and McFadden, P. L., The Magnetic Fieldof the Earth.
Paleomagnetism. the Core and the Deep Mantle. New Y o r k Academic Press (1996). Momson, G. and Fearn, D. R., “The influence of Rayleigh number, aximuthal
wavenumber and inner core radius on 2.5 D hydromagnetic dynamos,” Phys. Earth Planet. Inter. 117, 237-258 (2000). Munk, W. H. and Macdonald, G. F. J., Rotation of the Earth A Geophysical Discussion. Cambridge: University Press (1960). Olson, P., Christensen, U. and Glatzmaier, G. A., “Numerical modeling of the geodynamo: mechanisms of field generation and equilibration,” J. Geophys. Res. 104, 10,383- 10,404 (1999). Roberts, P. H. and Glatrmaier, G. A., “A test of the frozen flux approximation using geodynamo simulations,” Philos. Trans. R . SOC.London Ser. A 358, 1109-I121
(2000). Sakuraba, A. and Kono, M., “Effect of the inner core on the numerical solution of the
magnetohydrodynamic dynamo,” Phys. Earth Planer. Inter. 111, 105- 121 (1999). Sarson,G. R. and Jones, C. A., “A convection driven geodynamo reversal model,” Phys.
Earth Planer. Inter. 111, 3-20 (1999). Stacey, F. D., Physics of the Earth. Brisbane: Brook6eld Press (1992). Stevenson, D. J., Spohn, T. and Schubert, G., “Magnetism and thermal evolution of the
terrestrial planets,” I c a w 54, 466-489 (1983). Taylor, J. B., “The magneto-hydrodynamics of a rotating fluid and the Earth’s dynamo
problem.” Proc. R . SOC.Land. A 274, 274-283 (1963). Verhoogen, J., “Heat balance in the Earth’s core,” Geophys. J. 4, 276-281 (1961).

Downloaded by [New York University] at 16:58 07 January 2015

THE GEODYNAMO

75

APPENDIX A: CONSTRUCTION OF REFERENCE MODELS

In what follows we aim to develop three mutually consistent reference states, in the sense that the S state could evolve into the G state and the G state into the L state. There are considerable uncertainties in some of the parameters we require. Having assumed values for these parameters in the G model, we use values for the S and L models that are consistent with them. This has meant that we have retained more decimal places than is geophysically justifyable, but only because we wish to ensure consistency between the models.
Our reference state is an adaptation of that used in our GR simulations. It was assumed in GR1 that density in the FOC is (see also p. 80 of BR)

$ = c7(1 - cy),

(Al)

4 4 where and are constants chosen to agree with PREM densities at

CMB and ICB:

4 = 12484kgm-’, C ~ C ; I $ . ,=~ 2581 kgm-’.

(A21

We assume the mass of the inner core is M s ~ c= 9.84 x kg, which gives the acceleration, g, due to gravity at the ICB correctly (gICB= 4.40m sP2).This suffices to determine g throughout the FOC
as

where x =r/rCMB and
4 4 4 = 0.0250 m s-~, = 12.14 m s-~, = 0.1240. (A4)
The resulting gCMBdiffers from that of PREM by only 0.05%.
It is of course impossible with as simple a law as (Al) to obtain detailed agreement with the p(r) of PREM, but the errors are not great. For example, at r =2400 km,3 is 11259kg m - according to (Al), which differs from that of PREM by less than 0.3%. The mass
M F for~the fluid core is 1.837 x lo2*kg, which is less than that of

76

P. H. ROBERTS AND G . A. GLATZMAIER

PREM, but by under 0.2%. The velocity of sound, US, departs from that of PREM more significantly. Using hydrostatic balance in the form

Downloaded by [New York University] at 16:58 07 January 2015

we have

from which the sound velocities at CMB and ICB are 8437m s- and 10140ms-’, which exceed those of PREM by 5% and -2% respectively. We consider that (A6) provides an acceptable compromise, and that the constancy of r-‘dp/dr = - 2 4 4 implied by (Al) does not grossly violate the facts.
As in (Al), we model density in the SIC by

and since r-ld p/dr is almost the same on each side of the ICB, we take
cf4= 44.where cf is chosen to give Msrc correctly. We then have

.f = 13080kgmP3, C~C$&,.,~ = 2581 kgm-3.

(A8)

The central density is 13080kgm-3, which differs from the PREM

value by less than 0.07%. The density of the SIC at the ICB is

ZcBzcB 12762kgmP3which differs from that of PREM by less than 0.02%.

The density jump Ap = - m 6 gm

agrees with PREM.

The axial moments of inertia of the SIC and FOC are is, accord-

ing to (A7) and (Al), Zs~c= 5.856 x kg m2 and ZFW=

9.044 x kgm2.

In modifying this model for the past and future Earth, we make a

number of simplifying assumptions. We ignore changes of material

properties with time, and in particular we suppose that r-ldp/dr is
4 independent not only of r as above but also of 1. The constants and
cf are however necessarily time-dependent. We now write:

THE GEODYNAMO

77

’. where C =2581kg m - The implied masses of FOC and ICB are

Downloaded by [New York University] at 16:58 07 January 2015

where X =TICB/ICMB. Since the total mass of the core is constant, we therefore have

+ (1 - X3)cf X31$ = K ,

(A121

where, from today’s value of X =0.351 =X, (say), the constant K e I2510kgm-3.
Consider next the light constituent. BR introduced the “rejection factor”

which measures how much of the light constituent is rejected from the fluid as it freezes onto the ICB. We shall assume that this is constant and adopt the value 0.41, as suggested by BR. BR also introduce another significant parameter, bHL,which measures, peL,the effective density of light constituent when dissolved in the iron:
~ H =L -P-Fe 1,
PeL
where me is the density of iron. We assume that ~ H Lis constant:
6HL =0.7. This is much larger if Oxygen is the alloying element than it is if Silicon or Sulfur are: 0.95 as compared with 0.66 or 0.68. We ignore the Oxygen possibility. Also basic to the BR theory are

78

P. H. ROBERTS A N D G. A. GLATZMAIER

where a is the coefficient of volume expansion and cp is the specific heat at constant pressure; p is chemical potential. BR argued that, to a good approximation,

where the second expression follows from (D50)of BR. Here
A$ [ = f a y (say)] is that part of the discontinuity, AT, in density at
the ICB that can be attributed to the discontinuity, At, in composition; the remainder (1 -f)Ap is the result of contraction on freezing and (like E ) creates a general contraction of the core, which we ignore, SUppOSing throughout that rCMB = 3.48 X 1O6m.
The second of (A17) gives

Equating this to the value implied by (A9), we have

Downloaded by [New York University] at 16:58 07 January 2015

T Consider next conservation of the light constituent. Since convec-
tion thoroughly mixes the FOC, is uniform in space, though a

function of t. It is therefore clear that

M : =~J M~ F ~ .

(A211

The SIC is, on the present model, unmixed so that it is layered, the
p composition at radius r being the value of at the time that layer was p laid down, so that is a monotonically increasing function of r. We
have

THE GEODYNAMO

79

where the overdot denotes time derivatives. Since the total mass of light constituent is constant, we now obtain from (A13), (A21) and
(A22)

Downloaded by [New York University] at 16:58 07 January 2015

The equation may be integrated as

where ti, and M F X ,a~re the values of -EF ,-(and MFOCfor PDE.
We now equate the two expressions (AlO) and (A21) for M F E :

A(tF)-'""

=

K( 1

-

X3)

-

3
-

C

( 1

-

A')

5

- x3(1 - x ~ ) ( K- CX2)cFGHLrFs

f(1 +E F h ) -fFhLrFS(l - x3)>

(A25)

where A x 101.78kgmP3is chosen so that f F = 0.15 for PDE. The
values of for the ancient and future Earth are 0.145 and 0.200. Solutions of (A25) for other values of X are given in Table Al. It is clear that the SICSof the S and G models occupy such a small part of
tF the core that, if we took to be the same for each, the error made
would be less than those inherent in the uncertainties of other parameters. This is not so obviously true for the L model.

TABLE A1 Values of EF for different A

x

€F

x

EF

0.00

0.1448

0.50

0.1608

0.05

0.1449

0.55

0.1669

0.10

0.1450

0.60

0.1749

0.15

0.1452

0.65

0.1856

0.20

0.1458

0.70

0.2000

0.25

0.1467

0.75

0.2202

0.30

0.1480

0.80

0.2500

0.35

0.1500

0.85

0.2979

0.40

0.1526

0.80

0.3876

0.45

0.1561

0.95

0.6293

80

P. H. ROBERTS AND G. A. GLATZMAIER

It is now easy to evaluate quantities such as glCBandPICBw, here we
adopt for PCMBthe PREM value of 135.75GPa. Following BR, we
take

Downloaded by [New York University] at 16:58 07 January 2015

where the subscript m refers to liquidus values; y is the Griineisen
constant, which we take to be 1.27at the ICB, and KTis the isothermal
compressibility, which we take to be 1300GPa, as for PDE. By Taylor expansion, we therefore have approximately

+ Tm= To TT,P - 4375g.

(A271

where T = 1.44 x IO-l'Pa-'. The dependence of T, on tFis very
uncertain, the 700°K in (A26) being little more than an educated guess but, provided it is not grossly in error, the dependence of T,,, on p is much more significant,and we retain that alone, choosing the constant Toin (A27) to be 2795°K so that 7 1ag~ree~s with the value of 5300°K assumed for PDE by BR and GRI. Thus, we replace (A27) by the simpler

T , = 2795"K/(1 - T P ) ,

(A281

and so obtain TICBfrom pICBW. e estimate T ( r )elsewhere in the FOC
by integrating

We adopt for E / q an average value of 1.62 x lo-* kg/J, as used by BR. Equation (A29) is a consequence of hydrostatic equilibrium; see BR (3.7).
BR show that to a good approximation

which can be used both in determining the average rate of advance of the ICB and the rate of advance as a function of 0 and C#J on that surface. Here rZs and r y are constants and, taking values for the

Downloaded by [New York University] at 16:58 07 January 2015

THE GEODYNAMO

81

physical quantities given in Appendix E of BR, we have 2.148 x IO5"K
rx = rICB EICB '
rx = 1.931 x 1012m-20Ks-*
~ I C B~ I C BTICB

(A31a) (A31b)

To simplify the expression for r q slightly, we have used a repre-
7 . sentative value of 0.145 for Conservation of the light constituent
implies that

-
('=: a rICB
rICB

a = 5.1522 ~~CGLEL
M FOC

When combined with (A30), this gives [see Eq. (6.39) of BR]

To relate the rate of advance of the ICB to the heat withdrawn from the core by the mantle, BR show [see their (7.32)] that, in the absence of radioactive sources

Q = eE+es+e",

(A341

where Q5 represents energy released by gravitational differentiation, Qs that arising from the decreasing entropy of the core through cooling, and QNthat owing to the latent heat source:

(A35)

'. where h is the latent heat, which we take to be lo6J kg-
By (A15) and hydrostatic balance, @/dr = &g = -&dU/dr, where
17is gravitational potential. Since 55 is constant in the reference state

82

P. H. ROBERTS A N D G. A. GLATZMAIER

by (A17), this can be simply integrated to give

Downloaded by [New York University] at 16:58 07 January 2015

The evaluation of Qs requires a model of heat transport in the SIC. The following, though approximate, should not introduce significant errors, even in the L model. We suppose (see Section 1) that the SIC does not convect but, since its thermal conductivity is not negligible,
its entropy will not be frozen in; S will obey the heat conduction
equation, which is approximately

@F

as - at =

v

.

(

P

V

P

) +

&,

where E & ~is the Joule heating per unit volume and K is thermal
gIC, conductivity; see BR (5.19). Given T l c ~and (A39) can be solved
to give ? and hence the heat flux QsIc from the SIC into the ICB. We
can then calculate part of (A36):

gIc where is the net Joule heat loss in the SIC, which is partially offset

a&, by the greater heat flux QsIcarising from that internal heating. The

right-hand side of (AM) is unlikely to be less than

the adiabatic

heat flux to the ICB from an isentropic SIC. Following GR we assume

further, for simplicity, that that flux is continuous with the adiabatic

heat flux from the ICB into the FOC. In other words, the advance of

the ICB is principally due to convection in the FOC, as is also seen in

the first of the two boundary conditions at the ICB employed here and

in GRl -3:

THE GEODYNAMO

83

where X is the turbulent diffusivity of both heat and composition and
AS = ~ / T I CiBs the entropy release on freezing; see (6.29b) of BR. We
therefore replace (A40) by

Q&= &c.

('443)

The remaining part of Qsis due to the FOC and by (A36) is

Downloaded by [New York University] at 16:58 07 January 2015

Finally

By (A34), (A38), (A43), (A44) and (A45), we now have

z., 6 , The relative sizes of qs and qN indicate the relative importance of

each source in driving the convection. Since they and

can be

computed from the model (see Tab. A2), (A46) provides the means of

assessing the rate of growth, i . 1 ~ o~ f. the ICB from the heat flux from

the core to the mantle; see Table 1.

TABLE A2 Prouerties of Dast. Dresent and future models

Small

Generic

Large

Unit

0.147 0.620 0.0156 1.933 9.121 9930
12490
13090 600 174 0.0602

0.150 0.620 0.984 1.837 9.044 9903
12166
12762 597 189 0.0615

0.179 0.620 7.369 1.198 7.298 9710
11027
11665 638 235 0.0733

84

P.H.ROBERTS A N D G. A. GLATZMAIER

TABLE A2 (Continued)

Small

Generic

Large

1.118 357.3 5757 4183
0.00515 0.629 982 0.0146 0.00220 0.0171 0.00209 0.00468

4.401 328.2 5300 4Ooo
0.275 0.04oO 67.8 0.979 0.0504 0.352 0.1 12 0.292

8.287 238.8 4260 3627
1.33 0.0106 22.4
13.0 0.445
1.754 0.506 2.119

~
Unit
m/52
GPa
"K "K
Tw
kg"K/J
-
0.01 -
1030~
1030 J 1030 J

Downloaded by [New York University] at 16:58 07 January 2015

ZiB, It may be noted that the adiabatic heat flux,

from core to

mantle is 5.4TW for the G model according to Table I, but 5.2TW

according to GRl -3. This is because a polytropic model was adopted

for T in GRl, and this gives a slightly different value from the direct

integration of (A29) used here.