Summary

We model the thermal evolution of the core and mantle using a parametrized convection scheme, and calculate the entropy available to drive the geodynamo as a function of time. The cooling of the core is controlled by the rate at which the mantle can remove heat. Rapid core cooling favours the operation of a geodynamo but creates an inner core that is too large; slower cooling reduces the inner core size but makes a geodynamo less likely to operate. Introducing potassium into the core retards inner core growth and provides an additional source of entropy. For our nominal model parameters, a core containing ≈ 400 ppm potassium satisfies the criteria of present-day inner core size, surface heat flux, mantle temperature and cooling rate, and positive core entropy production.

We have identified three possibilities that may allow the criteria to be satisfied without potassium in the core. (1) The core thermal conductivity is less than half the generally accepted value of 50 W m−1 K−1. (2) The core solidus and adiabat are significantly colder and shallower than results from shock experiments and ab initio simulations indicate. (3) The core heat flux has varied by no more than a factor of 2 over Earth history.

All models we examined with the correct present-day inner core radius have an inner core age of <1.5 Gyr; prior to this time the geodynamo was sustained by cooling and radioactive heat production within a completely liquid core.

Introduction

The evolution of the Earth's core is one of the most fundamental topics in geophysics, relating as it does to the processes of planetary accretion and differentiation (e.g. Stevenson 1990); the history of mantle and core temperatures (e.g. Stacey & Loper 1984; Breuer & Spohn 1993); and the generation of the Earth's magnetic field. The latter topic has been investigated both analytically (e.g. Backus 1975; Loper 1978; Gubbins et al. 1979; Stevenson et al. 1983; Braginsky & Roberts 1995; Buffett et al. 1996) and through sophisticated numerical simulations (e.g. Glatzmaier & Roberts 1995; Kuang & Bloxham 1997; Olson et al. 1999; Kutzner & Christensen 2002).

Progress has recently been made in three areas relevant to core and geodynamo evolution. First, there is now a much better understanding of how to parametrize (Solomatov & Moresi 2000) or numerically model (Tackley 2000; Montague & Kellogg 2000) planetary heat transfer when plate tectonics, compositional contrasts, or strongly temperature-dependent fluids are involved. Secondly, ab initio simulations of core properties (Alfe et al. 2002a) are now in reasonable agreement with experimental determinations (Brown & McQueen 1986; Yoo et al. 1993), although some disagreements remain (Boehler 1993). Finally, there is now experimental evidence (Gessmann & Wood 2002; Lee & Jeanloz 2002; Murthy et al. 2003) to support older geochemical arguments (e.g. Goettel 1974) that radioactive elements may be present within the core, acting as an extra power source for the geodynamo.

In this paper, the thermal history of the core and mantle are modelled simultaneously using a parametrized convection scheme, and the entropy available to drive the geodynamo as a function of time calculated. A suite of different models with different core and mantle parameters are explored. Successful models must: (1) successfully reproduce the inferred present-day core and mantle temperature and viscosity structure; (2) successfully reproduce the present-day heat flux; and (3) generate enough entropy within the core to allow a geodynamo to function over the last 3 billion years (McElhinny & Senanayake 1980). These criteria are readily satisfied if the core contains a few hundred ppm potassium; successful models without core potassium can be constructed, but require more extreme parameter choices. In agreement with other models (e.g. Labrosse et al. 2001; Buffett et al. 1996), we find that the inner core is unlikely to exceed ≈1.5 Gyr in age; at earlier times, radioactive heat production and cooling of an entirely liquid core generates sufficient entropy to maintain the geodynamo.

Section 2 discusses previous contributions to the field, and compares the different approaches to that adopted in this paper. Section 3 lays out the theoretical framework, and Section 4 discusses the parameters adopted in this study. Section 5 outlines the results obtained for the nominal parameter set. The effects of the uncertainties in these parameters and comparisons with previous work are discussed in Section 6, and Section 7 summarizes the results.

Previous Work

Although the magnetohydrodynamic equations governing geodynamo behaviour are hard to solve numerically (e.g. Glatzmaier 2002), a necessary (though not sufficient) condition for geodynamo operation may be obtained by simply considering the energy or entropy production in the core. Braginsky (1964) and Gubbins et al. (1979), among others, give the basic equations for the energy and entropy changes in a cooling core. Stevenson et al. (1983) considered only energy, not entropy, but coupled the core thermal evolution to that of the mantle, and investigated the consequences for various terrestrial planets. Schubert et al. (1988) did the same for Mercury. Mollett (1984) used a similar coupled core–mantle approach, and also calculated the present-day entropy production. Although some more recent mantle thermal evolution models have included a core (e.g. Yuen et al. 1995; Honda & Iwase 1996; Tackley & Xie 2002), they do not generally include effects, such as entropy generation or latent heat release, relevant to the geodynamo. The approach of this paper is similar to that of Mollett (1984), although the details and the results are significantly different.

While the core cooling rate is ultimately determined by the rate at which heat can be extracted by the mantle, recent core calculations have tended to neglect the mantle thermal history and assume a core cooling rate by fiat. Buffett et al. (1996) obtained analytical expressions for the contributions to the core energy budget as a function of core growth; Roberts et al. (2003) took a similar approach and stressed the role of radioactivity in entropy production. Labrosse et al. (1997, 2001) concluded that the inner core was unlikely to be older than about 1.7 Gyr unless the core contained significant radioactive elements. Roberts & Glatzmaier (2001) performed a numerical investigation of likely terrestrial geodynamo behaviour for different inner core sizes, but did not explicitly examine the history of the core. Buffett (2002) calculated the core cooling rate based on a parametrized approach similar to the one adopted here, and also concluded that the core thermal history required radioactive elements within the core; he did not, however, calculate the entropy production with time.

Although previous approaches have often been analytical, the state of the Earth's inner core is well known, allowing a numerical approach to be used. Gubbins (2003a,b, hereafter G1 and G2, respectively) used seismological (e.g. Dziewonski & Anderson 1981) and mineralogical (e.g. Alfe et al. 2002a) data to peform numerical integrations giving the present-day rate of energy and entropy change. They did not, however, consider the time evolution of the core. This paper uses the same equations as G1 and G2, but adopts various approximations which allow analytical expressions to be obtained, in a similar manner to Buffett et al. (1996), Labrosse et al. (1997, 2001) and Roberts et al. (2003). The advantages of this approach are threefold: first, it allows the effects of uncertainties in the various parameters to be readily quantified; secondly, it allows the time-evolution problem to be investigated; and thirdly, the same approach can be applied to other planets (e.g. Schubert et al. 1988; Stevenson et al. 1983; Nimmo 2002) for which the relevant parameters are less well known.

A major focus of this work will be to explore the role of potassium in the core. Core potassium would provide an extra source of entropy (e.g. Mollett 1984) and a mechanism for slowing inner core growth (e.g. Labrosse et al. 2001), and its role will be examined in some detail below. Recently, three experimental groups (Murthy et al. 2003; Gessmann & Wood 2002; Lee & Jeanloz 2002) have independently argued that significant [O(100 ppm)] amounts of potassium may have partitioned into the core during its formation, contrary to previous expectations (e.g. Chabot & Drake 1999). The results are preliminary, and there is considerable uncertainty in extrapolations to the pressures and temperatures of equilibration, themselves uncertain. It is therefore of interest to see whether theoretical models of the kind presented here also support the presence of potassium in the core. Although, owing to uncertainties in other parameters, we cannot conclude that potassium is required, we will conclude that its presence in the core is highly probable.

Theory

The thermal evolution of a planet may be calculated if the sources of energy and the rates at which energy is transferred are both known. This section examines these two factors separately. In the first part, the sources of energy and entropy are derived, following the formulation in papers G1 and G2. In the second part, the rates of energy transfer are derived, and the expressions for planetary thermal evolution given. The final part gives the method for calculating the state of the inner and outer core.

Sources of entropy and energy

The calculations of G1 and G2 depend on a particular model of the density and temperature structure of the core (PREM: Dziewonski & Anderson 1981). It is therefore of interest to extend their calculations to a more general case. We will adopt the analytical density and temperature expressions obtained by Labrosse et al. (2001) to calculate entropy and energy. It will be shown below that the results obtained are very similar to those obtained by G1 and G2. An advantage of the analytical approach is that the effects of uncertainties in parameters are more easily investigated.

Preliminaries

The density within the Earth's core increases by about 35 per cent from the core–mantle boundary (CMB) to the centre of the planet. Following Labrosse et al. (2001), we assume that the variation of density ρ with radial distance r from the centre of the Earth is given by 

1
formula
where ρcen is the density at the centre of the Earth and L is a lengthscale given by 
2
formula
Here K0 and ρ0 are the compressibility and density at zero pressure, respectively, and G is the universal gravitational constant. Roberts et al. (2003) take a similar approach to the one here, but assume a more complicated form for the density variation. From eq. (1), the mass of the core Mc is given by 
3
formula
where R is the core radius. Similarly, the acceleration due to gravity g is given by 
4
formula
The error introduced by neglecting the density jump Δρ across the inner core boundary (ICB) is of order R4i/L4≪ 1, where Ri is the inner core radius (Labrosse et al. 2001). Note that we incorporate the effect of the density jump when considering compositional convection (see below).

The adiabatic temperature T within the core is given by 

5
formula
where Tcen is the temperature at the centre of the Earth and D is a lengthscale given by 
6
formula
Here Cp is the specific heat capacity and αc the thermal expansivity. Note that eq. (5) assumes that the ratio αc/Cp is constant.

Finally, the pressure is given by 

7
formula
where pc is the pressure at the core–mantle boundary (CMB).

Papers G1 and G2 assume an outer core that is well-mixed, close to hydrostatic and adiabatic outside thin boundary layers. The solidification of the inner core releases light material into the outer core, which may drive compositional convection; the cooling of the outer core drives thermal convection. Under these assumptions, the equations for changes in entropy E and energy Q during core cooling were obtained in G1 and G2. With the variation of ρ, g and T within the core described by eqs (1), (4) and (5), respectively, analytical solutions to these equations are derived below. Hereafter, equation numbers from paper G1 and G2 are denoted by I-(xx) and II-(xx), respectively.

Specific heat (Qs, Es)

Eq. I–(31), graphic , is simply given by (Labrosse et al. 2001)  

8
formula
where  
9
formula

The contributions Qs and Es can thus be derived: 

10
formula
where Mc is the mass of the core, and Tc the core temperature at the CMB.

Radioactive heating (QR, ER)

Eq. I–(26), graphic is more complicated and depends on the relative sizes of D and L . For D>Lwe have  

11
formula
where  
12
formula

For L > D we have 

13
formula
where 
14
formula
and Sn(B, R) is proportional to graphic and can be obtained using (Abramowitz & Stegun 1972) 
15
formula
These equations allow the contributions QR and ER to be determined: 
16
formula
where H is the heat production per unit mass within the core.

Latent heat (QL, EL)

The contributions QL and EL are easily obtained from I–(33) to I–(37): 

17
formula
where Lh is the latent heat of fusion, ρi is the density at the inner core boundary (ICB), which may be obtained using eq. (1), and Ti is the temperature at the ICB.

Gravitational contribution (Qg, Eg)

From (4), the potential ψ(r) relative to zero potential at the CMB is given by 

18
formula

Eq. II–(39) gives  

19
formula
where Moc is the mass of the outer core, ‘oc’ denotes the outer core, β c is a compositional expansion coefficient, Cr dTc / dt = dRi / dt , and  
20
formula
where χ is the weight concentration of element being expelled from the inner core during freezing. The compositional expansion coefficient β c is defined in II–(5) and Gubbins et al. (1979) and is given by  
21
formula
where Δρ is the density change across the ICB due to the presence of a light element in the outer core.

The integral graphic is given by 

22
formula
where 
23
formula

The mass of the outer core Moc can be obtained from (3) by changing the limits of integration. Eqs (19)–(23) allow Qg to be obtained; Eg is simply Qg/Tc.

Entropy of heat of solution (EH)

The quantity QH= 0. The contribution EH is given by II–(40): 

24
formula
where RH is the heat of reaction. The integral part of this equation can be derived using the same approach as that for radioactive heating (eqs 11–15). Note that this expression assumes that core contraction is negligible (see below).

Adiabatic contribution (Q k , E k)

The entropy change due to thermal diffusion is given by I–(67): 

25
formula
where k is the core thermal conductivity.

From (5) we have ∇T/T=−2r/D2; using (25) it can therefore be shown that 

26
formula
Qk at the CMB is simply given by 
27
formula

Other contributions

From the results of paper G2, the contribution of molecular diffusion and pressure changes are small compared to other terms and are ignored here. The rate of change of the core temperature is governed by 

28
formula
where u is the core contraction velocity. Paper G1 showed that core contraction is negligible at the present day, in which case DTc/Dt=dTc/dt.

Thermal evolution calculations

In carrying out thermal evolution calculations it is important to distinguish between potential and real temperatures. Hereafter the former will be denoted T′ and the latter T, with a subscript if applicable. The mantle potential temperature is simply the temperature obtained if the real temperature is projected along an adiabat to the surface: 

29
formula
where z is the depth and g, αm and Cpm are the acceleration due to gravity, thermal expansivity and specific heat capacity, all assumed constant through the mantle.

The thermal evolution of the core may be described by an equation similar to eq. (4a) of Schubert et al. (1988):  

30
formula
where QL, Qg, QR and Qs are defined above. graphic, and similarly for other quantities. The quantity QC is the heat extracted from the core by the mantle and is given by 
31
formula
where Fb is the heat flux across the core–mantle boundary.

In a similar fashion, the thermal evolution of the mantle is given by 

32
formula
where QM is the heat extracted from the mantle, Hm is the internal heat generation within the mantle, Mm and Cpm are the mantle mass and specific heat capacity, respectively, and Th is the temperature half-way through the mantle. The half-way point is used to represent the mean mantle temperature, but the results are insensitive to the depth chosen. The equation for QM is analogous to (31) and is given by 
33
formula
where Ft is the heat flux across the lithosphere and Rp is the planetary radius.

To determine the thickness of the top and bottom thermal boundary layers of the mantle, we use a local Rayleigh number approach (e.g. Howard 1964; Solomatov 1995). The top thermal boundary layer thickness δt is then given by (e.g. Nimmo & Stevenson 2000) 

34
formula
where Rac is the critical Rayleigh number, κt, ρm and αm the mantle thermal diffusivity, density and thermal expansivity, respectively, g is the surface acceleration and Ts the surface temperature. The viscosity at the mantle potential temperature ηt(Tm) is given by 
35
formula
where η0 is the viscosity at a reference temperature T0 and ζ is a quantity related to the activation energy (Solomatov 1995).

The heat flux Ft is then given by 

36
formula
where kt is the thermal conductivity at the top of the mantle and the potential temperature Tm is assumed to be a good approximation of the mantle temperature immediately beneath the top thermal boundary layer. Tm, Tm and Th may all be related through eq. (29). It will be shown below that this approach gives a reasonable description of plate tectonic heat fluxes at the present day.

Based on the experiments of Manga et al. (2001), we will assume that the heat flux across the bottom boundary is controlled by the formation of local instabilities (thermals) and is not greatly affected by descending slabs. In this case, expressions analogous to eqs (34)–(36) exist for the bottom thermal boundary layer δb and heat flux Fb:  

37
formula
 
38
formula
 
39
formula
Here κb and kb are the thermal diffusivity and conductivity at the base of the mantle, Ta is the mean of Tc and Tm (e.g. Manga et al. 2001), T1 is a reference temperature appropriate for the deep mantle and f is a constant which allows the deep mantle to have a higher viscosity than the shallow mantle, as is observed to be the case on the Earth.

The advantage of eqs (35) and (39) is that by choosing appropriate values for f, η0, T0 and T1 a realistic present-day mantle viscosity structure can be obtained (see Section 4 below). A similar approach is employed by Buffett (2002) and Yukutake (2000). Parametrized schemes such as the one adopted here accurately reproduce mean mantle temperatures calculated by numerical thermal evolution simulations for simple convection cases (e.g. Butler & Peltier 2000). However, it is not yet clear whether these parametrizations are appropriate to the case of plate tectonics, or how chemical layering (such as D′) will affect the results. We discuss this issue at some length in Section 6.2.3.

Given a set of initial conditions, the sources of energy can be calculated using Section 3.1, the rates of transfer of energy by eqs (31), (33)–(39), and the resulting temperature changes by eqs (32) and (30). Eqs (32) and (30) are integrated forwards in time using a constant timestep of 4 Myr; reducing the timestep to 1 Myr changes the results by one part in 105. The corresponding entropy production is calculated using eqs (10)–(27), assuming that all the material parameters are constant through time. Although in reality the density structure may change with time, for example due to the cooling of the core, such effects are likely to be small: for instance, a 1000 K change in core temperature will lead to a change in density of about 100 kg m−3, small compared with the density change across the core as a whole.

Temperature and solidification of core

The method followed here is essentially identical to that in Schubert et al. (1988). The melting curve of the core material as a function of pressure p is given by 

40
formula
where θ is a constant determined by the concentration of non-iron contaminants in the core, and their effect on the melting temperature. The constants Tm0, Tm1 and Tm2 are determined by specifying Tm0, Tma at one pressure and the gradient dTma/dp at another (see Section 4 below). The pressure at any depth is obtained using eq. (7). The constant θ is obtained by comparing the results for pure iron with those of iron alloyed with light elements such as sulphur, silicon and oxygen (see Section 4).

Following Schubert et al. (1988), the adiabatic temperature is given as a function of pressure by 

41
formula
where Ta1 and Ta2 are determined by fitting eq. (5) and relating depth to pressure via eq. (7). Eqs (40) and (41) may thus be used to determine the radius of the inner core for any Tc (see Schubert et al. 1988). The pressure at the centre of the Earth, pm, is used to check whether the core is entirely liquid.

Some of the light elements within the core are excluded from the inner core as it freezes. Schubert et al. (1988) assumed that this exclusion is complete. Here a more general case is implemented by specifying the partition coefficient Δ of the light element, where Δ is the ratio of the solid to liquid concentrations. For an initial concentration of χ0, the concentratio the outer core at subsequent times χ is given by 

42
formula
where ξ=Ri/R and the density variation within the core is here ignored. Because in general ξ3≪ 1, the results are insensitive to variations in Δ.

Parameters

Table 1 gives the parameters used in Section 3.1 to determine the energy and entropy production in the core. The values of K0 , ρ 0 and ρ cen are adopted so that the outer core density profile (eq. 1) agrees with the observations of PREM, and the value of ρ 0 agrees with the data of Anderson & Ahrens (1994).

Table 1

Values of quantities assumed for core calculations (Section 3.1).

Table 1

Values of quantities assumed for core calculations (Section 3.1).

The melting temperature for pure iron at ICB pressures (330 GPa) is 6300 ± 300 K, including a free-energy correction, and the gradient at 350 GPa is 9 K GPa−1 (Alfe et al. 2002a). Assuming that for h.c.p. iron Tm0= 1695 K, these data allow Tm1 and Tm2 to be determined using eqs (7) and (40). A temperature uncertainty of ±300 K results in a range for Tm1 of 9.9 to 12.0 × 10−12 Pa−1 and for Tm2 of −6.5 to −9.5 × 10−24 Pa−2. The ranges due to a gradient uncertainty of ±1 K GPa−1 are similar.

The density and temperature structure of the core is affected by the presence of one or more light elements. According to Alfe et al. (2002b), the density of the inner core can be explained the presence of about 8 molar per cent of S and/or Si, which reduces the temperature at the ICB by 700 ± 100 K. This temperature reduction implies θ= 0.11 ± 0.015 (eq. 40). Although θ is dependent on the concentration of contaminants, neither S nor Si are preferentially expelled during freezing (Alfe et al. 2002b), so θ is assumed to stay constant.

In addition to S and Si, the outer core contains oxygen, expelled as the inner core freezes, which further reduces the density of the outer core. It is the expulsion of oxygen from the inner core that provides the buoyancy to drive compositional convection. The amount of oxygen in the outer core may be inferred from the density drop across the ICB. Alfe et al. (2002b) used the density contrast from Shearer & Masters (1990) of 580 kg m−3 and inferred about 2.5 wt per cent O in the outer core. However, more recent normal-mode studies resulted in a higher density contrast of 820 ± 180 kg m−3 (Masters & Gubbins 2003). The higher value results in a range of O concentrations of 2.5–5.7 wt per cent. Here we assume an initial concentration χ0 of 4.2 wt per cent, and explore the effects of varying this amount below. The compositional expansion coefficient βc (eq. 21) obtained by Alfe et al. (2002b) differs by 10 per cent from that obtained by Roberts et al. (2003).

Based on the shock experiments of Matassov (1977), Stacey & Anderson (2001) concluded that the thermal conductivity of the core at the CMB was 46 W m−1 K−1, increasing to 63 W m−1 K−1 at the ICB. More recent shock experiments (Bi et al. 2002) suggest an electrical resistivity 50 per cent higher than that of Matassov (1977). Thus, the core thermal conductivity could be closer to 30 W m−1 K−1 at the CMB.

The thermal expansivity within the core varies by a factor of 2 (G1). However, adopting fixed values of αc and Cpc of 1.35 × 10−5 K−1 and 840 J kg−1 K−1, respectively, results in a similar adiabatic profile to those in G1. Varying the temperature at the ICB by ± 300 K gives a range in Ta1 and Ta2 of 3.1 to 3.9 × 10−12 Pa−1 and −1.2 to −2.4 × 10−24 Pa−2, respectively. The value of RH for the light element in the outer core is that quoted for oxygen in Paper G2; other core parameters are generally those used in papers G1 and G2.

Fig. 1 shows the resulting core temperature profile for the nominal parameters. The dashed line shows the solidus, including the effect of the light elements, using eq. (40) . The solid line is the adiabat from eq. (5) , and the dotted line that from eq. (41) . The two are identical within the likely error. The adiabatic gradient from eq. (5) is 7.5 K GPa−1 (= 0.38 K km−1) at the ICB, in agreement with recent calculations (paper G1). The relative slopes of the adiabat and the solidus determine the rate of inner core growth.

Figure 1

Present-day core thermal profile, using the nominal parameters of Tables 1 and 2 . The dashed line is the solidus calculated from eq. (40) , taking into account light elements. The solid line is the adiabat calculated from eq. (5) , and used in entropy and energy calculations. The dotted line is the adiabat calculated from eq. (41) , used in locating the inner core. ‘ICB’ is the inner core boundary, and ‘CMB’ is the core–mantle boundary. The temperature at the ICB is 5607 K; at the CMB the solidus temperature is 3562 K and the adiabat is 4161 K.

Figure 1

Present-day core thermal profile, using the nominal parameters of Tables 1 and 2 . The dashed line is the solidus calculated from eq. (40) , taking into account light elements. The solid line is the adiabat calculated from eq. (5) , and used in entropy and energy calculations. The dotted line is the adiabat calculated from eq. (41) , used in locating the inner core. ‘ICB’ is the inner core boundary, and ‘CMB’ is the core–mantle boundary. The temperature at the ICB is 5607 K; at the CMB the solidus temperature is 3562 K and the adiabat is 4161 K.

Table 2 gives the parameters used in Sections 3.2 and 3.3 . The projection of the mantle temperature of 1880 K at 670 km implies a CMB temperature of 2700 ± 200 K (Boehler 2000). Since the mantle potential temperature is about 1600 K (McKenzie & Bickle 1988 ; Stein & Stein 1992), the mean mantle thermal expansivity is 1.9–2.5 × 10−5 K−1 . We adopt here an intermediate value of 2.2 × 10−5 K−1 and discuss the effects of uncertainties below. The initial conditions are that Tc = Tm = 4800 K , which ensures that the core begins completely molten and that the mantle is close to its solidus (Stevenson 1990). It will be shown below that 1000 K variations in the initial conditions have negligible effects on the results.

Table 2

Nominal values of quantities assumed for thermal evolution calculations (Sections 3.2, 3.3).

Table 2

Nominal values of quantities assumed for thermal evolution calculations (Sections 3.2, 3.3).

Typical activation energies for mantle materials are 250–350 kJ mol−1 (Karato & Wu 1993; Yamazaki et al. 2000); for a mantle temperature of 1600–2500 K these energies imply values of ζ in the range 0.005–0.016 (Solomatov 1995), so a value of 0.01 is adopted here. The values of κt and κb are based on Hofmeister's (1999) calculations of the increase in conductivity with depth within the mantle; Manga & Jeanloz (1997) find similar or slightly higher values. Rac is generally assumed to be 600 (Nimmo & Stevenson 2000). The parameters Ta1, Ta2, Tm1 and Tm2 are obtained as described above. The heat generation rate within the mantle Hm is obtained from the radiogenic abundances of Sun & McDonough (1989).

The present-day viscosity of the mantle ranges from about 3–10 × 1020 Pa s in the upper mantle to 3–10 × 1021 Pa s in the lower mantle (e.g. Mitrovica & Forte 1997; Peltier & Jiang 1996). The behaviour of the bottom thermal boundary layer is determined by its mean temperature (Manga et al. 2001), which from the arguments above is about 3400 K. We adopt this value for T1 and 1573 K for T0. The nominal viscosity parameters adopted are shown in Table 2, but both the upper and lower mantle viscosities are varied, with effects which will be discussed below.

Results

The first part of this section compares the results of the analytical expressions above for present-day entropy production and power with the numerical results of papers G1 and G2. Having established that the analytical results do give a reasonable agreement, they are then used to explore various thermal evolution scenarios for the Earth, and to investigate their implications for geodynamo generation.

Comparisons with G1 and G2

Table 3 compares the results obtained in papers G1 and G2 by numerical integration with those obtained using the analytical methods outlined above. The pressure term (contributing <0.2 per cent of the total entropy budget) has been neglected. With the exception of EH (the smallest term), all terms agree to within 20 per cent or better, and the total heat and entropy production terms to within 10 per cent or better. Given the likely uncertainties in other parameters, we conclude that the analytical results reproduce the numerical ones sufficiently well, and we are thus justified in using the former to calculate entropy production and power over the core's lifetime.

Table 3

Comparison between Gubbins values and values obtained here.

Table 3

Comparison between Gubbins values and values obtained here.

Thermal evolution—nominal case

In order for a model to be judged successful, it must meet criteria regarding (1) the entropy generation, (2) the present-day inner core radius, and (3) the present-day mantle temperature, viscosity and heat flux.

The entropy available to drive the geodynamo, ΔE is given by 

43
formula
where Es, EL, EH and Eg depend on dTc/dt, ER depends on H, and Ek depends on the adiabat at the CMB. Although the effects of the geodynamo (i.e. the magnetic field) can be measured at the surface, the excess entropy production rate required to power it is unclear, for two reasons (Roberts et al., 2003). First, the measured magnetic field at the surface does not include the toroidal component, which may greatly exceed the poloidal component. Secondly, small lengthscale fluctuations in magnetic field and velocity may dominate Ohmic heating within the core, but are not observable at the surface. Based on numerical simulations, the excess entropy required is probably ∼ 100 MW K−1, but could lie anywhere within the range 0.1–1000 MW K−1 (G1, Roberts et al., 2003; Christensen et al. 1999; Kuang & Bloxham 1997). For the purposes of this work, any positive ΔE is assumed sufficient to drive the geodynamo, an assumption which is discussed further below. Three entropy criteria will be taken: first, ΔE at the present day must be positive; secondly, the mean of the entropy production over the last 3.1 Gyr graphic must also be positive; and finally, the minimum value of entropy production, ΔEmin, must be positive over the same period. The latter two criteria are a crude way of reproducing the approximately continuous geomagnetic record over that time interval (McElhinny & Senanayake 1980; Kroner & Layer 1992).

The other criteria are as follows. As discussed above, the current viscosities of the upper and lower mantle are likely to lie within the ranges 3–10 × 1020 Pa s and 3–10 × 1021 Pa s, respectively, and the present-day mantle potential temperature is 1330°C (McKenzie & Bickle 1988; Stein & Stein 1992). The present-day heat flux is approximately 42 TW (Sclater , 1980) and the dimensionless inner core radius at the present day is 0.35 (Dziewonski & Anderson 1981).

Fig. 2 shows the result of one thermal evolution calculation incorporating 400 ppm K in the core and using the nominal values of Tables 1 and 2 . It results in a positive value of Δ E and graphic , and a realistic present-day core and mantle structure. The results are tabulated in Table 4.

Figure 2

Thermal evolution and resulting entropy generation as a function of time. (a) Evolution of potential and real temperatures for the mantle (Tm, Tm) and core (Tc, Tc) with time. Ti is the temperature of the inner core boundary; prior to solidification, it represents the temperature at the centre of the planet. All parameters adopt the values given in Tables 1 and 2 , except the core heat generation rate, which is that due to 400 ppm potassium. The starting temperature of both mantle and core was 4800 K. (b) Evolution of heat fluxes (in TW) from the core (QC) and mantle (Qm) with time. The heat generated by radioactive decay in the mantle is given by HmMm (see eq. 32). (c) Solid and long-dashed lines give entropy contributions Es + EL + Eg + EH + ER and Ek , respectively, as a function of time. The available entropy (Δ E) is given by the difference between the two curves. The short-dashed line gives the dimensionless inner core radius ξ with time.

Figure 2

Thermal evolution and resulting entropy generation as a function of time. (a) Evolution of potential and real temperatures for the mantle (Tm, Tm) and core (Tc, Tc) with time. Ti is the temperature of the inner core boundary; prior to solidification, it represents the temperature at the centre of the planet. All parameters adopt the values given in Tables 1 and 2 , except the core heat generation rate, which is that due to 400 ppm potassium. The starting temperature of both mantle and core was 4800 K. (b) Evolution of heat fluxes (in TW) from the core (QC) and mantle (Qm) with time. The heat generated by radioactive decay in the mantle is given by HmMm (see eq. 32). (c) Solid and long-dashed lines give entropy contributions Es + EL + Eg + EH + ER and Ek , respectively, as a function of time. The available entropy (Δ E) is given by the difference between the two curves. The short-dashed line gives the dimensionless inner core radius ξ with time.

Table 4

Results of nominal thermal evolution calculation.

Table 4

Results of nominal thermal evolution calculation.

As expected, both the temperatures (Fig. 2a) and the heat fluxes (Fig. 2b) decline with time because of the decay of radioactive elements within the core and mantle. The present-day mantle potential temperature Tm is 1340° C, similar to the estimate of Stein & Stein (1992). Because of the way in which the model viscosity is chosen to depend on the temperature (eqs 32 and 36), a correct present-day temperature guarantees a realistic viscosity structure (see Table 4). The reduction in mantle temperature over the last 2.8 Gyr is 174 K, which is compatible with the petrological estimates of Abbott et al. (1994). The initial mantle temperature adjusts to its steady-state value within about 500 Ma because the mantle time constant is short at these high temperatures. The kink in the inner core temperature Ti at 3500 Ma is due to the initiation of core solidification; the present-day value of Ti is 5581 K, which agrees with the likely terrestrial value, as does the size of the inner core. The present-day surface heat flux (Fig. 2b) is about 42 TW, similar to the likely value. The thickness of the bottom boundary layer is 140 km, similar to estimates of D′ thickness. The total current heat flux from the core is 9 TW, which lies just within the 4.5–9 TW range recently proposed by Anderson (2002) and is similar to the heat fluxes deduced by Buffett (2002) and Yukutake (2000), but exceeds earlier estimates of around 2–3 TW (Sleep 1990; Stacey & Loper 1984). This total heat flux exceeds that carried by the adiabat, 6.2 TW, so that the outer core is likely to be convecting.

The entropy production (Fig. 2c) is high [O(1000 MW K−1)] early on, but declines steeply due to the rapid decay of potassium within the core and the reduction in the core cooling rate. These results strongly imply that an early geodynamo can be driven by purely thermal convection, without an inner core. Similar results were obtained by Mollett (1984) and Buffett et al. (1996). Just prior to core solidification, ΔE reaches its minimum value of 134 MW K−1. When the core begins to solidify, the available entropy increases again, because of the contributions from latent heat and gravity.

A model (not shown) identical to Fig. 2 but with no potassium in the core results in a core with a radius 80 per cent too large and a reduction in present-day core entropy production of 45 per cent. The role of potassium is therefore important in determining the history of the geodynamo and the inner core.

Varying the core cooling rate

Both the rate of inner core growth and the production of entropy are governed by the rate at which the core cools. A core that cools slowly will give rise to an inner core that is small at the present day; the slow cooling will mean a low rate of entropy production. The importance of potassium in the core is that it provides an extra source of entropy to power the dynamo, while at the same time retarding inner core growth. Fig. 2 shows that a model with 400 ppm potassium in the core can satisfy all available observational constraints. In this section, we show that it is hard to construct a model in which the observational constraints can be satisfied without having some potassium in the core. In the absence of potassium, it is difficult to satisfy the requirements of inner core size and entropy production simultaneously.

Fig. 3 shows the effect of different core cooling rates (controlled by the reference mantle viscosity) on the present-day available entropy production and inner core radius. Increasing the mantle reference viscosity reduces the core cooling rate, and thus the core size. The rate of entropy production depends on both the core size and the rate at which it is cooling. For mostly solid cores (ξ > 0.75) a reduction in cooling rate actually leads to an increase in entropy production, because of the concomitant increase in outer core volume. At smaller inner core sizes, however, a reduction in cooling rate decreases the rate of entropy production, because most entropy terms depend on dTc / dt . Figs 3(a) and (b) show the relationship between core size and entropy production when there is no potassium in the core. Obtaining an inner core of the correct size requires unrealistically large mantle viscosities (Fig. 3b). In these cases, the present-day entropy production is small (<100 MW K −1 ; Fig. 3a) and Δ Emin is negative (unfilled symbols), suggesting a cessation of dynamo activity in the past. Thus, for the nominal parameters, it is not possible to satisfy simultaneously the inner core radius and geodynamo requirements in the absence of potassium.

Figure 3

Variation in present-day entropy generation Δ E and inner core radius ξ with viscosity structure. The horizontal axis is the log of the reference viscosity η 0 . The lower mantle reference viscosity is given by f η 0 , where f is in the range 1–10 and different values of f are denoted by different lines. (a) Entropy generated with no potassium in the core. Solid symbols have positive Δ Emin ; open symbols have negative Δ Emin . (b) Present-day dimensionless inner core radius with no potassium. The horizontal dotted line denotes the observed value. (c) As (a) but with 400 ppm K in the core. (d) As (b) but with 400 ppm K in the core.

Figure 3

Variation in present-day entropy generation Δ E and inner core radius ξ with viscosity structure. The horizontal axis is the log of the reference viscosity η 0 . The lower mantle reference viscosity is given by f η 0 , where f is in the range 1–10 and different values of f are denoted by different lines. (a) Entropy generated with no potassium in the core. Solid symbols have positive Δ Emin ; open symbols have negative Δ Emin . (b) Present-day dimensionless inner core radius with no potassium. The horizontal dotted line denotes the observed value. (c) As (a) but with 400 ppm K in the core. (d) As (b) but with 400 ppm K in the core.

Figs 3(c) and (d) show the result for a core containing 400 ppm potassium. The radioactive decay in the core reduces the core cooling rate, which reduces ξ. The corresponding reduction in entropy production, however, is offset by the radioactive decay itself, which provides an additional entropy source (ER , eq. 16). Thus, for a given inner core radius, a potassium-rich core has more entropy available to drive the geodynamo. Figs 3(c) and (d) show that, for reasonable upper-mantle reference viscosities, the correct inner core radius can be produced at the same time as several hundred MW K −1 of entropy production with 400 ppm potassium in the core.

Fig. 4 shows the trade-off between inner core size and entropy production explicitly. The results of models using different viscosity structures are plotted and show that there is an almost linear relationship between graphic and ξ for ξ < 0.75 , irrespective of the particular viscosity structure adopted. High graphic requires large ξ, and vice versa; this is the same trade-off as identified in G2. For the zero potassium case, producing an inner core of the correct size requires a mean graphic of <100 MW K −1 and a negative rate of entropy production at earlier times. The effect of adding potassium to the core is to increase the mean entropy production for any inner core size. Adding 400 ppm potassium allows the inner core size to be satisfied with a mean entropy production of 200–300 MW K −1.

Figure 4

Variation in mean entropy production graphic with dimensionless inner core radius for various viscosity combinations and various potassium concentrations. Solid symbols have positive Δ Emin ; open symbols have negative Δ Emin . The vertical line indicates the present-day inner core radius. The effect of adding potassium is to increase the entropy production for the same inner core size.

Figure 4

Variation in mean entropy production graphic with dimensionless inner core radius for various viscosity combinations and various potassium concentrations. Solid symbols have positive Δ Emin ; open symbols have negative Δ Emin . The vertical line indicates the present-day inner core radius. The effect of adding potassium is to increase the entropy production for the same inner core size.

The present-day inner core radius and the age of the inner core are obviously linked. For the range of mantle parameters shown in Fig. 4, we found that it was impossible to obtain an inner core older than 1.5 Gyr whilst simultaneously producing a correct present-day inner core radius. Similar conclusions have been reached by previous authors (Yukutake 2000; G2; Labrosse et al. 2001; Buffett et al. 1996).

Discussion

The main results above may be summarized as follows. First, the analytical expressions of Section 3 give a reasonable agreement with the numerical results of papers G1 and G2. Secondly, our model can only match both the entropy production and the inner core radius simultaneously if the core contains ≈ 400 ppm potassium (Fig. 4). The conclusion that potassium is required to make the models work is important, because it provides an additional argument for the presence of potassium in the core. We therefore devote the first part of this section to investigating the uncertainties in this conclusion, and to possible scenarios that do not require potassium in the core. We then compare the results of this investigation with those of other authors, to point out similarities and differences. Finally, we make suggestions for future work.

Uncertainties

As with any analysis involving large numbers of parameters, it is important to try to quantify the sensitivity of the results to the uncertainties in those parameters. Table 5 lists the amount by which each parameter in turn has to be varied in order to satisfy the present-day inner core radius without requiring any potassium in the core. While this is not an exhaustive test (e.g. it does not vary all parameters simultaneously), it does allow us to identify which parameters are the most important (see Section 6.2).

Table 5

Uncertainty analysis

Table 5

Uncertainty analysis

The rate of inner core growth is controlled by the slope of the solidus and adiabat. One would therefore expect that variables that change the solidus (e.g. Tm1) or adiabat (e.g. αc) would have a significant effect on the results, and Table 5 shows that this is indeed the case. The value of Tm1 required to fit the observations with no potassium corresponds to a reduction in Tc of 400 K and lies just outside the estimated bounds on Tm1. Furthermore, the resulting value of ΔEmin is negative, which suggests that this model is not appropriate. We return to this issue below.

The value of αc (and the corresponding values of Ta1 and Ta2) required to generate a correct inner core size in the absence of potassium does not fall within the likely uncertainty range. Similarly, the likely uncertainties in αm, θ, and κb (see Table 2) are all smaller than the values required from Table 5. All these parameters also produce negative values of ΔEmin.

Another way of reducing the core growth rate is to reduce the heat flux out of the core. Table 5 shows that a reduction in ζ or an increase in η0 or f might be able to produce an inner core of the correct size without a requirement for potassium. However, the required values of f or η0 result in lower mantle viscosities in the range 2–6 × 1022 Pa s, which is probably too large. Furthermore, in both cases ΔEmin is negative, due to the trade-off between core growth rate and entropy production shown in Fig. 4. Reducing the value of ζ means that the core and mantle heat flux are less sensitive to variations in temperature, and provides a possible mechanism for matching the observations without requiring potassium in the core. We discuss this issue further below.

Reducing the core thermal conductivity or increasing the expansion coefficient βc or oxygen concentration χ increases the available entropy production, but has no effect on the rate of inner core growth. The effect is therefore to move points in Fig. 4 vertically upwards. Table 5 shows that a larger mantle reference viscosity can produce an inner core of the correct size, but that ΔEmin is negative. For this higher reference viscosity, increasing βc by 50 per cent still results in negative ΔEmin, and is a much larger change than the estimated error. Increasing χ0 to 5.7 wt per cent increases the present-day entropy production by 25 per cent but has no effect on ΔEmin. Reducing the conductivity to 20 W m−1 K−1 results in positive entropy production throughout (because Ek is reduced) and is discussed further below.

Other parameters are unlikely to affect the overall conclusions. For instance, the potassium-free model results are rather insensitive to the starting conditions. Increasing the initial temperature (Tm=Tc) by 1000 K results in changes of less than 8 per cent in graphic and 3 per cent in ξ, the dimensionless core radius. The lack of sensitivity is due to the short mantle time constant at high temperatures (see Fig. 2). Reducing Tm and Tc makes the problem of inner core growth worse and requires correspondingly higher core potassium values to fit the present-day observations. Similarly, changing Lh or Cp by 20 per cent results in a change in ξ of <5 per cent.

Alternative scenarios

Based on the results of Table 5, there are three possibilities that would allow the observations to be satisfied without requiring potassium in the core. We discuss these possibilities below, and summarize the results in Table 6.

Table 6

Alternative scenarios which do not require potassium in the core.

Table 6

Alternative scenarios which do not require potassium in the core.

Reduced thermal conductivity

Reducing the conductivity to 20 W m−1 K−1 results in an acceptable model, but at the expense of mantle viscosities, which are higher than the likely values (see Table 6). Furthermore, the required conductivity value is lower than current estimates of the likely range (see Section 4). We therefore consider this scenario unlikely.

Different core properties

It is possible that the adiabat and solidus of the core are better represented by the experimental results of Boehler (2000) than by the ab initio studies used above. The principal differences between the experiments and the numerical simulations are that the former give a solidus that has a lower temperature and a shallower slope, and the effect of the alloying element on melting is much reduced. The mean thermal expansivity implied by the results of Boehler (2000) is about half the nominal value used in this paper.

Using the parameters given in the footnote in Table 6, and keeping all other parameters at their nominal values, we find that a reasonable fit to the observations can be obtained without requiring potassium in the core. The principal difference from the nominal model is that more heat needs to be extracted from the core before solidification can begin; consequently, the mantle reference viscosity has to be significantly reduced. Although the positive contributions to entropy are reduced, the conductive contribution Ek is also reduced because of the lower slope of the adiabat, with the result that the net entropy production remains positive.

Different heat flux parametrization

An alternative is that our parametrization of the core heat flux, eq. (38), does not capture the physics of the CMB region. We accordingly examined the effects of a different parametrization, in which the core heat flux simply declines exponentially, irrespective of conditions in the mantle: 

44
formula
where t is time, τ is an adjustable parameter, and F0 was chosen so that the present-day core heat flux is 6.2 TW, equal to the adiabatic contribution Qk. Table 6 shows the results when τ= 6500 Myr, and demonstrates than an acceptable model can be obtained without requiring potassium. The principal difference from the nominal model is that the time-averaged heat flux out of the core is lower—the core heat flux at t= 0 is only a factor of 2 higher than that at the present day.

Whether this scenario is likely is as yet unclear. Our original parametrization (eq. 38) effectively assumes that descending subducted slabs have little effect on the heat transfer across the base of the mantle. Laboratory results (Schaeffer & Manga 2001) show that the frequency of rising plumes due to instabilities in the bottom boundary layer is a factor of 2–3 higher, and of similar temperature amplitude, to the rising plumes caused by descending cold plumes. Thus, in these experiments the principal cause of the heat flux across the bottom boundary layer is a local instability, similar to the model we adopt. Similarly, Labrosse (2002) suggests that the bottom heat flux scales as Ra1/4ΔT, where ΔT is the temperature drop across the bottom boundary layer and Ra is the Rayleigh number. While the exponents in this scaling are slightly different from those in eq. (38), the basic result (that the heat flux across the CMB will decrease markedly as the mantle cools) remains the same. Even if it is the descending slabs that dominate the heat flux across the CMB, scaling arguments suggest that their thickness will vary as Ra−1/3 (eq. 34). Requiring the CMB heat flux to vary by less than a factor of 2 thus places a stringent constraint on the variation of the mantle Rayleigh number over 4.5 Gyr. Given that the upper mantle potential temperature has apparently reduced by 137–187 K over the last 2.8 Gyr (Abbott et al. 1994), such a constraint seems unlikely to be satisfied.

Comparisons with other work

Table 4 shows that the entropy contributions are dominated by those from gravity (Eg) and latent heat release (EL), the former being about twice as large. Conversely, latent heat release (QL) and cooling (Qs) dominate the power contribution, while gravity (Qg) accounts for less. Thus, as expected, the gravitational contribution to geodynamo generation is much more efficient than the thermal contribution (G2; Gubbins et al. 1979 ; Buffett et al. 1996 ; Stevenson et al. 1983).

In the absence of potassium, both entropy production and inner core growth rate have a nearly linear dependence on the core cooling rate dTc/dt (G1; G2), so the ratio of the two quantities is almost independent of dTc/dt (Fig. 4). This independence explains why it is difficult to obtain a small inner core and an active geodynamo simultaneously, and has been noted by other authors (G2; Buffett 2002, Roberts et al., 2003). The approach of Roberts et al. (2003) is similar to that in this paper, although their choice of parameters alters the relative importance of the gravitational and latent heat terms, and they assume a core cooling rate rather than calculating one. Like us, they conclude that the presence of core potassium allows the constraints of a small (young) inner core and an ancient magnetic field to be satisfied simultaneously.

A long-standing problem in Earth history is that the planet appears to be losing heat at about twice the rate at which heat is being produced; that is, it has a Urey ratio of ≈ 0.5 (McKenzie & Richter 1981). Such a low ratio is not expected for a conventional single-layered mantle, and implies ancient mantle temperatures higher than those inferred from petrological observations (e.g. Abbott et al. 1994). The model shown in Fig. 2(a) is compatible with these observations, however, because the core heat contribution reduces the mantle cooling rate. The effect of core potassium on mantle cooling rates was previously noted by Breuer & Spohn (1993).

For core heat fluxes similar to that in Table 4, Labrosse et al. (2001) required 120 ppm potassium (with U and Th in chondritic ratios) to obtain an inner core age of 2.7 Gyr. Mollett (1984) obtained an inner core age of about 1.7 Gyr with a present-day core radiogenic heat flux of 2 TW, equivalent to 50 ppm potassium (with chondritic U and Th). Buffett (2002) obtained an inner core age of 1.8 Gyr with 1.4 TW of radiogenic heating. Using the Fig. 2 model but with 120 ppm core potassium and U and Th in chondritic proportions results in the correct inner core radius, an age of 1.4 Gyr, ΔE= 264 MW K−1 and 4.9 TW of heat production in the core. Although, as noted above, we can construct scenarios that do not require potassium in the core, there appears to be general agreement that 2–5 TW of heat production in the core makes it far easier to produce plausible core thermal histories. Strictly speaking, this heat production could be due to sources other than radioactivity. Tidal heating in the core is insignificant at the moment, though it may have been higher earlier on (Munk & MacDonald 1960); it is not clear what other energy sources might be present.

Future work

The greatest uncertainty in this work arises from the parametrization of the CMB heat flux adopted. As discussed above, it may be too simple; future models will need to investigate the effects on the CMB heat flux of chemical layering (Montague & Kellogg 2000; Jellinek & Manga 2002), mantle phase transitions (McKenzie & Richter 1981; Butler & Peltier 2000), and plate tectonics (Labrosse 2002).

A second source of uncertainty is in the correct parameters to use for the core, in particular its solidus behaviour and thermal conductivity. While there is agreement between shock measurements and ab initio simulations for the former, static experiments are still producing different results. Regarding the latter, the difference between the results of Matassov (1977) and Bi et al. (2002) is significant, and it will be important to investigate this further.

A final source of uncertainty is the entropy production rate required to drive a geodynamo. The 2 TW Ohmic heating estimate of Roberts et al. (2003) implies ΔE≈ 400 MW K−1, which would severely limit the acceptable models and make the presence of potassium in the core highly probable (see Fig. 4). A less extreme value of ≈100 MW K−1 makes it easier for models lacking potassium to drive a dynamo. It is obviously of great interest to place better theoretical bounds on the entropy production required to generate a geodynamo. It is also unclear what happens to the dynamo if the outer part of the core becomes stably stratified, as G2 suggest. Further numerical simulations of geodynamo activity (e.g. Kuang & Bloxham 1997; Glatzmaier 2002; Kutzner & Christensen 2002) may be helpful in resolving these questions.

Summary

A self-consistent model has been developed in which the entropy available to drive a geodynamo is calculated from the thermal evolution of the core and mantle. There is a trade-off between entropy production and present-day inner core size (Fig. 4). For our nominal model, sustaining a geodynamo and obtaining an inner core of the correct size is only possible with ≈ 400 ppm potassium in the core. Experimental results have recently shown that such concentrations are possible. There are, however, at least three alternatives to the requirement of potassium in the core. One is that the thermal conductivity of the core is less than half of the generally accepted value; this option also requires a mantle viscosity higher than the likely values. The second possibility is that the behaviour of the core solidus and adiabat are better described by the static experiments of Boehler (2000) than the ab initio simulations of Alfe et al. (2002a,b). A third possibility is that our description of the heat flux across the CMB is too simple. An acceptable model can be obtained using an alternative parametrization (eq. 44), but only if the CMB heat flux has varied by less than a factor of 2 over 4.5 Gyr. All our models require an inner core that is younger than about 1.5 Gyr BP; prior to that, the geodynamo was probably sustained by a mixture of radioactive decay and cooling in an entirely liquid core.

Acknowledgments

FN is supported by the Royal Society and NSF-EAR 0309218. We are grateful to Orson Anderson, Rama Murthy, Paul Roberts and Jerry Schubert for helpful discussions, and to two anonymous reviewers for their comments.

References

Abbott
D.L.
Burgess
L.
Longhi
J.
Smith
W.H.F.
,
1994
.
An empirical thermal history of the Earth's upper mantle
,
J. geophys. Res.
 ,
99
,
13835
13850
.
Abramowitz
A.
Stegun
I.A.
,
1972
.
Handbook of Mathematical Functions
 ,
National Bureau of Standards
,
Washington DC
.
Alfe
D.
Price
G.D.
Gillan
M.J.
,
2002
.
Iron under Earth's core conditions: liquid-state thermodynamics and high-pressure melting curve from ab initio calculations
,
Phys. Rev. B.
 ,
65
, art. 045123.
Alfe
A.
Gillan
M.J.
Price
G.D.
,
2002
.
Composition and temperature of the Earth's core constrained by combining ab initio calculations and seismic data
,
Earth planet. Sci. Lett.
 ,
195
,
91
98
.
Anderson
O.L.
,
2002
.
The power balance at the core–mantle boundary
,
Phys. Earth planet. Inter.
 ,
131
,
1
17
.
Anderson
W.W.
Ahrens
T.J.
,
1994
.
An equation of state for liquid-iron and implications for the Earth's core
,
J. geophys. Res.
 ,
99
,
4273
4284
.
Backus
G.E.
,
1975
.
Gross thermodynamics of heat engines in deep interior of Earth
,
Proc. Nat. Acad. Sci. USA
 ,
72
,
1555
1558
.
Bi
Y.
Tan
H.
Jin
F.
,
2002
.
Electrical conductivity of iron under shock compression up to 200 GPa
,
J. Phys. condensed Matter
 ,
14
,
10849
10854
.
Boehler
R.
,
1993
.
Temperature in the Earth's core from melting-point measurements of iron at high static pressures
,
Nature
 ,
363
,
534
536
.
Boehler
R.
,
2000
.
High pressure experiments and the phase diagram of lower mantle and core materials
,
Rev. Geophys.
 ,
38
,
221
245
.
Braginsky
S.I.
,
1964
.
Magnetodhydrodynamics of the Earth's core
,
Geomag. Aeron.
 ,
4
,
989
916
.
Braginsky
S.I.
Roberts
P.H.
,
1995
.
Equations governing convection in Earth's core and the geodynamo
,
Geophys. astrophys. Fluid Dyn.
 ,
79
,
1
97
.
Breuer
D.
Spohn
T.
,
1993
.
Cooling of the Earth, Urey ratios, and the problem of potassium in the core
,
Geophys. Res. Lett.
 ,
20
,
1655
1658
.
Brown
J.M.
McQueen
R.G.
,
1986
.
Phase transitions, Grueneisen parameter, and elasticity for shocked iron between 77 and 400 GPa
,
J. geophys. Res.
 ,
91
,
7485
7494
.
Buffett
B.A.
Huppert
H.E.
Lister
J.R.
Woods
A.W.
,
1996
.
On the thermal evolution of the Earth's core
,
J. geophys. Res.
 ,
101
,
7989
8006
.
Buffett
B.A.
,
2002
.
Estimates of heat flow in the deep mantle based on the power requirements for the geodynamo
,
Geophys. Res. Lett.
 ,
29
, DOI: .
Butler
S.L.
Peltier
W.R.
,
2000
.
On scaling relations in time-dependent mantle convection and the heat transfer constraint on layering
,
J. geophys. Res.
 ,
105
,
3175
3208
.
Chabot
N.L.
Drake
M.J.
,
1999
.
Potassium solubility in metal: the effects of composition at 15kbar and 1900 degrees C on partitioning between iron alloys and silicate melts
,
Earth planet. Sci. Lett.
 ,
172
,
323
335
.
Christensen
U.
Olson
P.
Glatzmaier
G.
,
1999
.
Numerical modelling of the geodynamo: a systematic parameter study
,
Geophys. J. Int.
 ,
138
,
393
409
.
Dziewonski
A.M.
Anderson
D.L.
,
1981
.
Preliminary reference Earth model
,
Phys. Earth. planet. Inter.
 ,
25
,
297
356
.
Gessmann
C.K.
Wood
B.J.
,
2002
.
Potassium in the Earth's core?
,
Earth planet. Sci. Lett.
 ,
200
,
63
78
.
Glatzmaier
G.A.
,
2002
.
Geodynamo simulations - how realistic are they?
,
Ann. Rev. Earth planet. Sci.
 ,
30
,
237
257
.
Glatzmaier
G.A.
Roberts
P.H.
,
1995
.
A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle
,
Phys. Earth planet Int.
 ,
91
,
63
75
.
Goettel
K.A.
,
1974
.
Potassium in the Earth's core: evidence and implications
, in
Physics and Chemistry of Minerals and Rocks
 , pp.
479
489
, ed.
Strews
R.G.
,
John Wiley Interscience
,
New York
.
Gubbins
D.
Masters
T.G.
Jacobs
J.A.
,
1979
.
Thermal evolution of the Earth's core
,
Geophys. J. R. astr. Soc.
 ,
59
,
57
99
.
Gubbins
D.
Alfe
D.
Masters
T.G.
Price
D.
Gillan
M.J.
,
2003
.
Can the Earth's dynamo run on heat alone?
,
Geophys. J. Int.
 ,
155
,
609
622
.
Gubbins
D.
Alfe
D.
Masters
T.G.
Price
D.
Gillan
M.J.
,
2003
.
Gross thermodynamics of 2-component core convection
,
Geophys. J. Int.
 , in press.
Hofmeister
A.M.
,
1999
.
Mantle values of thermal conductivity and the geotherm from phonon lifetimes
,
Science
 ,
283
,
1699
1706
.
Honda
S.
Iwase
Y.
,
1996
.
Comparison of the dynamic and parameterized models of mantle convection including core cooling
,
Earth planet. Sci. Lett.
 ,
139
,
133
145
.
Howard
L. N.
,
1964
.
Convection at high Rayleigh numbers
, in ed.
Görtler
H.
,
Proc. 11th Int. Conf. on Applied Mechanics
 , pp.
1109
1115
,
Springer-Verlag
, New York.
Jellinek
A.M.
Manga
M.
,
2002
.
The influence of a chemical boundary layer on the fixity, spacing and lifetime of mantle plumes
,
Nature
 ,
418
,
760
763
.
Karato
S.
Wu
P.
,
1993
.
Rheology of the upper mantle: A synthesis
,
Science
 ,
260
,
771
778
.
Kroner
A.
Layer
P.W.
,
1992
.
Crust formation and plate motion in the early Archean
,
Science
 ,
256
,
1405
1411
.
Kuang
W.
Bloxham
J.
,
1997
.
An Earth-like numerical dynamo model
,
Nature
 ,
389
,
371
374
.
Kutzner
C.
Christensen
U.R.
,
2002
.
From stable dipolar towards reversing numerical dynamos
,
Phys. Earth planet. Inter.
 ,
131
,
29
45
.
Labrosse
S.
,
2002
.
Hotspots, mantle plumes and core heat loss
,
Earth planet. Sci. Lett.
 ,
199
,
147
156
.
Labrosse
S.
Poirier
J.-P.
Le Mouel
J.-L.
,
1997
.
On cooling of the Earth's core
,
Phys. Earth. planet. Inter.
 ,
99
,
1
17
.
Labrosse
S.
Poirier
J.-P.
Le Mouel
J.-L.
,
2001
.
The age of the inner core
,
Earth planet. Sci. Lett.
 ,
190
,
111
123
.
Lee
K.K.
Jeanloz
R.
,
2002
.
Equation of state of K+Fe at high pressure: K in the Earth's core?
,
EOS Trans. Am. geophys. Un.
 ,
83
, Fall Meet. Suppl., Abstract MR62B–1065.
Loper
D.E.
,
1978
.
The gravitationally powered dynamo
,
Geophys. J. R. astr. Soc.
 ,
54
,
389
404
.
Manga
M.
Jeanloz
R.
,
1997
.
Thermal conductivity of corundum and periclase and implications for the lower mantle
,
J. geophys. Res.
 ,
102
,
2999
3008
.
Manga
M.
Weeraratne
D.
Morris
S.J.S
,
2001
.
Boundary-layer thickness and instabilities in Benard convection of a liquid with temperature-dependent viscosity
,
Phys. Fluids
 ,
13
,
802
805
.
Masters
T.G.
Gubbins
D.
,
2003
.
On the resolution of density within the Earth
,
Phys. Earth planet. Inter.
 ,
140
,
159
167
.
Matassov
G.
,
1977
.
The electrical conductivity of iron-silicon alloys at high pressures and the Earth's core
,
PhD thesis
 , Lawrence Livermore Laboratory, University of California, CA.
McKenzie
D.P.
Bickle
M.J.
,
1988
.
The volume and composition of melt generated by extension of the lithosphere
,
J. Petrol.
 ,
29
,
625
679
.
McKenzie
D.
Richter
F.
,
1981
.
Parameterized convection in a layered region and the thermal history of the Earth
,
J. geophys. Res.
 ,
86
,
11667
11680
.
McElhinny
M.W.
Senanayake
W.E.
,
1980
.
Paleomagnetic evidence for the existence of the geomagnetic field 3.5 Ga ago
,
J. geophys. Res.
 ,
85
,
3523
3528
.
Mitrovica
J.X.
Forte
A.M.
,
1997
.
Radial profile of mantle viscosity: results from the joint inversion of convection and postglacial rebound observables
,
J. geophys. Res.
 ,
102
,
2751
2769
.
Mollett
S.
,
1984
.
Thermal and magnetic constraints on the cooling of the Earth
,
Geophys. J. R. astr. Soc.
 ,
76
,
653
666
.
Montague
N.L.
Kellogg
L.H.
,
2000
.
Numerical models of a dense layer at the base of the mantle and implications for the geodynamics of D"
,
J. geophys. Res.
 ,
105
,
11101
11114
.
Munk
W.H.
MacDonald
G.J.F.
,
1960
.
The Rotation of the Earth: a Geophysical Discussion
 ,
Cambridge Univ. Press
,
Cambridge
, p.
323
.
Murthy
V.R.
Van Westrenen
W.
Fei
Y.
,
2003
.
Experimental evidence that potassium is a substantial radioactive heat source in planetary cores
,
Nature
 ,
423
,
163
165
.
Nimmo
F.
,
2002
.
Why does Venus lack a magnetic field?
,
Geology
 ,
30
,
987
990
.
Nimmo
F.
Stevenson
D.
,
2000
.
The influence of plate tectonics on the thermal evolution and magnetic field of Mars
,
J. geophys. Res.
 ,
105
,
11969
11980
.
Olson
P.
Christensen
U.
Glatzmaier
G.A.
,
1999
.
Numerical modeling of the geodynamo: mechanisms of field generation and equilibration
,
J. geophys. Res.
 ,
104
,
10383
10404
.
Peltier
W.R.
Jiang
X.
,
1996
.
Mantle viscosity from simultaneous inversion of multiple data sets pertaining to postglacial rebound
,
Geophys. Res. Lett.
 ,
23
,
503
506
.
Roberts
P.H.
Glatzmaier
G.A.
,
2001
.
The geodynamo, past, present and future
,
Geophys. astrophys. Fluid Dyn.
 ,
94
,
47
84
.
Roberts
P.H.
Jones
C.A.
Calderwood
A.R.
,
Energy fluxes and ohmic dissipation in the Earth's core
, in
Earth's Core and Lower Mantle
 , eds
Jones
C.A.
Soward
A.M.
Zhang
K.
, Book Series The Fluid Mechanics of Astrophysics and Geophysics (Series ed. A.M. Soward & Ghil,) Vol.
11
.
Taylor & Francis
(in press).
Schaeffer
N.
Manga
M.
,
2001
.
Interaction of rising and sinking mantle plumes
,
Geophys. Res. Lett.
 ,
28
,
455
458
.
Schubert
G.
Ross
M.N.
Stevenson
D.J.
Spohn
T.
,
1988
.
Mercury's thermal history and the generation of its magnetic field
, in
Mercury
 , pp.
429
460
, eds
Vilas
F.
Chapman
C.R.
Matthews
M.S.
,
Univ. Ariz. Press
,
Tucson, AZ
.
Sclater
J.G.
Jaupart
C.
Galson
D.
,
1980
.
The heat flow through oceanic and continental crust and the heat loss of the Earth
,
Rev. Geophys. Space Phys.
 ,
18
,
269
311
.
Shearer
P.M.
Masters
T.G.
,
1990
.
The density and shear velocity contrast at the inner core boundary
,
Geophys. J. Int.
 ,
102
,
491
498
.
Sleep
N.H.
,
1990
.
Hot spots and mantle plumes: some phenomenology
,
J. geophys. Res.
 ,
95
,
6715
6736
.
Solomatov
V. S.
,
1995
.
Scaling of temperature- and stress-dependent viscosity convection
,
Phys. Fluids
 ,
7
,
266
274
.
Solomatov
V.S.
Moresi
L.-N.
,
2000
.
Scaling of time-dependent stagnant lid convection: application to small-scale convection on Earth and other terrestrial planets
,
J. geophys. Res.
 ,
105
,
21795
21817
.
Stacey
F.D.
Anderson
O.L.
,
2001
.
Electrical and thermal conductivities of Fe-Ni-Si alloy under core conditions
,
Phys. Earth planet. Inter.
 ,
124
,
153
162
.
Stacey
F.D.
Loper
D.E.
,
1984
.
Thermal histories of the core and mantle
,
Phys. Earth planet. Inter.
 ,
36
,
99
115
.
Stein
C.A.
Stein
S.
,
1992
.
A model for the global variation in oceanic depth and heat-flow with lithospheric age
,
Nature
 ,
359
,
123
129
.
Stevenson
D.J.
,
1990
.
Fluid dynamics of core formation
, in
Origin of the Earth
 , pp.
231
250
, eds
Newsom
H.E.
Jones
J.E.
,
Oxford Univ Press
,
Oxford
.
Stevenson
D.J.
Spohn
T.
Schubert
G.
,
1983
.
Magnetism and thermal evolution of the terrestrial planets
,
Icarus
 ,
54
,
466
489
.
Sun
S. S.
McDonough
W.F.
,
1989
.
Chemical and isotopic systematics of oceanic basalts: implications for mantle composition and processes
,
Geol. Soc. Spec. Pub.
 ,
42
,
313
345
.
Tackley
P.J.
,
2000
.
Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations Part 1: Pseudo-plastic yielding
,
Geochem., Geophys., Geosys.
 ,
1
, Paper number 2000GC000036.
Tackley
P.J.
Xie
S.
,
2002
.
The thermo-chemical structure and evolution of Earth's mantle: constraints and numerical models
,
Phil. Trans. R. Soc. Lond.
  A.,
360
,
2593
2609
.
Yamazaki
D.
Kato
T.
Yurimoto
H.
Ohtani
E.
Toriumi
M.
,
2000
.
Silicon self-diffusion in MgSiO3 perovskite at 25 GPa
,
Phys. Earth planet. Inter.
 ,
119
,
299
309
.
Yoo
C.S.
Holmes
N.C.
Ross
M.
Webb
D.J.
Pike
C.
,
1993
.
Shock temperatures and melting of iron at Earth core conditions
,
Phys. Rev. Lett.
 ,
70
,
3931
3934
.
Yuen
D.A.
Balachandar
S.
Steinbach
V.C.
Honda
S.
Reuteler
D.M.
Smedsmo
J.J.
Lauer
G.S.
,
1995
.
Non-equilibrium effects of core cooling and time-dependent internal heating on mantle flush events
,
Nonlinear Processes Geophys.
 ,
2
,
206
221
.
Yukutake
T.
,
2000
.
The inner core and the surface heat flow as clues to estimating the initial temperature of the Earth's core
,
Phys. Earth planet. Inter.
 ,
121
,
103
137
.