## On the Brink of Melting

The considerable pressures and temperatures of Earth's iron-rich inner core make it a challenge to compare measurements made in experimental systems with observed seismic data. Computational simulations of core materials may reconcile any apparent differences. **Martorell et al.** (p. 466, published online 10 October) used ab initio simulations to predict the elastic properties of iron at core pressures. As temperatures approached the melting point of pure iron, the material was predicted to weaken to the point that seismic waves would be slowed considerably. An inner core with a small percentage of light elements like oxygen and silicon near its melting temperature would correspond well with measured seismic velocities.

## Abstract

The observed shear-wave velocity* V*_{S} in Earth’s core is much lower than expected from mineralogical models derived from both calculations and experiments. A number of explanations have been proposed, but none sufficiently explain the seismological observations. Using ab initio molecular dynamics simulations, we obtained the elastic properties of hexagonal close-packed iron (hcp-Fe) at 360 gigapascals up to its melting temperature *T*_{m}. We found that Fe shows a strong nonlinear shear weakening just before melting (when *T/T*_{m} > 0.96), with a corresponding reduction in *V*_{S}. Because temperatures range from *T/T*_{m} = 1 at the inner-outer core boundary to *T/T*_{m} ≈ 0.99 at the center, this strong nonlinear effect on *V*_{S} should occur in the inner core, providing a compelling explanation for the low *V*_{S} observed.

Earth’s inner core is predominantly made of iron (Fe), but it is commonly assumed to contain 5 to 10% Ni (*1*) and also light elements such as Si, C, and S, ~2 to 3 weight percent in total (*1*, *2*). Seismic wave velocities through the inner core are known, but at present, seismological and mineralogical models for the inner core do not agree (*3*–*9*). A major discrepancy between the observed seismic data and current mineralogical models derived from ab initio calculations is that these mineralogical models predict a shear-wave velocity *V*_{S} that is up to 30% greater than the seismically observed values (*4*, *9*, *10*). The addition of small quantities of Ni under these conditions does not reduce *V*_{S} by a sufficient amount to explain this (*9*), and although the effect of light elements on the velocities of Fe is not totally clear at inner-core conditions (*11*, *12*), all studies show that light-element effects are too small [<5% in *V*_{S} for 7% molar fraction in Si at 5000 K and 13,000 kg m^{−3} (*11*)] to solve the discrepancy.

Another possible cause of the discrepancy between mineralogical models and seismic data is that the elastic constants of Fe may soften drastically and nonlinearly very near to its melting point *T*_{m}, as has been observed in other metals. For instance, the shear modulus of Sn has been experimentally and theoretically shown to decrease by more than 50% at temperatures within ~1% of its melting point (*13*, *14*). According to ab initio simulations, the melting point of pure Fe at the conditions of the inner core is in the range 6200 to 6900 K (*15*–*17*) according to phase coexistence calculations (solid and liquid), with upper limit estimates up to 7500 K (*18*) when only the solid phase is heated until melting. The highest temperature for which the elastic properties of hcp-Fe have been obtained computationally is 6000 K (*5*); however, relative to the melting point of this simulation, *T/T*_{m} is probably ~0.8 and so this temperature is likely to be too low to reveal any strong elastic shear weakening just below melting.

To examine whether premelting effects in the elasticity of Fe can resolve the discrepancy in *V*_{S} between mineralogy and seismology, we simulated the effect of temperature on *V*_{S} at 360 GPa for hcp-Fe up to its melting point. We performed periodic ab initio calculations based on density functional theory (DFT) derived from quantum mechanics, coupled with molecular dynamics to obtain the elastic properties of hcp-Fe at finite temperatures (*19*).

Simulations of hcp-Fe at 360 GPa and at temperatures of 6600, 7000, 7250, and 7340 K (*19*), together with previous simulations at the same pressure and lower temperatures (*9*), indicated that the behavior of the elastic constants up to ~6600 K (Fig. 1) is very similar to that found in our earlier work on hcp-Fe at ~315 GPa and temperatures up to 5500 K (*20*). Specifically, elastic constants [defined as the ratio of applied stress on a material to the strain produced (*19*)] c_{11}, c_{33}, and c_{44} decreased with temperature, and c_{12} and c_{13} slightly increased (Fig. 1). However, one important difference between the simulations at 360 GPa and 315 GPa is that at higher pressure c_{33} is always larger than c_{11}; this finding suggests a large pressure dependence for the c_{11}-c_{33} crossover.

Above 6600 K, our calculations show that all of the elastic constants decreased with temperature, with some of them displaying very strong temperature dependence (Fig. 1). In particular, c_{44}, c_{12}, and c_{11} dropped by 46%, 19%, and 32%, respectively, from 7000 to 7340 K. This pronounced drop for a temperature increase of only 340 K indicates that the calculations above 7000 K are approaching the melting point of the simulated system. Analysis of the radial distribution functions and the root-mean-square displacements of the atoms, however, confirmed that the system remained completely solid during the simulation at 7340 K (fig. S1). A simulation at 8000 K melted completely after 16 ps (*19*).

The temperature dependence of the shear modulus (*G*) reveals an almost linear decrease up to 7000 K, followed by an abrupt drop beyond this point (fig. S2). Most models for *G* versus *T* describe only the linear region [e.g., the mechanical threshold stress model or the Steinberg-Cochran-Guinan model (*21*, *22*)] and do not describe its behavior close to the melting temperature. For this reason, Nadal and Le Poac [NP (*13*)] implemented a new model, based on Lindemann melting theory, that accounts for both the linear region and the region close to the melting temperature. Using this model (*19*), we obtained a Lindemann coefficient of *f* = 0.112 and a melting temperature of 7350 K for hcp-Fe at 360 GPa. The *f* coefficient is a material-dependent parameter and is normally between 0.1 and 0.3 (*23*); hence, our value falls in a reasonable range.

We note that the melting temperature obtained in this way is ~850 K higher than that expected from previous ab initio simulations (*16*, *17*), which used the phase coexistence method. This reflects the fact that the goal of the present work is not to obtain an accurate estimate for the melting temperature of hcp-Fe at 360 GPa, but rather to investigate the behavior of its elastic constants (and therefore the seismic velocities) very close to melting. To obtain the elastic constants, our simulations needed to be performed on a system with no preexisting surface or defects (such as the solid-liquid interface required for the phase coexistence approach), and it is well known that melting temperatures obtained in such a homogeneous system (mechanical melting) are substantially higher than the true thermodynamic (or heterogeneous) melting temperature (*24*, *25*). The melting temperature obtained here using the Nadal–Le Poac model is about 15% higher than that obtained using free energies or phase coexistence methods and is in accord with previous work showing that homogeneous melting temperatures are about 20% higher than the heterogeneous melting temperatures (*26*–*31*).

Although the strong decrease in elastic moduli in Fe observed here is seen close to the homogeneous melting temperature, there is good evidence that this also happens in a real heterogeneous sample. First, the experimentally measured elastic constants of Sn show a strong weakening at *T/T*_{m} ≈ 0.99, where *T*_{m} is the true heterogeneous melting temperature. Second, this decrease has also been observed in atomistic simulations in body-centered cubic (bcc) vanadium at *T/T*_{m} ≈ 0.98 (*24*). Third, the strong elastic weakening is associated with a rapid increase in defects (defined as over- or undercoordinated atoms), and this occurs at both the homogeneous and heterogeneous melting temperatures (*27*–*31*). The only difference in the heterogeneous case is that surface and preexisting defects can propagate into the bulk at a lower temperature than in a homogeneous solid. The number of atomic defects in our simulation jumped from 34% at 7340 K to 70% at 8000 K (*19*) (fig. S3). This is in very good agreement with previous simulations on much larger systems (*27*–*31*). Thus, it is relative temperature (*T/T*_{m}) rather than absolute temperature that is important.

Both compressional-wave (*V*_{P}) and *V*_{S} velocities decreased almost linearly with temperature up to ~7000 K at 360 GPa, with a substantial drop beyond this point (Fig. 2). The temperatures at which the velocities from the Preliminary Reference Earth Model (PREM) and our NP-like (like NP model but for velocities instead of *G*) model agree (7130 K for *V*_{P} and 7250 K for *V*_{S}) are at *T/T*_{m} values of 0.971 and 0.988, respectively, relative to the melting temperature of the simulation. Using an adiabatic geotherm (*32*), we find that the center of Earth should be ~200 K hotter than the inner-core boundary (ICB), whereas the melting line at the center of the inner core is 280 K above the temperature at the ICB (*17*); this leads to a value at the center of the inner core of *T/T*_{m} = 0.988. Thus, the core does indeed lie in a range of *T/T*_{m} where the velocities might be expected to be strongly decreased near melting.

Our results show that *V*_{P} and *V*_{S} for the inner core can be fitted with pure Fe for a physically sensible value of *T/T*_{m}. However, our simulated density of pure Fe is ~3% too high, and so the presence of light elements is still required to match inner-core values (table S1), but we would expect Fe with a few percent light elements to also show a strong shear softening near the melting temperature. If we assume that the light elements reduce the melting temperature of the Fe alloy, then the softening will occur at lower temperatures than in pure Fe, putting it in a more reasonable range of likely core temperatures. However, further investigations into multicomponent systems are essential to fully understand their effect on the elastic properties of the core. Overall, our results demonstrate that the inner core is likely to be in the strongly nonlinear regime; hence, there is no need to invoke special circumstances such as strong anelasticity, partial melts, or combinations of crystalline phases in order to match the observed seismic velocities and densities of the inner core.

## Supplementary Materials

www.sciencemag.org/content/342/6157/466/suppl/DC1

Materials and Methods

Supplementary Text

Figs. S1 to S3

Table S1

References (*33*–*43*)

## References and Notes

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
See supplementary materials on
*Science*Online. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
**Acknowledgments:**Supported by Natural Environment Research Council grant NE/H003975/1 (L.V.). Calculations were performed in the HECTOR supercomputer facility. Computer code VASP is available at www.vasp.at. The data presented in this paper are given in the main text and in the supplementary materials. B.M. performed research, analyzed data, and wrote the paper. L.V. designed research, analyzed data, and wrote the paper. J.B. analyzed data and wrote the paper. I.G.W. designed the research, analyzed data, and wrote the paper.