Thermal conductivity of Fe-Si alloys and thermal stratification in Earth's core.

Light elements in Earth's core play a key role in driving convection and influencing geodynamics, both of which are crucial to the geodynamo. However, the thermal transport properties of iron alloys at high-pressure and -temperature conditions remain uncertain. Here we investigate the transport properties of solid hexagonal close-packed and liquid Fe-Si alloys with 4.3 and 9.0 wt % Si at high pressure and temperature using laser-heated diamond anvil cell experiments and first-principles molecular dynamics and dynamical mean field theory calculations. In contrast to the case of Fe, Si impurity scattering gradually dominates the total scattering in Fe-Si alloys with increasing Si concentration, leading to temperature independence of the resistivity and less electron-electron contribution to the conductivity in Fe-9Si. Our results show a thermal conductivity of ∼100 to 110 W⋅m-1⋅K-1 for liquid Fe-9Si near the topmost outer core. If Earth's core consists of a large amount of silicon (e.g., > 4.3 wt %) with such a high thermal conductivity, a subadiabatic heat flow across the core-mantle boundary is likely, leaving a 400- to 500-km-deep thermally stratified layer below the core-mantle boundary, and challenges proposed thermal convection in Fe-Si liquid outer core.

The geodynamo of Earth’s core is thought to be mainly driven by compositional (chemical) convection associated with the crystallization and light-element release of the inner core as well as thermal convection driven by a superadiabatic heat flow across the core–mantle boundary (CMB). The relative importance of these energy sources to the geodynamo, however, remains uncertain (1). The magnitudes of these energy sources can change throughout the evolution of the planet. The thermal gradient across the CMB can be constrained from both heat flow of the core and mantle, where a subadiabatic heat flow out of the core may hinder thermal convection and cause a thermally stratified layer at the top of the outer core (2). A global nonadiabatic structure at the top of the core has been inferred from seismic observations and geomagnetic fluctuations (3, 4), where the mechanisms for the origin rely on accurate determinations of the CMB heat flow and the core conductivity. Based on seismological observations and high-pressure and -temperature (P-T) mineral physics results, Earth’s outer and inner core are mainly composed of Fe (∼85 wt %) alloyed with Ni (∼5 wt %) and ∼8 to 10 wt % and 4–5 wt % of light elements, respectively, such as Si, O, S, C, and H (5–10). The effects of the candidate light elements on the electrical resistivity (ρ) and thermal conductivity (κ) of iron and their partitioning between the inner and outer core at relevant P-T conditions are thus of great importance for understanding the thermal state of the core as well as the generation and evolution of Earth’s magnetic field (2, 9, 11, 12). The thermal conductivity of the constituent core alloy controls the heat flow of the core, while the electrical resistivity of the constituent Fe alloy determines the ohmic dissipation rate of the magnetic field.

Extensive studies on iron’s transport properties have been conducted via experiments and calculations (e.g., refs. 13–21), and recent studies report a thermal conductivity of ∼100 W⋅m−1⋅K−1 at conditions near the CMB (22, 23). Such a high thermal conductivity reduces the amount of heat that can be transported by convective flow (11) and raises a question as to what powered the convection prior to inner core growth over Earth’s history [the so-called new core paradox (24)]. Thus far, several hypotheses have been proposed to reconcile this paradox, including a possible large conductivity reduction due to nickel and light elements (25–28), a rapid core cooling rate (29), or exsolution of chemically saturated species from the core to the lowermost mantle, such as MgO, SiO2, or FeO (e.g., refs. 30–32). The general consensus is that incorporation of light element(s) depresses high P-T thermal conductivity of iron by impurity scattering (12); this effect was assumed in our previous modeling of the high P-T transport properties of Fe-Ni alloyed with 1.8 wt % Si (25). The lowered thermal conductivity implies that thermal convection is easier to maintain. The rapid core cooling model would imply a young inner core and requires a hidden core heat source, such as radioactivity, which is not supported by geochemical evidence (29). The exsolution mechanism would offer an additional energy source to drive an early compositionally driven geodynamo (32), although some experiments find exsolution to be unlikely (33). The viability of each of these scenarios depends sensitively on the transport properties of iron alloyed with a significant amount of light element(s) (∼8 to 10 wt %) at core P-T conditions. Information on these electrical and thermal transport properties of iron alloys remain uncertain due to the sparsity of experimental and theoretical data.

Here we focus on the geodynamic consequences of the transport properties of iron alloyed with 4 to 10 wt % silicon, which is considered to be one of the major light element candidates in the Earth’s core due to its geo- and cosmochemical abundance (5), high solubility in solid and liquid iron (34), and isotopic evidence (35). Fe-Si alloys have been the subject of previous studies focused on understanding the structural and physical properties of the core material, including its high P-T phase diagram (36–39), elasticity (40–44), melting behavior (36, 45, 46), and transport properties (25, 47–49). The observed density discontinuity of ∼4 to 5% across the inner-core boundary (ICB) indicates that excess light elements partition into the outer core during inner-core solidification (6, 50). We should note that the concentration of Si in the core remains uncertain. While some experiments have shown that Fe alloyed with ∼9 wt % Si can satisfy the density profile of the outer core and Fe alloyed with ∼4 wt % Si for the inner core, respectively (37, 40, 41, 51, 52), other studies indicate that a dominant Si light alloying component is unable to reproduce both the density and sound velocity distribution in the outer core (53, 54).

High P-T diamond anvil cell (DAC) experiments had been previously conducted to constrain the electrical and thermal conductivity of Fe-Si alloys (28, 47, 55, 56), specifically their T-dependent resistivity and thermal conductivity at core pressures. The thermal conductivity of Fe-8 wt % Si (hereafter Fe-8Si) was measured using a high-P ultrafast optical pump probe and high P-T flash-heating methods (28). The results showed that 8 wt % silicon in solid hexagonal close-packed (hcp) Fe can strongly reduce the conductivity of pure iron by a factor of ∼2, i.e., giving ∼20 W⋅m−1⋅K−1 at ∼132 GPa and 3,000 K. However, the electrical resistivity of solid Fe-6.5Si at ∼99 GPa and 2,000 K was recently measured to be ∼73 µΩ⋅cm (56), which is higher than that of pure iron (22) by ∼60% at comparable conditions. The results imply a thermal conductivity of ∼66 W⋅m−1⋅K−1 using the Wiedemann–Franz law (TL = ρ) assuming an ideal Sommerfeld Lorentz number (L = L0: 2.44 × 10−8 W⋅m−1⋅K−2). Meanwhile, another recent study reported a moderate thermal conductivity of 50 to 70 W⋅m−1⋅K−1 for an Fe-5Ni-8Si alloy near CMB P-T conditions (∼140 GPa and 4,000 K) modeled from the measured resistivity of Fe-10Ni and Fe-1.8Si alloys using the four-probe van der Pauw method in laser-heated DACs (25). The results on Fe-10Ni and Fe-1.8Si alloys reveal a linear relationship between resistivity and temperature at a given high pressure, which is very similar to that of hcp Fe (22), over the range of measurements. In contrast, density functional theory (DFT)-based molecular dynamics simulations predict a small negative T dependence of the resistivity at high pressure when liquid Fe is alloyed with a significant amount of light elements (e.g., ∼13 wt % Si) (27). These experimental and computational results raise the possibility that the high P-T thermal transport behavior and its temperature dependence in Fe-Si alloys with a few wt % Si (e.g., 2 wt %) and a larger wt % Si (e.g., 8 to 10 wt %) can be quite different, making it difficult to evaluate the light element effects on the energetics of the core.

In this study, we directly measured the electrical resistivities of polycrystalline hcp Fe-4.3 wt % Si (Fe-4.3Si, or Fe0.92Si0.08) and Fe-9 wt % Si (Fe-9Si, or Fe0.84Si0.16) alloys to ∼136 GPa and 3,000 K. We also computed the electrical resistivity and thermal conductivity of these Fe-Si alloys in solid and liquid phases using first-principles molecular dynamics (FPMD) and dynamical mean field theory (DMFT) calculations. The calculations include contributions from scattering off of Si as well as both electron–phonon (e-ph) and electron–electron (e-e) scattering. Our results are used to evaluate the Si impurity effects on the transport properties of Fe-Si alloy at P-T conditions of the topmost outer core. Assuming Si is the sole light element in the core, our results are used to constrain core thermal conductivity, which is in turn used to assess core heat flux, thermal state, and energy sources driving the geodynamo through geodynamical modeling.

Results

The electrical resistivity of polycrystalline Fe-4.3Si and Fe-9Si alloys was measured in flat-top laser-heated DACs using the four-probe van der Pauw method (see and ). Initial measurements of Fe-9Si up to ∼136 GPa at room temperature () are consistent with previous studies on Fe-Si alloys up to ∼90 GPa (47, 55). The well-known body-centered cubic (bcc)–hcp (α-ε) phase transition observed at ∼15 GPa causes an increase in the resistivity (). With further increase in pressure, the resistivity of hcp Fe-9Si gradually diminishes and the pressure coefficient of resistivity decreases. The hcp-structured Fe-Si alloys have higher resistivities than pure hcp Fe due to the impurity scattering. We find that the resistivity in the Fe-Si alloys increases quasilinearly with increasing silicon concentration up to ∼9.0 wt % at room temperature (47).

At ∼100, 120, and 136 GPa, electrical resistivities of the polycrystalline Fe-Si alloys were measured with increasing temperature up to ∼3,000 K (Figs. 1 and 2 and ). In situ powder X-ray diffraction (XRD) patterns of the samples show well-distributed intensity, indicating that the polycrystalline samples in the hcp structure did not experience significant lattice-preferred orientation due to possible deviatoric stresses and/or high temperature-induced recrystallization over the investigated P-T conditions (). The use of strong cubic boron nitride (cBN) as the gasket insert and two flat layers of SiO2 as the pressure medium in our experiments provided that the samples maintained their initial crossed-heart shapes and even thicknesses that are critical for the success of the four-probe van der Pauw method. The laser heating experiments were conducted up to ∼2,500 to 3,000 K, which is below the expected fast-recrystallization threshold needed to avoid the possible influence of recrystallization on conductivity () (57, 58). The hcp phase was observed in Fe-4.3Si and Fe-9Si alloys at these high P-T conditions, consistent with the previously determined stability field of the hcp phase () (59). Analysis of the quenched Fe-9Si samples after the high P-T measurements showed compositional homogeneity of Fe with 9.6(1.8) wt % Si (see ). These results indicate that the samples were not chemically contaminated and the hcp Fe-9Si alloy was the predominant phase during the high P-T measurements ().

Fig. 1.

Temperature-dependent electrical resistivity of hcp Fe and Fe-Si alloys at ∼100 GPa. The experimentally measured resistivities of hcp Fe-4.3Si (solid red squares) and Fe-9Si alloys (solid and open blue squares) up to ∼2,500 K at pressures of 90 to 105 GPa show a decreasing temperature dependence with increasing Si concentration. The resistivities of Fe-4.3Si and Fe-9Si alloys are compared with our DFT/DMFT calculations, where the open diamonds represent the DFT/DMFT calculations contributed by e-ph, e-e, and impurity scattering, and the semiopen diamonds represent contributions only by e-ph and impurity scattering using DFT. The data for hcp Fe-4.3Si are generally consistent with the literature results of hcp Fe-4Si (open black triangles) and Fe-6.5Si (open black inverted triangles) alloys at ∼99 GPa in an internally heated DAC (56). The open circles are for the resistivities of hcp Fe by DMFT calculations including the contributions from e-ph and e-e scattering (23). The colorrf dashed and dashed-dotted lines represent the resistivity–temperature relation in hcp Fe-4.3Si and Fe-9Si fitted using the Bloch–Grüneisen formula. Literature resistivity results of hcp Fe (22) and Fe-1.8Si alloy (25) are fitted using the Bloch–Grüneisen formula (Fe: black dashed line, Fe-1.8Si: black dotted-line).

Fig. 2.

Measured and calculated electrical resistivities of hcp Fe-4.3Si and Fe-9Si alloys up to ∼3,000 K at 120 to 136 GPa. The open and solid blue squares represent the measured resistivity of hcp Fe-9Si and the solid red squares represent the measured resistivity of hcp Fe-4.3Si. The open diamonds represent the resistivities of the Fe-Si alloys using DFT/DMFT calculations contributed by e-ph, e-e, and impurity scattering, and the semiopen diamonds represent the calculated resistivities of the Fe-Si alloys using DFT calculations contributed only by e-ph and impurity scattering. The measured high P-T resistivity of hcp Fe (22) is also plotted for comparisons (open circles). The black dashed line represents the resistivity of hcp Fe fitted using the Bloch–Grüneisen formula. The colored dashed and dashed-dotted lines represent the resistivity–temperature relation in hcp Fe-4.3Si and Fe-9Si fitted using the Bloch–Grüneisen formula.

Our results show that the resistivity of hcp Fe-4.3Si alloy increases linearly with increasing temperature at high pressures over the measured range (). However, the resistivity-temperature slope in hcp Fe-4.3Si is significantly lower than that of literature values for hcp Fe and Fe-1.8Si alloy (22, 25). Our measured resistivity of hcp Fe-4.3Si is also consistent with the previous results using a four-probe method in an internally heated DAC, where it is located between those of hcp Fe-4Si and Fe-6.5Si up to ∼3,000 K at ∼100 GPa (56) (Fig. 1). For simplicity, we used the Bloch–Grüneisen formula to fit the measured resistivity data ():where ρ is the residual electrical resistivity of Fe-Si alloys when the phonons of the alloys are frozen and ρ, n, and D(V) are fitted constants based on the measured high P-T resistivities. Fitting results show a weakening temperature dependence of the resistivity in hcp Fe-Si alloys with increasing concentrations of silicon. The fitted n parameter in the Bloch–Grüneisen formula that describes the T dependence of the resistivity decreases from 1 to 3 in hcp Fe (22) to 0.6 to 0.7 in hcp Fe-1.8Si (25), then down to 0.15 to 0.18 in hcp Fe-4.3Si at similar P-T conditions (). Furthermore, the resistivity of the Fe-9Si alloy appears to be nearly constant with increasing temperature (Fig. 1). For example, the resistivity of hcp Fe-9Si varies from 81(3) µΩ⋅cm at room temperature to 84(4) µΩ⋅cm at 2,120(90) K at ∼100 GPa, so there is no visible temperature effect on the resistivity within uncertainties. As a result, the Bloch–Grüneisen n parameter is close to 0 in hcp Fe-9Si (Fig. 1 and ). Note that the Bloch–Grüneisen formula is derived only for e-ph scattering, so it is not surprising that the fitted n parameters are anomalous for Fe-Si alloys where impurity scattering dominates. At relatively low temperatures, the resistivity of hcp Fe is much lower than those of Fe-Si alloys due to the impurity scattering; however, they approach each other with increasing temperature due to the lower temperature coefficient of resistivity in the hcp Fe-4.3Si and Fe-9Si alloys (Fig. 1).

We computed the electrical resistivity of the hcp Fe-Si alloys at P-T conditions similar to those in the present experiments using FPMD/DFT/DMFT calculations that include e-ph, chemical impurity, and e-e scatterings (see ). Overall, the calculated high P-T resistivities of hcp Fe-4.3Si and Fe-9Si alloys are consistent with our experimental results (open red and blue diamonds for Fe-4.3Si and Fe-9Si, respectively, Figs. 1 and 2). These calculations do not assume Matthiessen’s rule (that scattering is additive) or the Wiedemann–Franz law. Both experiments and theory consistently show that the T dependence of the resistivity decreases with increasing the concentration of silicon. The calculated resistivity of hcp Fe-Si alloys is also compared with the FPMD/DFT calculations including only e-ph and impurity scattering (semiopen diamonds in Figs. 1 and 2). The results show that the e-e contribution to the resistivity of hcp Fe-4.3Si and Fe-9Si alloys is less than ∼10% at high temperatures and pressures of 90 to 140 GPa (Figs. 1 and 2). Specifically, the e-e contribution increases with increasing temperature but it decreases with increasing Si concentration in hcp Fe-Si alloys at high pressures (). For example, the e-e contribution in hcp Fe-9Si could be less than ∼6% near the CMB conditions (∼4,000 K and 136 GPa; Fig. 2), in contrast to the significant contribution to hcp iron where it is ∼30% ().

Our theoretical results indicate that impurity scattering is the main contribution to the high P-T resistivity of Si-doped Fe alloys (4.3 to 9.0 wt %Si) (Figs. 1 and 2). The resistivities of hcp Fe and Fe-Si alloys with low Si concentrations such as Fe-1.8Si exhibit a significant T dependence (22, 23, 25) because the e-ph and e-e scattering from iron then dominates the total scattering and is T-dependent. With increasing Si concentration in Fe-Si alloys, the T-independent Si impurity scattering gradually takes over the total scattering, leading to T independence and less e-e contribution in the resistivity of Fe-9Si (Figs. 1 and 2). The low T dependence in the resistivity of Fe-Si alloys was predicted and previously explained by the chemically induced “resistivity saturation” in some early DFT studies (47). However, the effect documented here is not resistivity saturation but rather is due to the impurity scattering’s being nearly T-independent and dominant, as well as a breakdown of Matthiessen’s rule.

The electrical resistivity and thermal conductivity of hcp and molten Fe-Si alloys at P-T conditions relevant to the core were further calculated using DFT/DMFT (Fig. 3). We assumed that the Fe-Si alloys remain in the hcp phase near CMB pressure below ∼4,000 K and are molten above 6,000 K, consistent with the determined melting relation of Fe-Si alloys (36, 38, 45). Based on our experimental and theoretical results, the resistivity of hcp Fe-4.3Si and Fe-9Si is estimated to be 80(6) and 90(6) µΩ·cm at 140 GPa and 4,000 K, respectively (Fig. 3). The computations also predict that the resistivity of Fe-Si alloys does not significantly change upon melting above 6,000 K (Fig. 3 and ). The minor effect of melting can be understood as a result of the dominant electronic contribution on the resistivity in both hcp Fe and liquid phase and the fact that under high compression the local structures of liquid and solid are similar. The thermal conductivity of hcp Fe-Si alloys was then determined from both the calculations and the measured resistivity via the Wiedemann–Franz law using our computed Lorentz number. The calculated Lorentz number (L) varies with temperature at the pressure of the topmost core (Fig. 3). At P-T conditions of the CMB (135 to 150 GPa and ∼4,000 K), the calculated Lorentz number of hcp Fe-Si alloys (2.33 to 2.39 × 10−8 W⋅Ω⋅K−2) does not deviate appreciably from the Sommerfeld value L (2.44 × 10−8 W⋅m−1⋅K−2), but they are significantly larger than that of hcp Fe (2.1 × 10−8 W⋅Ω⋅K−2) (22). Our results show that the thermal conductivity of hcp and liquid Fe-Si alloys increases linearly with increasing temperature at the pressure near the CMB (Fig. 3). Further, the theoretical results predict little change on the melting due to electronic contribution’s being dominant in liquid iron, which has a local atomic structure similar to the hcp structure (60).

Fig. 3.

Calculated electrical resistivity (ρ), thermal conductivity (κ), and Lorentz number (L) of Fe-Si alloys in solid and liquid states at the relevant P-T conditions of the topmost outer core. (A) Calculated electrical resistivities of Fe-Si alloys using DFT/DMFT as a function of temperature along isochores. We assume that the Fe-Si alloys are in the solid hcp phase at temperatures below 4,000 K and molten above 6,000 K. Error bar has the size of ρ − ρ, where s = x, y, z, and ρ is the average value over three directions. The anisotropy disappears for the liquid state so the SE bar decreases. (B) Calculated thermal conductivities of Fe-Si alloys as a function of temperature along isochores. The thermal conductivity of the Fe-Si alloys increases almost linearly with increasing temperature. (C) Calculated Lorentz number (L) of the Fe-Si alloys as a function of temperature along isochores. The blue dotted line represents the measured resistivity of hcp Fe-9Si and its thermal conductivity was derived using the Wiedemann–Franz relation. The calculated resistivities and thermal conductivities of liquid Fe-6.7Si by DFT are also plotted for comparisons (open squares) (14). The thermal conductivities of Fe-2Si and Fe-8Si were measured using the transient heating laser technique in LHDAC (28). Pressures (and temperatures) along the isochore in the DMFT computations for hcp Fe-9Si (54.6 bohr3 per atom) are 136 GPa (2,400 K), 150 GPa (4,000 K), 183 GPa (6,000 K), 203 GPa (8,000 K), respectively; for hcp Fe-4.3Si (55.6 bohr3 per atom) they are 122 GPa (2,500 K), 135 GPa (4,000 K), 169 GPa (6,000 K), 188 GPa (8,000 K), respectively; for hcp Fe (54.9 bohr3 per atom) they are 140 GPa (2,500 K), 145 GPa (3,000 K), and 149 GPa (3,500 K), respectively (23).

The thermal conductivity of hcp Fe-9Si is found to be ∼100 to 110 W⋅m−1⋅K−1 at ∼140 GPa and 4,000 K based on both experimental and theoretical results, which is unexpectedly close to that of pure hcp Fe at similar conditions (100 ± 10 W⋅m−1⋅K−1) (22). Therefore, adding 9 wt % Si into iron does not reduce its thermal conductivity near CMB conditions. Surprisingly, the thermal conductivity of hcp/liquid Fe-9Si alloy is likely larger than that of Fe with further increasing temperature above 4,000 K (Fig. 3). The incorporation of a large amount of Si (e.g., 9 wt %) in Fe makes Fe-Si alloy an even better thermal conductor at very high temperatures because the thermal conductivity continuously increases but resistivity is almost constant with increasing temperature. Although we did not consider the phonon contribution to the thermal conductivity of Fe-Si alloys, it has been previously evaluated to be ∼3 W⋅m−1⋅K−1 in Fe alloy (61). The result indicated a negligible contribution from phonons compared to the electronic contribution of 100 to 110 W⋅m−1⋅K−1 in Fe-Si alloys.

Compared with the recently modeled thermal conductivity of Fe-Si alloys obtained from high P-T measurements of heat-pulse propagation in samples in laser-heated DACs (28), our estimated thermal conductivity of Fe-(4.3 to 9.0)Si alloy (∼78 W⋅m−1⋅K−1) at 120 to 140 GPa and ∼2,500 K are slightly higher than that of Fe-2Si (∼64 ± 16 W⋅m−1⋅K−1, inverted triangle, Fig. 3), but approximately four times larger than that of Fe-8Si alloys (∼20 ± 10 W⋅m−1⋅K−1, triangle, Fig. 3) from the previous experiments (Fig. 3). The large discrepancy of high P-T thermal conductivity in Fe-8Si alloy obtained from the transient heating laser technique may arise from the large temperature uncertainty in the latter, leading to an indistinguishable temperature effect on the thermal conductivity in Fe-8Si (62). Another possibility is that Fe-8Si underwent a phase transition from hcp to mixed hcp with Si-enriched bcc in the experiments of ref. 28 because of higher temperatures (e.g., up to ∼3,800 K), where the hcp–bcc mixed phase has lower thermal conductivity because of the lower conductivity in the bcc phase. Since our calculations predict a negligible effect of melting on these transport properties of Fe-Si alloys, our determined resistivity and thermal conductivity of liquid Fe-4.3Si and Fe-9Si alloys at the P-T conditions of the outer core are also generally consistent with those of liquid Fe-6.7Si previously predicted using DFT techniques (Fig. 3 and ) (14, 27). Altogether, the results indicate that the e-ph and e-e scattering contribute less than the Si impurity scattering to the transport properties in these high P-T Fe-Si alloys.

Discussion

The experimental and theoretical results described above have important geodynamic implications. Although some distributions of other light elements are present in the core, we assume Si is dominant for purposes of the analysis (35). The high thermal conductivity of ∼100 to 110 W⋅m−1⋅K−1 of the Fe-9Si alloy at P-T conditions corresponding to the top of the liquid core (κ) indicates an adiabatic heat flow of ∼15 TW (, r is the diameter of the outer core surface), assuming the temperature gradient () as ∼0.9 to 1.0 K/km at the top of the outer core (11, 63). The present-day CMB heat flow across the lowermost mantle (Q) is constrained to be ∼10 to 12 TW based on the estimated thermal conductivity (∼10 W⋅m−1⋅K−1) of the mantle minerals and temperature gradient in the lowermost mantle (64, 65). Thus, if ∼9.0 wt % Si is the sole light element concentration in the core, the adiabatic heat conduction limit of the core (Q) is most likely higher than Q and implies a subadiabatic temperature gradient in the core, thus inhibiting thermal convection.

The geodynamo is maintained by fluid convection in the outer core (66). It is believed that the fluid convection is mainly driven by thermal and/or compositional buoyancy, where the former requires a higher Q than the Q based on the Schwarzschild criterion (67), and the latter occurs with light-element release from inner-core solidification (68). The high thermal conductivity in liquid Fe-9Si makes the core Q > Q and thus hinders thermal convection. However, chemical buoyancy released at the ICB by inner core growth could still drive convection at depth. As a result, the present-day geodynamo would be mainly driven by compositional convection associated with the release of incompatible light elements and latent heat upon inner-core solidification (Fig. 4). Additionally, the present study shows a linearly increasing thermal conductivity in Fe-9Si alloy with temperature, which causes Q increased at higher temperatures. Therefore, the early hotter liquid core had a higher thermal conductivity, making it even harder to drive a thermal dynamo.

Fig. 4.

Thermal stratification at the topmost outer core and energy sources for geodynamo influenced by the presence of light alloying elements. (A) Thermal stratification thickness (z) at the top of the outer core as a function of the CMB heat flow (Q) for an adiabatic heat flow (Q) of ∼15 TW out of the core (with the thermal conductivity of ∼100 W⋅m−1⋅K−1 for liquid Fe-9Si alloy). The thermal stratification may have a thickness in the order of a few hundred kilometers depending on the magnitude of the Q, for example, taking the Q of 10 to 13 TW gives a thickness of 700 to 300 km. The presence of the stratified layer may explain the radial variation via the seismic observation in this region. The recent outer core wavespeed models AK135 (69), IASP91 (70), and KHOCQ10 (4) compared to preliminary reference Earth model (PREM) (V) vs. stratified thickness are also plotted for comparison. (B) The schematics shows that the topmost outer core is thermally stratified due to the subadiabatic heat flux across the CMB (Q < Q) as shown in A. The radial convection is prevented in the stratified layer. The subadiabatic CMB hinders thermal convection; therefore, the geodynamo is dominantly driven by compositional (chemical) convection below the stratified layer that is related to light element(s) partitioning during the inner core crystallization in B. With the inner-core nucleation, the excess light elements (e.g., O or S) continuously release from the solid inner core to the liquid outer core at the ICB. The outer core contains ∼8 to 10 wt % light elements and the CMB temperature is ∼4,000 K.

Our results indicate that the uppermost core is subadiabatic and stable to thermal convection (11), which may result in thermal stratification at the top of the core unless the present-day Q is underestimated (Fig. 4). As a result, convection may be thermally stratified below the CMB and convectively unstable above the ICB in the core, with a transition radius r, or the thickness z = R − r (R is the radius of the core) that we refer to as the depth of a thermally stratified layer (see ). To examine this effect of the thermal stratification in the outer core, we characterize the thermal inversion by investigating the Brunt–Väisälä frequency (N) (see and ) and modeling the light-element concentration perturbation, composition buoyancy, and total thermal and compositional buoyancy as a function of the radius (). Convection occurs where the Brunt–Väisälä buoyancy frequency N is imaginary, where the N is defined byg is the gravity, ρ is the density, r is the radius, and dρ/dr is the radial derivation (3). Therefore, N = 0 corresponds to the location of a transition from a convective to nonconvective region. In the model, the computed thermal and compositional buoyancy profiles depend on Q and the total buoyancy gradient crosses zero at the stratification radius r (). The thickness of the stratified layer (z) is then computed as a function of the Q assuming our determined Q of ∼15 TW (see ) (Fig. 4). Taking the Q of ∼12 TW (64, 65), the thickness of the stratified layer is constrained to be ∼400 to 500 km. The estimated thermal stratification layer overall overlaps the regions of the observed seismological velocity gradient at the top 100 to 700 km of the outer core from SmKS (2 ≤ m ≤ 4) analyses (4, 69, 70) (Fig. 4), but a recent seismic model requires no such an anomalous layer (71). We should also note that the thickness of the thermal stratification layer may show regional heterogeneity due to the proposed CMB heat-flux variations (72). Further, the thermal stratification combined with possible chemical stratification likely causes the slight reduction of the seismic velocity atop Earth’s outer core (73, 74).

Besides Si, the other leading candidate light elements present in the core are S, O, C, and H (5). At high P-T deep-Earth conditions, Si and S substitute for iron to form Fe alloys, while O, C, and H can be incorporated into iron in its interstitial position. Hydrogen has long been considered as a light element (8) and in only small concentrations can have a significant effect on the core’s properties (75). Recent studies indicate the presence of the other light elements chemically limit the H abundance (76). Different alloying may cause different temperature responses of the resistivity at high pressure according to a recent prediction (27), so electrical and thermal conductivity for other Fe-light element alloys at high P-T conditions need to be further evaluated by both experiments and theories. On the other hand, silicon is readily soluble in both solid and liquid iron at core P-T conditions (46). Thus, some other (multiple) light elements with different partitioning behaviors across the ICB [e.g., O and C (77, 78)] likely coexist in the core, which may partition between liquid and solid inner core to compositionally drive convection. Evidence for the presence of the geomagnetic field goes back to at least ∼3.5 Ga (79), which may predate inner-core nucleation. This observation suggests that an Fe-9Si-rich core may have prevented an early thermal geodynamo before inner-core nucleation, raising again the core paradox (24). We therefore suggest that Si may not be the sole light element alloyed with iron in the core considering the negligible effect of Ni on the conductivity of iron (25), and a significant amount of other light elements discussed may coexist in the liquid outer core. These light elements may have precipitated from the early molten core and accumulated in the lowermost mantle with the cooling of the core, consistent with a recent simulation indicates that 1 to 2 wt % light element or oxide precipitation may have been sufficient to power the early dynamo (32). Future studies on the effects of all candidate light elements on the thermal transport properties of Fe-Ni alloy as well as Q and Q values at relevant high P-T conditions are critically needed to reliably evaluate the possible occurrence and thickness of the thermal stratification layer at the top of the outer core.

In summary, we have performed electrical resistivity measurements coupled with in situ XRD of Fe-Si alloys with 4.3 to 9.0 wt % Si up to ∼140 GPa and 3,000 K in laser-heated DACs. The temperature-dependent resistivity of Fe-Si alloys has been determined by experiments and first-principles simulations. The overall consistent agreement between experiments and theory shows a greatly decreasing temperature coefficient of resistivity in Fe-Si alloys with increasing Si concentration. Electron–phonon scattering and impurity scattering dominate the contributions to the resistivity of the Fe-Si alloys, whereas electron–electron scattering is much weaker (e.g., compared to pure iron). The thermal conductivity of liquid Fe-9Si is constrained to be ∼100 to 110 W⋅m−1⋅K−1 near CMB conditions. Such a high thermal conductivity hinders thermal convection and requires compositional convection in an Fe-Si core. If a thermally stratified layer extends 400 to 500 km below the CMB, as has been inferred by several seismic models (Fig. 4), our model that uses Si as the major light element constrains the CMB heat flow to be 10 to 12 TW. Such a thick stratified layer, however, is in contrast to dynamo models that Earth’s magnetic field is more consistent with a layer of 100 km thick, or less (80), which would imply a CMB heat flow of ∼15 TW. Further studies of the transport properties in other Fe alloys (e.g., Fe-O/S/C/H and ternary Fe-Si-O/S/C/H alloys, etc.) at core P-T conditions could shed new light on the mechanism of the core convection over Earth’s history.

Materials and Methods

Starting Materials.

We used polycrystalline Fe-4.3Si and Fe-9Si alloys with bcc phase as starting materials. Fe-9Si alloy was purchased from Goodfellow Corporation, the same one in ref. 39. We synthesized Fe-4.3Si alloy from a mixture of iron and silicon powder with a starting mass ratio of 95.7 and 4.3 by an arc-melting method in an arc furnace full of pure argon atmosphere to avoid any oxidation (43). The synthesized Fe-4.3Si sample was well-characterized using an electron probe microanalyzer (Shimadzu, EPMA-1720H) at Sichuan University, China. The back-scattered electron image and electron probe microanalyzer analyses show a chemically homogeneous composition in Fe-4.3Si alloy (). The samples were compressed to thin foils 1.5 to 2 μm thick using a pair of diamond anvils with 600-μm culet diameter. The polycrystalline foils were shaped into a uniform Greek cross shape with a diameter of ∼6 μm at the center () by a focused ion beam (FIB, FEI Versa 3D). The samples were loaded in DACs with 100- to 300-μm beveled culets with a beveled angle of 9°. Two dried SiO2 layers were used as thermal insulator (). We compressed the Fe-Si alloys to targeted pressures in DACs and heated them using a flat-top double-sided laser heating system ().

Electrical Resistivity Measurements in Laser-Heated DACs.

High P-T electrical resistivity and in situ XRD on Fe-4.3Si and Fe-9Si alloys were performed at the 13-IDD beamline station, GeoSoilEnviro Consortium for Advanced Radiation Sources, Advanced Photon Source. The resistivities of Fe-Si alloys were measured at high pressures from ∼90 to 136 GPa and temperatures up to ∼3,000 K by a van der Pauw method, where a multimeter source (Keithley model 6221) and an ultralow voltmeter (Keithley model 2182A) were used. The high P-T phase diagrams of the Fe-Si alloys have been well-determined by previous laser-heated DAC experiments (38, 39, 59) (). The resistance (R) at high P-T was calculated by Ohm’s law according to the averaged voltage. The resistivity was obtained from the equation of ρ = ρ(R/R)(V/V)1/3, where ρ is the resistivity at room temperature and high P, R, and V are the resistance and volume at high pressure and room temperature and R and V are the measured resistance and voltage at high P-T conditions. (). The pressure in Fe-Si samples was obtained from the measured lattice parameters by in situ XRD, which were compared with the previously determined thermal equation of states in Fe-4Si (59) and Fe-9Si (37, 39) alloys, respectively. The chemical composition of the quenched Fe-Si samples was analyzed using scanning electron microscopy and energy-dispersive X-ray spectroscopy in FIB (). The details of the experimental process can be found in previous literature (22, 25).

Computational Details for Electrical and Thermal Transport Properties.

FPMD simulations for Fe-Si were performed using Quantum ESPRESSO (81). All the FPMD calculations used the ground-state Perdew–Burke–Ernzerhof (82) generalized gradient approximation for the exchange-correlation functional. Supercells containing 216 atoms in the hcp structure were prepared for random placement of the Si atoms in the hcp lattice with the tool Alloy Theoretic Automated Toolkit (ATAT) using the Special Quasi-random Structures (SQS) algorithm (83, 84). The lattice constants from the XRD measurements in this study were used. We used scalar-relativistic Garrity–Bennett–Rabe–Vanderbilt (GBRV) ultrasoft pseudopotentials for Fe and Si, which were generated using smaller cutoff radii. The pseudopotentials were tested against linearized augmented plane wave (85) and standard GBRV pseudopotentials. The plane-wave cutoff energy was 40 Ry. For each temperature, the ionic temperature was regulated by the stochastic-velocity rescaling thermostat (86) in the N-V-T ensemble and identical electronic temperatures were set with the Fermi–Dirac smearing function. The Brillouin zone is sampled at the Γ-point in the FPMD. The FPMD step was 1 fs and the simulations were longer than 10 ps to ensure the equilibrium state was reached. All bands with occupations larger than 10−5 were included. Macroscopic quantities, such as temperature and pressure, were obtained by averaging over 3,000 or more equilibrium steps.

In the transport calculations of Fe-Si alloys, we employed relativistic KKR-DMFT as implemented in a modified spin-polarized relativistic KKR (SPRKKR) package (87, 88) to selected snapshots. These snapshots were separated from each other by more than the velocity correlation time under each condition. For each one, the calculations were separated into three steps: 1) DFT self-consistent-field (SCF) calculation, 2) SCF calculation at the DMFT level, and 3) transport calculations with potentials from 1 and 2 using the Kubo–Greenwood method (see ) (89, 90). The converged potential of selected ionic configurations was then used as the starting potential for the next step. Additional details are presented in .