# Electrical resistivity and thermal conductivity of liquid Fe alloys at high *P* and *T*, and heat flux in Earth’s core

See allHide authors and affiliations

Edited by Peter L. Olson, Johns Hopkins University, Baltimore, MD, and approved January 11, 2012 (received for review July 21, 2011)

## Abstract

Earth’s magnetic field is sustained by magnetohydrodynamic convection within the metallic liquid core. In a thermally advecting core, the fraction of heat available to drive the geodynamo is reduced by heat conducted along the core geotherm, which depends sensitively on the thermal conductivity of liquid iron and its alloys with candidate light elements. The thermal conductivity for Earth’s core is very poorly constrained, with current estimates based on a set of scaling relations that were not previously tested at high pressures. We perform first-principles electronic structure computations to determine the thermal conductivity and electrical resistivity for Fe, Fe–Si, and Fe–O liquid alloys. Computed resistivity agrees very well with existing shock compression measurements and shows strong dependence on light element concentration and type. Thermal conductivity at pressure and temperature conditions characteristic of Earth’s core is higher than previous extrapolations. Conductive heat flux near the core–mantle boundary is comparable to estimates of the total heat flux from the core but decreases with depth, so that thermally driven flow would be constrained to greater depths in the absence of an inner core.

The generation of the Earth’s magnetic field is directly coupled to the thermal evolution of the liquid outer core, the cooling of which is modulated by heat flux from the core into the base of the mantle (1⇓⇓⇓⇓–6). In a simple thermally driven core, as, for example, in the absense of a crystallizing inner core and associated latent heat release and chemical buoyancy, heat transported by conduction is not available to drive the geodynamo. Knowledge of its relative contribution to thermal transport in the core is therefore critical to understanding the long term stability of Earth’s magnetic field, which has been present as far back as 3.45 billion years ago (7), when the core was likely too hot for a solid inner core to crystallize (1, 4).

Existing estimates of thermal conductivity (*k*_{el}) and electrical resistivity (*ρ*_{el}) of Earth’s outer core are based on extrapolations (8, 9) of resistivity measurements in shock-compressed Fe and Fe–Si alloys (10⇓–12) to core temperatures and pressures. These extrapolations assume the direct proportionality of electrical resistivity to temperature, its invariability along and across the Fe liquidus, and adherence to the Wiedemann–Franz law, which relates electrical resisitivity and thermal conductivity for metals through the Lorenz number (13) *λ*_{0} = *k*_{el}*ρ*_{el}/*T* = 2.44 × 10^{-8} WΩ/K^{2}. No data are available for Fe alloys of the other candidate light elements that have been proposed to account for the seismically observed density deficit of Earth’s core relative to pure Fe (14). Previous high-pressure studies of the electronic transport properties of Fe at high pressure (15⇓⇓–18) were limited to low temperatures. Clearly, there is a need for direct determination of electrical resistivity and thermal conductivity of ferro-metallic liquids at pressures and temperatures characteristic of Earth’s outer core.

We compute *k*_{el} and *ρ*_{el} for Fe, Fe_{7}Si, Fe_{3}Si, Fe_{7}O, and Fe_{3}O liquid alloys (6.7, 14.3 wt% Si; 3.9, 8.7 wt% O) from first principles, using density functional theory and the Mermin functional to determine finite temperature equilibrium charge density and electronic structure (19⇓–21). First-principles molecular dynamics (FPMD) simulations are performed in the canonical ensemble for temperatures of 2,000–8,000 K and volumes corresponding to pressures of 0–360 GPa. Electronic transport properties are subsequently computed for a series of uncorrelated snapshots from the FPMD simulations using the Kubo–Greenwood equation (22, 23) (see *Methods*), which expresses the electronic Onsager coefficients *L*_{ij} directly in terms of the expectation values of the electronic velocity operator (24).

Computed resistivities for pure Fe liquid (Fig. 1) closely agree with the shock compression measurements of Keeler (10) and the lowest pressure point of Bi et al. (11), at pressures where Hugoniot temperatures (25) are comparable to those in our simulations. Similarly, values for Fe_{3}Si liquid are in agreement with the shock measurements of Matassov (12) for the same composition. Lower pressure shock compression measurements are at progressively lower temperatures; the 18 GPa measument of Keeler (10) is at 320 K, and agrees well with room temperature static measurements in hcp Fe (17, 18), adding confidence to the shock compression resistivity data, and our results. In light of the large scatter in the Bi et al. (11) dataset, and the serious disagreement between the data of Bi et al. (11) and Keeler (10) above 120 GPa, we feel that comparison of our results with the measurements of Keeler (10) is more appropriate. Computed low pressure resistivity for pure Fe liquid is somewhat smaller than experimental values (26, 27), yet computed ambient pressure thermal conductivity does agree closely with experimental estimates for liquid Fe at temperatures above melting (1,810 K) (28).

In drawing comparisons to the shock compression resistivity measurements, we assume that the resistivity of liquid and solid metal phases are similar, which is known to be the case for Fe at low pressure (26, 27). In our comparison of thermal conductivities, we further assume the electronic contribution to heat transport to be much larger than that due to transport by phonons alone. In Fe liquid at ambient pressure, the latter is estimated to be about 3 W/m K (8), much smaller than the experimental total conductivity of 40.3 W/m K (28).

In contrast with previous assumptions for Earth’s core (8, 9), we find that the validity of the Wiedemann–Franz relation depends strongly on temperature and composition. At core temperatures, computed Lorenz numbers for Fe and Fe–Si liquids are in the range of 2.2–2.4 × 10^{-8} WΩ/K^{2}; alloying of liquid Fe with O results in a notable decrease in *λ*, with values as low as 1.8 × 10^{-8} WΩ/K^{2} for Fe_{3}O liquid at high pressure and temperture. This deviation from the Wiedemann–Franz relation suggests that electron scattering in Fe–O liquid alloys is strongly anelastic, resulting in a breakdown of the simple relaxation time picture of electronic migration.

We find electrical resistivity to vary linearly with temperature, consistent with the prediction of the Bloch–Grüneisen equation for systems where electrons are primarily scattered by phonons. However, *ρ*_{el} is not directly proportional to *T* (Fig. 2), as is assumed in the oft-used extrapolations of experimental measurements to core conditions (8, 9). Augmenting the Bloch–Grüneisen formalism (29), we construct model descriptions for *k*_{el} and *ρ*_{el} as a function of volume and temperature for the different liquid phases considered (Table 1; see *Methods*). Combining these models with a pressure–volume–temperature equation of state for conductive liquids (30), we derive *ρ*_{el} and *k*_{el} values along a set of candidate adiabatic thermal profiles for Earth’s core, derived from the range of Fe-melting temperatures (25, 31, 32).

Our *k*_{el} values for the outer core are notably higher than previous extrapolations [Fig. 3; 30–60 W/m K (8, 9)], reflecting the incorrect temperature dependence of resistivity assumed in earlier studies (Fig. 2). As a consequence, our predicted conductive heat flux at the top of the core is 14–20 TW, larger than the 5–15 TW of total heat flux across the core–mantle boundary inferred from intraplate volcanism (5). A lower heat flux at the top of the core would require a smaller *k*_{el}, which can be obtained through a larger concentration of light elements in this region. Indeed, an anomalous light element fraction at the top of the core has been proposed on dynamical grounds (33) and is supported by seismic observations (34). However, because this concentration of light elements likely results from inner core crystallization and the associated chemical buoyancy, the difficulty posed by a large conductive heat flux remains, especially in the absence of a crystallizing inner core.

Model conductive heat flux decreases rapidly with depth in the core, primarily due to the spherical geometry. Excluding the possible effects of radiogenic heating and latent heat released through inner core crystallization, conservation of energy requires total heat flux to remain constant throughout the core. The fraction of the total heat flux transported via convection would therefore increase with depth, suggesting that thermally driven flow would occur mainly within the deeper portion of the liquid core.

## Methods

### First-Principles Simulations.

The method used to compute electronic transport properties is similar to that of refs. 22 and 23. Born–Oppenheimer molecular dynamics (FPMD) is performed using the VASP code (35). To test for finite size effects, we consider system sizes of 64, 128, 144, and 192 atoms (see below). The exchange-correlation potential is represented in the generalized gradient approximation (GGA-PBE) (36), with valence electrons represented as planewaves with a cutoff of 300 eV in the projector augmented wave (PAW) formalism (37, 38). The Brillouin zone is sampled at the Γ point only. Simulations are performed in the *NVT* ensemble for volumes *V*/*V*_{X} = 1.0, 0.7, 0.65, and 0.6, where *V*_{X} = 7.121 cm^{3}/mol atom, and temperatures of 2,000, 3,000, 4,000, 6,000 and 8,000 K, and cover at least 20 ps of simulation time. The time dependent mean square displacement is used to check that systems are indeed in the liquid state.

From each FPMD-generated phase trajectory, we extract atomic configuration snapshots every 1,000 fs (i.e., 20 per *P*–*T* point) for which we compute the electrical resistivity (*ρ*_{el}) and thermal conductivity (*k*_{el}). Velocity autocorrelation functions (39) for our simulations all decay within 250 fs, indicating that a 1,000-fs time separation is sufficient for individual snapshots to be uncorrelated. In this way a representative sampling of the liquid structure at each *P*–*T* point is obtained.

Electronic transport properties *ρ*_{el} and *k*_{el} are computed using the Kubo–Greenwood equation, as implemented in the Abinit code (22, 40). The equation, which follows from the electronic current autocorrelation function via Kubo’s linear response formalism, is [1][2][3]In Eqs. **1**–**3** *ϵ*_{F} is the Fermi energy; *ψ*_{k}, *ϵ*_{k}, and *f*(*ϵ*_{k}) are the wave function, eigenvalue, and Fermi–Dirac occupation of eigenstate *k*, respectively; is the velocity operator; and *V*_{cell} is the simulation cell volume. For a given snapshot, the self-consistent electronic relaxation is performed for electronic temperature equal to the ionic temperature via the Mermin functional (21). *ψ*_{k} and *ϵ*_{k} are represented by the Kohn–Sham eigenfunctions and eigenvalues for each given snapshot, while is computed from the Hamiltonian gradient, .

Cells of 144 atoms for Fe and 128 atoms for Fe–Si and Fe–O alloys are used in production runs; test simulations with 192 atoms at *V*/*V*_{X} = 0.6 and 1.0 for *T* = 8,000 K yielded transport coefficients within 1% of the values determined with the smaller systems. To avoid core charge overlap in the linear response calculations at high degrees of compression, we constructed GGA-PAW atomic potentials with small cutoff radii (0.9 Å for Fe and Si, 0.53 Å for O) (36, 41). A planewave basis set cutoff of 400 eV is found to yield converged electronic transport coefficients. Production runs sample the Brillouin zone only at the Γ-point. This choice is appropriate in large simulation cells, where the Brillouin zone edge is effectively folded into its center. Test values computed using a 2 × 2 × 2 Monkhorst–Pack **k**-point grid (42) at *V*/*V*_{X} = 0.6 and 1.0 for *T* = 8,000 K vary by no more than 5% from the single **k**-point results.

It should be noted that the Kubo–Greenwood equation determines the electronic transport properties directly from the self-consistent electronic structure within the Born–Oppenheimer approximation. Ionic and electronic scattering are implicitly included in the computation, but electron–phonon coupling is not described. The approach is therefore well suited to characterizing high temperature electronic transport coefficients, especially in liquid metals.

*P*-*T* Model for *ρ*_{el} and *k*_{el}.

We represent the temperature and volume dependence of electrical resistivity via the Bloch–Grüneisen equation, with volume dependence included through a power law dependence of the prefactor (29) [4]where *a* and *b* are constants that include volume dependence of the vibrational frequencies. The electronic thermal conductivity is represented self-consistently through the Lorenz number using [5]by fitting the computed Lorenz numbers to [6]

For the Fe liquid alloys we consider, reference volume and temperature are chosen as *T*_{R} = 3,000 K and *V*_{R} = 7.0 cm^{3}/mol atom. Values for *ρ*_{0R}, *ρ*_{1R}, *λ*_{R}, *a*, *b*, *c*, and *d* for each of the simulated systems is given in Table 1.

## Acknowledgments

We thank Henri Samuel for helpful discussions. This work was made possible by support from the Deutsche Forschungsgemeinschaft under Contract KO3958/2-1 in the focus program Planetary Magnetism. Computing facilities were provided by the Leibniz Supercomputing Centre of the Bavarian Academy of Sciences and Humanities.

## Footnotes

- ↵
^{1}To whom corresopondence should be addressed. E-mail: nico.dekoker{at}uni-bayreuth.de.

Author contributions: N.d.K. and G.S.-N. designed research; N.d.K. and V.V. performed research; N.d.K. analyzed data; and N.d.K. and G.S.-N. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- Glatzmaier GA,
- Coe RS,
- Hongre L,
- Roberts PH

- ↵
- ↵
- ↵
- ↵
- ↵
- Tarduno JA,
- et al.

- ↵
- ↵
- ↵
- Caldirola P,
- Knoepel H

- Keeler RN

- ↵
- ↵
- Matassov G

- ↵
- Ashcroft NW,
- Mermin ND

- ↵
- ↵
- ↵
- Sha X,
- Cohen RE

- ↵
- ↵
- Balchan A,
- Drickamer H

- ↵
- ↵
- ↵
- Mermin ND

- ↵
- Recoules V,
- Crocombette J.-P

- ↵
- Holst B,
- French M,
- Redmer R

- ↵
- Ziman JM

- ↵
- ↵
- Zytveld JBV

- ↵
- ↵
- Touloukian YS,
- Powell RW,
- Ho CY,
- Klemens PG

- ↵
- ↵
- ↵
- Ma YZ,
- et al.

- ↵
- Alfè D

- ↵
- Buffett BA,
- Garnero EJ,
- Jeanloz R

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Allen MP,
- Tildesley DJ

- ↵
- ↵
- ↵

*P*and

*T*, and heat flux in Earth’s core

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Earth, Atmospheric, and Planetary Sciences