## 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

_{cen}is the density at the centre of the Earth and

*L*is a lengthscale given by Here

*K*

_{0}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

*M*

_{c}is given by where

*R*is the core radius. Similarly, the acceleration due to gravity

*g*is given by The error introduced by neglecting the density jump Δρ across the inner core boundary (ICB) is of order

*R*

^{4}

_{i}/

*L*

^{4}≪ 1, where

*R*

_{i}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

*T*

_{cen}is the temperature at the centre of the Earth and

*D*is a lengthscale given by Here

*C*

_{p}is the specific heat capacity and α

_{c}the thermal expansivity. Note that eq. (5) assumes that the ratio α

_{c}/

*C*

_{p}is constant.

Finally, the pressure is given by

where*p*

_{c}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 (Q_{s}, E_{s})

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

The contributions *Q*_{s} and *E*_{s} can thus be derived:

*M*

_{c}is the mass of the core, and

*T*

_{c}the core temperature at the CMB.

#### Radioactive heating (Q_{R}, E_{R})

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

*S*(

_{n}*B*,

*R*) is proportional to and can be obtained using (Abramowitz & Stegun 1972) These equations allow the contributions

*Q*

_{R}and

*E*

_{R}to be determined: where

*H*is the heat production per unit mass within the core.

#### Latent heat (Q_{L}, E_{L})

The contributions *Q*_{L} and *E*_{L} are easily obtained from I–(33) to I–(37):

*L*

_{h}is the latent heat of fusion, ρ

_{i}is the density at the inner core boundary (ICB), which may be obtained using eq. (1), and

*T*

_{i}is the temperature at the ICB.

#### Gravitational contribution (Q_{g}, E_{g})

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

Eq. II–(39) gives

where*M*

_{oc}is the mass of the outer core, ‘oc’ denotes the outer core, β

_{c}is a compositional expansion coefficient,

*C*

_{r}dT_{c}/

*dt*=

*dR*

_{i}/

*dt*, and 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 where Δρ is the density change across the ICB due to the presence of a light element in the outer core. where

The mass of the outer core *M*_{oc} can be obtained from (3) by changing the limits of integration. Eqs (19)–(23) allow *Q*_{g} to be obtained; *E*_{g} is simply *Q*_{g}/*T*_{c}.

#### Entropy of heat of solution (E_{H})

The quantity *Q*_{H}= 0. The contribution *E*_{H} is given by II–(40):

*R*

_{H}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):

where*k*is the core thermal conductivity.

From (5) we have ∇*T*/*T*=−2*r*/*D*^{2}; using (25) it can therefore be shown that

*Q*

_{k}at the CMB is simply given by

#### 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

where*u*is the core contraction velocity. Paper G1 showed that core contraction is negligible at the present day, in which case

*DT*

_{c}/

*Dt*=

*dT*

_{c}/

*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:

*z*is the depth and

*g*, α

_{m}and

*C*

_{pm}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):

*Q*

_{L},

*Q*

_{g},

*Q*

_{R}and

*Q*

_{s}are defined above. , and similarly for other quantities. The quantity

*Q*

_{C}is the heat extracted from the core by the mantle and is given by where

*F*

_{b}is the heat flux across the core–mantle boundary.

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

where*Q*

_{M}is the heat extracted from the mantle,

*H*

_{m}is the internal heat generation within the mantle,

*M*

_{m}and

*C*

_{pm}are the mantle mass and specific heat capacity, respectively, and

*T*

_{h}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

*Q*

_{M}is analogous to (31) and is given by where

*F*

_{t}is the heat flux across the lithosphere and

*R*

_{p}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)

*Ra*

_{c}is the critical Rayleigh number, κ

_{t}, ρ

_{m}and α

_{m}the mantle thermal diffusivity, density and thermal expansivity, respectively,

*g*is the surface acceleration and

*T*

_{s}the surface temperature. The viscosity at the mantle potential temperature η

_{t}(

*T*′

_{m}) is given by where η

_{0}is the viscosity at a reference temperature

*T*

_{0}and ζ is a quantity related to the activation energy (Solomatov 1995).

The heat flux *F*_{t} is then given by

*k*

_{t}is the thermal conductivity at the top of the mantle and the potential temperature

*T*′

_{m}is assumed to be a good approximation of the mantle temperature immediately beneath the top thermal boundary layer.

*T*

_{m},

*T*′

_{m}and

*T*

_{h}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 *F*_{b}:

_{b}and

*k*

_{b}are the thermal diffusivity and conductivity at the base of the mantle,

*T*

_{a}is the mean of

*T*

_{c}and

*T*

_{m}(e.g. Manga

*et al.*2001),

*T*

_{1}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}, *T*_{0} and *T*_{1} 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 10^{5}. 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

*T*

_{m0},

*T*

_{m1}and

*T*

_{m2}are determined by specifying

*T*

_{m0},

*T*

_{ma}at one pressure and the gradient

*dT*

_{ma}/

*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

*T*

_{a1}and

*T*

_{a2}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

*T*

_{c}(see Schubert

*et al.*1988). The pressure at the centre of the Earth,

*p*

_{m}, 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

*R*

_{i}/

*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 *K*_{0} , ρ _{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).

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 *T*_{m0}= 1695 K, these data allow *T*_{m1} and *T*_{m2} to be determined using eqs (7) and (40). A temperature uncertainty of ±300 K results in a range for *T*_{m1} of 9.9 to 12.0 × 10^{−12} Pa^{−1} and for *T*_{m2} 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 *C*_{pc} 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 *T*_{a1} and *T*_{a2} of 3.1 to 3.9 × 10^{−12} Pa^{−1} and −1.2 to −2.4 × 10^{−24} Pa^{−2}, respectively. The value of *R*_{H} 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.

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 *T*_{c} = *T*_{m} = 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.

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. *Ra*_{c} is generally assumed to be 600 (Nimmo & Stevenson 2000). The parameters *T*_{a1}, *T*_{a2}, *T*_{m1} and *T*_{m2} are obtained as described above. The heat generation rate within the mantle *H*_{m} is obtained from the radiogenic abundances of Sun & McDonough (1989).

The present-day viscosity of the mantle ranges from about 3–10 × 10^{20} Pa s in the upper mantle to 3–10 × 10^{21} 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 *T*_{1} and 1573 K for *T*_{0}. 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 *E*_{H} (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.

### 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

*E*

_{s},

*E*

_{L},

*E*

_{H}and

*E*

_{g}depend on

*dT*

_{c}/

*dt*,

*E*

_{R}depends on

*H*, and

*E*

_{k}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 must also be positive; and finally, the minimum value of entropy production, Δ

*E*

_{min}, 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 × 10^{20} Pa s and 3–10 × 10^{21} 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 , and a realistic present-day core and mantle structure. The results are tabulated in Table 4.

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 *T*′_{m} 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 *T*_{i} at 3500 Ma is due to the initiation of core solidification; the present-day value of *T*_{i} 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 *dT*_{c} / *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 Δ *E*_{min} 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.

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 (*E*_{R} , 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 and ξ for ξ < 0.75 , irrespective of the particular viscosity structure adopted. High 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 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}.

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).

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. *T*_{m1}) 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 *T*_{m1} required to fit the observations with no potassium corresponds to a reduction in *T*_{c} of 400 K and lies just outside the estimated bounds on *T*_{m1}. Furthermore, the resulting value of Δ*E*_{min} 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 *T*_{a1} and *T*_{a2}) 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 Δ*E*_{min}.

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 × 10^{22} Pa s, which is probably too large. Furthermore, in both cases Δ*E*_{min} 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 Δ*E*_{min} is negative. For this higher reference viscosity, increasing β_{c} by 50 per cent still results in negative Δ*E*_{min}, 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 Δ*E*_{min}. Reducing the conductivity to 20 W m^{−1} K^{−1} results in positive entropy production throughout (because *E*_{k} 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 (*T*_{m}=*T*_{c}) by 1000 K results in changes of less than 8 per cent in 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 *T*_{m} and *T*_{c} makes the problem of inner core growth worse and requires correspondingly higher core potassium values to fit the present-day observations. Similarly, changing *L*_{h} or *C*_{p} 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.

#### 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 *E*_{k} 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:

where*t*is time, τ is an adjustable parameter, and

*F*

_{0}was chosen so that the present-day core heat flux is 6.2 TW, equal to the adiabatic contribution

*Q*

_{k}. 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 *Ra*^{1/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 (*E*_{g}) and latent heat release (*E*_{L}), the former being about twice as large. Conversely, latent heat release (*Q*_{L}) and cooling (*Q*_{s}) dominate the power contribution, while gravity (*Q*_{g}) 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 *dT*_{c}/*dt* (G1; G2), so the ratio of the two quantities is almost independent of *dT*_{c}/*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

*PhD thesis*

*The Fluid Mechanics of Astrophysics and Geophysics*(Series ed. A.M. Soward & Ghil,) Vol.