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.
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.
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.
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.
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 byeq. (1), the mass of the core Mc is given by Labrosse et al. 2001). Note that we incorporate the effect of the density jump when considering compositional convection (see below). eq. (5) assumes that the ratio αc/Cp is constant.
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)
Radioactive heating (QR, ER)
Eq. I–(26), is more complicated and depends on the relative sizes of D and L . For D>Lwe haveAbramowitz & Stegun 1972)
Latent heat (QL, EL)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
Eq. II–(39) givesGubbins et al. (1979) and is given by
Entropy of heat of solution (EH)eqs 11–15). Note that this expression assumes that core contraction is negligible (see below).
Adiabatic contribution (Q k , E k)
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
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:(31) and is given by
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)Solomatov 1995). 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: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 bySection 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 byeq. (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
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).
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.
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.
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.
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.
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.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, Δ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 , 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 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.
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 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).
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.
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. 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 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.
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 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: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.
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.
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.
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.