07_Introduction (part 2: Geochemistry) - Yiwei666/08_computional-chemistry-learning-materials- GitHub Wiki
1. Trace element partitioning between silicate melts – A molecular dynamics approach
1. Introduction
The use of trace elements as tracers of magmatic processes in the Earth’s interior is a standard tool on the geochemist’s workbench and as such, these elements have received much attention throughout the last decades (Goldschmidt, 1937, Ringwood, 1966, Blundy and Wood, 2003, Shaw, 2006). The distribution of an element between phases involved in a geochemical reaction is conventionally quantified with the partition coefficient D. D is defined as the concentration ratio of element X between two phases and (e.g. a mineral and a melt) (Beattie et al., 1993): (1)
Knowledge of these partition coefficients is essential for quantitative modelling of mineral-melt reactions and the corresponding experimental studies are numerous (e.g. Green, 1994, Stalder et al., 1998, Klein et al., 2000, Corgne and Wood, 2002, Pertermann et al., 2004, Sun and Liang, 2011, Beyer et al., 2013). D values are not constants but may strongly depend on pressure, temperature and compositions of both mineral and melt. The most widely used model to predict mineral-melt partition coefficients is based on the lattice strain theory (Blundy and Wood, 1994). The model predicts partition coefficients by evaluating the mismatch in atomic radius and charge when a trace element is substituted for a major element in a regular crystal lattice. It has been applied and adjusted successfully in many studies (Wood and Blundy, 1997, van Westrenen et al., 2001, van Westrenen and Draper, 2007, Sun and Liang, 2011, Corgne et al., 2012, Davis et al., 2013, Sun and Liang, 2013).
By definition, the lattice strain model is focused on crystal structure and treats a potential influence of the corresponding melt indirectly. There is some debate about what control melt structure exerts on partitioning. Many recent studies address this issue with conclusions varying from insignificant influence to severe influence resulting in several orders of magnitude difference in D values (Kohn and Schofield, 1994, O’Neill and Eggins, 2002, Mysen, 2004, Prowatke and Klemme, 2005, Huang et al., 2006, Schmidt et al., 2006, Evans et al., 2008, Girnis et al., 2013). Results of such studies are even harder to interpret as, in contrast to crystals, the structure of most melts and the local environment of trace elements therein is usually not well constrained. Recently, Simon et al. (2013) addressed this issue and performed extended X-ray absorption fine structure (EXAFS) measurements of yttrium in silicate and aluminosilicate glasses with compositions taken from previous partitioning experiments (Prowatke and Klemme, 2005). They found a correlation between the preferential partitioning of Y into depolymerized melts and changes in Y coordination as well as Y–O distances in the glass. Molecular simulations provide another potentially powerful approach to investigate trace element partitioning. In particular, they offer a unique insight into structural properties of melts and crystals alike, while simultaneously providing access to many of their physicochemical and thermodynamic properties. Up till now most simulation studies have investigated element partitioning by statically calculating defect and solution energies in crystals and relating these properties to lattice strain theory (Purton et al., 1997, Purton et al., 2000, van Westrenen et al., 2000, Allan et al., 2003). On the other hand, studies explicitly focusing on melt properties are rare. Zhang and Yin (2012) performed two-phase coexistence first principles molecular dynamics (MD) simulations to study the partitioning of carbon and other light elements between metal and silicate melts. In an exploratory study, our research group used an approach based on classical MD to directly calculate trace element partitioning between a melt pair in the CaO-Al2O3-SiO2 system (Haigis et al., 2013). The MD simulations were combined with a thermodynamic integration scheme (Frenkel and Smit, 2002) to obtain the free energy difference of an element exchange reaction between the two melts, which directly relates to partitioning. Most recently, we applied a first principles MD approach combined with thermodynamic integration to explore the underlying mechanisms of Ni partitioning between metal and silicate melts (Künzel et al., 2017).
In this contribution we aim to further investigate the potential of this method by calculating trace element partition coefficients for several experimentally well known systems. The goals of this study are (a) to predict partition coefficients from first-principles and classical model systems where experimental data are available and (b) to elucidate the effect of melt structure on the predicted partitioning behavior. The first goal also serves as a benchmark for the capability of this method to make reliable predictions for systems where experimental data are lacking. We use molecular dynamics simulations in conjunction with thermodynamic integration (Marx and Hutter, 2012) to simulate partitioning of minor elements between several silicate melts. We explicitly focus on the effect of melt composition because the element exchange with a crystal usually involves coupled substitutions. Taking these into account would further complicate the thermodynamic integration approach we aim to test and needs to be discussed in future simulations involving this method. Therefore, our model systems are a gabbroic and granitic melt, similar in composition to coexisting melts investigated by Schmidt et al. (2006), as well as a Ti bearing silicate melt pair, resembling melts studied by Prowatke and Klemme (2005). Note that we by no means suggest that the two Ti-bearing melts are coexisting in nature. We rather evaluate the effect of melt composition on mineral-melt partitioning by assuming a virtually unchanged crystal as opposed to a changing corresponding melt (see original experimental study, Prowatke and Klemme, 2005, for details). We chose to investigate the trace elements Y, La and As for two reasons. Firstly, experimental data for these elements are readily available for the investigate melt systems and secondly, all three elements have the same formal charge (3+) but varying ionic radii (Y 1.04 Å, La 1.17 Å, As 0.58 Å, Shannon, 1976), making it possible to use the same substitution for Al3+ during thermodynamic integration (see below). The atomic interactions are either described by first-principles methods or classical force fields to compare their suitability for the thermodynamic integration approach. Also, the first-principles simulations may be used as reference points for end members of classical integrations (see Section 2.2).
First-principles MD is based on quantum mechanics, usually in the form of density functional theory (DFT, Hohenberg and Kohn, 1964, Kohn and Sham, 1965) and provides an essentially parameter-free description of the total energy of any given set of atoms. In many cases it has been shown to be highly accurate and predictive, also in investigating the structure and properties of melts over a range of pressures and temperatures (Alfè et al., 2002, de Koker et al., 2013, Muñoz Ramo and Stixrude, 2014). However, first-principles methods are limited by their comparatively high computational cost, allowing simulations of typically up to about a hundred atoms for MD simulations of complex silicate melts. Classical simulations, on the other hand, rely on force fields to describe atomic interactions. This makes the simulation of significantly larger systems and longer simulation times computationally affordable and thus enables the study of processes such as diffusion or the structure of melts and fluids beyond nearest neighbor distances. However, classical force fields can only be as accurate as the training set used for its parametrization. Here, we use an advanced ionic interaction potential that has been successfully applied to simulate oxides and silicate crystals, melts and glasses (e.g. Jahn and Madden, 2007, Adjaoud et al., 2011, Finkelstein et al., 2014).
2. Ni partitioning between metal and silicate melts: An exploratory ab initio molecular dynamics simulation study
1. Introduction
Element partitioning is a fundamental geochemical process that has occurred throughout the whole history and in all parts of the Earth. Scientific interest to understand element partitioning, e.g., for constraining the conditions during Earth accretion and its further evolution has been prevalent for decades (Goldschmidt, 1937, Ringwood, 1966). Besides measuring element concentrations in natural samples, experimental approaches have been used widely to determine element partition coefficients. In particular, many of these experimental studies were focussed on systems related to partitioning processes during Earth and planetary formation, e.g. Boujibar et al., 2014, Chi et al., 2014, Mann et al., 2012, O’Neill et al., 1995, Righter and Ghiorso, 2011, Rubie et al., 2003, Siebert et al., 2012 (Siebert et al., 2012, Siebert et al., 2011), Thibault and Walter, 1995, Wade and Wood, 2005, Wade et al., 2012 , and Rai and van Westrenen (2014). One of the biggest challenges in recent years has been to access the extreme conditions of pressure and temperature during accretion and deep inside planetary bodies and exciting progress has been made in this respect. For instance, metal–silicate melt partition coefficients have been determined up to conditions of the base of a presumed terrestrial magma ocean, from which unique information about Earth's core formation could be obtained (e.g. Bouhifd and Jephcoat, 2011, Siebert et al., 2012 ).
Knowledge of the structural properties of the phases between which elements are distributed can be instrumental in obtaining a deeper understanding of the processes involved. While most of the relevant crystal structures are relatively well known, the structure of melts (or fluids) and especially the local atomic coordination environment of the fractionating elements are not well constrained in many cases. Yet, structural information about melts and fluids on the atomic scale is even more challenging to derive from experiments, especially for minor and trace elements. Recently, silicate melt structures have been studied up to extreme conditions, including those of Fe2SiO4 (Sanloup et al., 2013) or aluminosilicate melts (Drewitt et al., 2015).
In recent years, computational performance and methods have expanded greatly, so that it is now possible to use ab initio methods to model element exchange reactions and to receive structural information at the same time. These methods are based on quantum-mechanics, often in the framework of density functional theory (DFT). In combination with the molecular dynamics method (Marx and Hutter, 2012), melt structures as well as their physical and thermodynamic properties have been predicted in a wide range of pressures and temperatures (e.g. Alfè et al., 2002; de Koker et al., 2008, Muñoz Ramo and Stixrude, 2014; de Koker et al., 2013). Many thermodynamic properties require the evaluation of the Helmholtz or the Gibbs free energy, which is not directly accessible by standard molecular dynamics methods. One way of calculating free energy differences is to combine molecular dynamics simulations with the thermodynamic integration method (Frenkel and Smit, 2002). For instance, the free energy difference of element exchange between silicate melts with different degrees of polymerization was examined recently using the thermodynamic integration method of alchemical transmutation in conjunction with classical molecular dynamics simulations (Haigis et al., 2013).
The goal of this study is to explore the possibility to predict element partition coefficients between metal and silicate melts using an ab initio approach. We will combine ab initio molecular dynamics simulations with the thermodynamic integration method to model the metal–silicate partitioning of Ni at ambient pressure, where both structural and partitioning data is available from experiments (e.g. Thibault and Walter, 1995; Farges and Brown, 1996, Schenk et al., 2002, Drewitt et al., 2013; Sanloup et al., 2013). Ni is assumed to be the second most abundant element in the Earth's core and Ni partitioning between metal and silicate melts has been studied extensively (Thibault and Walter, 1995; Bouhifd and Jephcoat, 2011; Siebert et al., 2012), e.g., to constrain the depth of the magma ocean. Thus, if eventually accurate enough, this computational approach has the potential to provide unprecedented insight into partitioning behavior also at extreme conditions and for other important elements such as Co and Cr.
3. Carbon and other light element contents in the Earth’s core based on first-principles molecular dynamics
The high solubility of C in liquid iron (1, 2) and the existence of graphite and cohenite (Fe,Ni)3C in iron meteorites (3) suggest that C could account for the density deficit of the Earth’s core (4–6). Carbon content also greatly influences siderophile and chalcophile element partitioning between metal and silicate, and hence their distribution in the Earth’s core and mantle (7, 8). In particular, its influence on Pb partitioning would affect the inference of the age of the Earth (9). Therefore, it is critically important to constrain the C content of the Earth’s core.
There are several ways to infer the C content of the core. The first method is a solubility and phase diagram study in the Fe–C system, which generally gives a high C content in the Earth’s core, on the order of 2–4 wt% (1). The implicit assumption is that C is supersaturated during the Earth’s accretion process, and so this study will generally give the maximum C content of the core. The second method involves density and bulk modulus measurements of solid Fe3C and Fe7C3 under the inner core pressures and comparison with the seismic data of the inner core. This kind of study results in a C content of 1–1.5 wt% (10, 11) and assumes that C is the only light element in the inner core, thus giving only an upper limit for the C content in the inner core. Application of the density and bulk modulus measurements of solid and liquid Fe3C to the outer core exclude C as a major alloying element in the Earth’s outer core (12, 13). The third method is based on carbonaceous chondrite composition, the element volatility curve of the bulk Earth, and the mass balance between primitive mantle and core. This method gives the smallest C content of the core: 0.2 wt% (14).
Thus, the C content of the core has proven to be difficult to constrain, and current estimates differ by a factor of ∼20. Some limits on the geochemical role of C in the core-forming process can be provided by P, which becomes lithophile during C saturation (7). Because P is depleted rather than enriched in the mantle, core formation probably did not occur at C saturation (7). Thus, we contend that the starting bulk composition of the Earth may have contained a lower C content than some of the experiments outlined above indicated.
Carbon is thought to have entered the core through core–mantle segregation processes in the early history of the Earth during the magma ocean stage. Inferring the amount of C involved thus requires determining the partition coefficient of C between liquid iron and silicate melt experimentally. However, as yet no such data exist, as pointed out by Dasgupta and Walker (15). Many partitioning experiments for siderophile elements use C as a sample encapsulating material, thus greatly affecting the behavior of these elements. In the present study, we use the two-phase first-principles molecular dynamics (FPMD) method (16) to determine the partition coefficient of C and a few other light elements between liquid iron and silicate melt. This partition coefficient can then be combined with the composition of the bulk Earth or the primitive mantle (i.e., bulk silicate Earth) to infer the C content of the Earth’s core.
4. Partitioning of Si and O between liquid iron and silicate melt: A two-phase ab-initio molecular dynamics study
1. Introduction
[2] The Earth's early accretion events and processes, in particular, the giant impact and the magma ocean process have a great influence on the evolution and the current state of the Earth [Stevenson, 2008]. For example, light elements of the core, with their identities and contents largely determined by the magma ocean process, may help sustain the geodynamo through compositional convection [Braginsky, 1963], and influence inner core growth through melting point depression [Jeanloz, 1990].
[3] To understand the magma ocean process and its effect on the mantle and the core composition, many element partitioning experiments have been performed to derive the temperature, the pressure, and the oxygen fugacity at which the magma ocean process occurs [e.g., Li and Agee, 1996; O'Neill et al., 1998; Walter et al., 2000; Righter, 2003; Wade and Wood, 2005; Corgne et al., 2008; Wood et al., 2008, and references therein]. The inferred temperatures and pressures range from 2000 to 4200 K and 25 to 40 GPa, reflecting both experimental difficulties and complexity of the magma ocean process [Rubie et al., 2007; Wood, 2008].
[4] Ab-initio methods have been applied to study element partitioning between the inner and the outer core by calculating chemical potentials of candidate elements in solid and liquid iron [Alfè et al., 2000, 2002]. The two-phase method, which makes two phases in direct contact in a molecular dynamics (MD) simulation cell, is employed in ab-initio studies of the melting temperatures of iron and MgO [Alfè, 2005, 2009] and the solubility of argon in liquid iron [Ostanin et al., 2006].
[5] In the present study, we extend the ab-initio two-phase method so that compositions of the two segregated phases can be calculated. Starting from a random distribution of elements (O, Mg, Fe, and Si) in the simulation cell, with the MD progressing, the system segregates into two phases, one is liquid iron and the other is silicate melt. Compositions of the two phases are calculated by using concepts and methods in computational geometry and compared with experimental data. The method is then applied to infer the amounts of Si and O that could be added by the deep magma ocean process to the Earth's core.
5. Magnesium partitioning between silicate melt and liquid iron using first-principles molecular dynamics: Implications for the early thermal history of the Earth's core
1. Introduction
Fluid convection in the Earth's core generates magnetic field through the geodynamo process. The energy for driving the geodynamo may be derived from thermal and compositional buoyancy, Earth's tides and precession, and/or radioactivity (e.g. Nimmo, 2015, Andrault et al., 2016). Recent first-principles calculations find that iron and iron alloys have a very high thermal conductivity (de Koker et al., 2012; Pozzo et al., 2012), which are further supported by high-pressure and high-temperature experiments in diamond-anvil cells (Ohta et al., 2016, Sata et al., 2010, Seagle et al., 2013). High values of thermal conductivity, when combined with the very early existence of Earth's palaeomagnetic field (Tarduno et al., 2010, Tarduno et al., 2015) imply that the Earth's core temperature has to start very high, at around 7800 K at the time of 4.2 Ga ago (Labrosse, 2015), and then needs to drop gradually. Even if radioactive element such as potassium exists in Earth's core at the amount far beyond its experimentally allowed value (Bouhifd et al., 2007, Corgne et al., 2008, Murthy et al., 2003, Watanabe et al., 2014), the temperature at the core-mantle boundary is still required to be as high as 5720 K (Labrosse, 2015). However, due to the expected fast cooling of magma ocean, the initial high temperature of the core must have cooled rapidly (Andrault et al., 2016, Monteux et al., 2016), making the thermal energy insufficient to drive the geodynamo.
Compositional buoyancy created by the release of light elements at the inner core boundary (ICB) is thought to be another principal energy source for the geodynamo (Stacey, 1992). This mechanism would require an early crystallization of the inner core. However, Recent palaeomagnetic study reveals that the Earth's inner core may have formed very late (Biggin et al., 2015, Smirnov et al., 2016), making the compositional buoyancy due to inner core crystallization an unlikely energy source in the very early history of the Earth.
To reconcile these difficulties, precipitation of magnesium is proposed recently as a dominant energy source for driving the geodynamo in the early history of the Earth (O'Rourke and Stevenson, 2016, Badro et al., 2016). Magnesium normally has a very low solubility in liquid iron. However, Mg could be dissolved into the core effectively at a much higher temperature achieved through violent accretion processes, such as multiple giant Impacts, including the last one that forms the Earth's Moon. Consequently, magnesium content of the core can be at much higher level than previously anticipated. Because of its very low solubility, magnesium would precipitate very early during the subsequent cooling and provide the needed buoyancy force to drive the geodynamo (O'Rourke and Stevenson, 2016; Badro et al., 2016, Badro et al., 2018; O'Rourke et al., 2017).
To study how much magnesium can be dissolved into Earth's core during the Giant Impact processes and at which temperature magnesium will be precipitated from the core, knowledge of magnesium partitioning between liquid iron and silicate melt is essential. However, equilibrium constant of magnesium partitioning is still under debate. Some previous studies find or simply assume (O'Rourke and Stevenson, 2016; Badro et al., 2016, Badro et al., 2018) that the equilibrium constant is dependent on temperature, whereas another experimental study (Du et al., 2017) argue that it is independent of temperature. In the latter case, the magnesium content of the core would be a constant as the core temperature decreases, and the Mg exsolution would have difficulty to drive the geodynamo.
In the present study, we first apply first-principles molecular dynamics simulations to calculate the equilibrium constant of magnesium partitioning between liquid iron and silicate melt and to verify its dependence on temperature. Then, we combine the equilibrium constant with the thermal evolution of the earth's core and palaeomagnetic evidences to discuss whether oxide or silicate can precipitate and provide enough entropy to sustain the geodynamo before the inner core nucleation. Several circumstances, mainly concerned with the initial composition of the core, are considered. These include the Moon-forming impactor's size, which is the major factor in controlling the initial magnesium content of the core, and the so-called “Grand Tack” and other accretion scenarios which mainly influence the initial silicon and oxygen contents of the core. Finally, we present a brief history of the early evolution of the Earth's magnetic field based on the thermal evolution of the core and paleomagnetic observation
6. Molecular dynamics simulations of Y in silicate melts and implications for trace element partitioning
1. Introduction
In the presence of two or more coexisting phases in thermodynamic equilibrium, a minor or trace element will, in general, not be distributed equally among these phases, but will be incorporated preferentially into some chemical environments at the expense of others. The resulting distribution of an element i between two phases α and β is quantified by the Nernst partition coefficient Diα/β = ciα/ciβ, with ciα denoting the concentration (mass fraction) of element i in phase α. The molar partition coefficient Di ∗α/β = xiα/xiβ is defined in terms of mole fractions x instead of concentrations and can easily be converted to Diα/β (see Beattie et al. (1993), for terminology). The partition coefficient of element i depends, in general, on temperature, pressure, chemical composition and structure of the involved phases. Conversely, if this dependence is known, either from a compilation of experimental data or from a suitable theory, then the distribution of trace elements in, e.g., rock samples can provide information about the petrogenetic history, and hence constitute a valuable tool for petrologists and geochemists (e.g. Shaw (2006)).
Trace element partitioning between coexisting crystal and melt phases can be understood, at least partially, in terms of the local environment of the incorporated cation in the crystal: if the trace element fits available crystal sites well by size and charge, then it is enriched in the crystal, otherwise it partitions into the melt. It has long been known that crystal-melt partition coefficients of a series of isovalent cations, plotted logarithmically as a function of ionic radius, form near-parabolic patterns which peak at an ideal radius (Onuma et al., 1968). These patterns have been described quantitatively by Blundy and Wood (1994) by means of the lattice strain model, based on work by Brice (1975). It translates the strain in the crystal lattice, induced by the mismatch of the incorporated cation, to a free energy penalty for cation incorporation, which in turn governs the partitioning.
The lattice strain model, when suitably parameterized, successfully describes observed crystal-melt partitioning behavior in terms of crystal chemistry alone, without explicitly taking into account melt properties. If the latter are important, their influence on partitioning is hidden in the adjustable model parameters and cannot be predicted nor explained by the original model (although Wood and Blundy (1997) extended the model by taking into account the Mg/(Mg + Fe) ratio in the melt). However, there is broad evidence that melt composition can indeed have a strong effect on trace element partitioning. Prowatke and Klemme (2005), in a series of experiments, measured partition coefficients of several trace elements between titanite (CaTiSiO5) and a range of coexisting silicate melts of different compositions. Although the crystal chemistry was virtually constant in all experiments, partition coefficients varied by two orders of magnitude for several rare-earth elements (REE) and Th. Moreover, the partition coefficients were found to depend systematically on melt polymerization, quantified by the molar ratio of Al2O3/(Na2O + K2O + CaO). In particular, all the REE probed in the study showed an increasing tendency to partition into the melt the more the melt was depolymerized.
More evidence for the influence of melt composition (and thus melt structure) on element partitioning comes from experiments with immiscible silicate melts (Watson, 1976, Ryerson and Hess, 1978, Schmidt et al., 2006). In the latter study, partition coefficients between coexisting gabbroic (highly depolymerized) and granitic (highly polymerized) melts were determined for a large set of elements. A strong preference of the REE for the depolymerized melt was found, with DREEgabbro/granite ≃ 10. According to the authors’ interpretation, the abundance of non-bridging oxygen in the depolymerized melt facilitates the formation of the preferred (i.e. energetically favorable) coordination polyhedra of REE and thus favors the observed distribution. In a similar study, Veksler et al. (2006) found the same partitioning trend between immiscible pairs of Fe-rich and Si-rich melt. Mysen (2004) suggested to understand the influence of melt composition on element partitioning in terms of Qn (0 ≤ n ≤ 4) species whose abundance and proportions vary with bulk composition, thus offering varying amounts of energetically favorable “sites” for trace element incorporation.
It is the aim of this computational study to elucidate the mechanisms by which melt composition influences the distribution of trace elements between crystal and melt or between two melts. We chose Y as an exemplary REE whose partition behavior was shown to depend strongly on melt properties (Prowatke and Klemme, 2005, Schmidt et al., 2006). From a computational point of view, Y is a more convenient element than the (chemically similar) lanthanides whose strongly correlated 4f electrons pose notorious problems for theoretical descriptions. Our approach is based on molecular dynamics, a method which provides simultaneous access to the atomic structure and dynamics and to the thermodynamic variables of a system. Taking four different model silicate melts, we first investigated how melt composition influences the local coordination environment of Y ions in the melt. We then translated these structural changes to differences in free energy, which in turn determine the partitioning behavior of Y.
7. Coordination of Zr4+/Hf4+/Nb5+/Ta5+ in silicate melts: insight from first principles molecular dynamics simulations
1. Introduction
As strategic metals, Zr, Hf, Nb and Ta have been attracting more and more attentions because of their indispensable role in emerging industries (U.S. Department of the Interior and U.S. Geological Survey, 2017; European Union, 2018). It has been accepted that magmatic processes especially carbonatite and highly fractionated granitic magmas are responsible for their ore formation (Pollard, 1995; Linnen et al., 2014). The granites related to high field strength elements (HFSE) deposits are usually rich in fluorine (Linnen and Cuney, 2005; Černý and Ercit, 2005) and it is proposed that fluorine may have played an important role in their transport and enrichment (e.g., Li et al., 2015; Demartis et al., 2014; Charoy and Raimbault, 1994).
The experiments by Keppler (1993) showed that the solubilities of zircon (ZrSiO4), hafnon (HfSiO4), manganocolumbite (MnNb2O6) and manganotantalite (MnTa2O6) increased with the fluorine content in water-saturated haplogranitic melts. He proposed two possible mechanisms for this phenomenon: (1) metal ions formed complexes with F−, and (2) F− increased the number of non-bridging oxygens (NBO) in the melts. However, the later experiments of Fiege et al. (2011) and Aseri et al. (2015) in water-saturated natural granitic and haplogranitic melts with F content up to 11 wt% indicated that fluorine increased the solubilities of zircon and hafnon, but did not influence those of manganocolumbite and manganotantalite. They attributed their different results of manganocolumbite and manganotantalite to the possible nonequilibrium in the experiment of Keppler (1993), and proposed a new mechanism that Zr4+/Hf4+ could form complexes with F− while Nb5+/Ta5+ could not. The key to verify this conjecture is to determine the local coordination structures of Zr4+, Hf4+, Nb5+ and Ta5+ in silicate melts.
XAFS (X-ray absorption fine structure) has been used to probe the coordinations of Zr4+, Nb5+ and Ta5+ in silicate melts and quenched glasses (Farges et al., 1991; Farges, 1996; Piilonen et al., 2006; Wilke et al., 2012; Louvel et al., 2013; Mayanovic et al., 2013). These studies indicated that the metal cations are all mainly 6-fold coordinated with the average Zr4+—O2− bond length of about 2.1 Å and the average Nb5+/Ta5+—O2− bond lengths of about 2.0 Å. In contrast, Farges et al. (1991) proposed an 8-fold coordination for Zr4+—O2− in highly polymerized glasses. XAFS did not detect significant spectroscopic difference in the local structures of Zr4+, as well as Nb5+, between F-free and F-bearing silicate melts or quenched glasses (Farges et al., 1991; Farges, 1996; Louvel et al., 2013; Piilonen et al., 2006). However, these results did not rule out the Zr4+—F− complexation, e.g. the ab initio spectra calculations by Louvel et al. (2013) indicated that the XAFS analysis could not distinguish Zr4+—F− complexation from Zr4+—O2− because of the similar spectra. Therefore, the coordination chemistry of Zr4+, Hf4+, Nb5+ and Ta5+ in F-bearing silicate melts is still inconclusive, which impedes the understanding of the effect of fluorine on the ore-forming processes of these elements.
As an implement of experimental techniques, first principles molecular dynamics (FPMD) has been widely applied to study the speciation of metal ions in hydrothermal fluids (Brugger et al., 2016). Jahn et al. (2015) investigated the complexation of Zr4+ and Hf4+ in supercritical aqueous fluids by using FPMD. Numerous FPMD simulations of various melts including oxides (e.g., Karki et al., 2007; Verma et al., 2011), silicates (e.g., Karki and Stixrude, 2010a; Ni and de Koker, 2011; Sun et al., 2019) and carbonates (e.g., Vuilleumier et al., 2014; Zhang and Liu, 2015) have been reported. These studies focused on the physical properties of melts, but less attention has been paid to the geochemical behavior of ore-forming metals in melts (e.g., Wagner et al., 2017, Wagner et al., 2017b). To the best of our knowledge, the coordination chemistry of Zr, Hf, Nb and Ta in silicate melts has not been investigated by using FPMD.
In this study, we performed FPMD simulations of Zr4+, Hf4+, Nb5+ and Ta5+ in anhydrous and hydrous F-free/F-bearing albite melts to investigate the chemical interaction between these metal ions and F−. We derived the coordination structures of these metal ions, and found that they did not form complexes with F−. F− bound with Si and Al, which generated NBO in silicate melts. Based on the results, we proposed a mechanism to explain the different effects of fluorine on the solubilities of zircon/hafnon and manganocolumbite/manganotantalite in silicate melts (Fiege et al., 2011; Aseri et al., 2015), and discussed the geological implications.
8. Ab Initio Prediction of Potassium Partitioning Into Earth's Core
1 Introduction
The generation of Earth's magnetic field, convection in the mantle and outer core, and the age of the inner core are intimately governed by energy transferred from the Earth core. The energy sources are mainly attributed to the gravitational energy accumulated when the proto-Earth was accreted, the latent heat released at the inner core crystallization, and the gravitational energy released from the light element exclusion that is associated with the inner core solidification (Buffett et al., 1996; Gibbins et al., 2004). Another possible contribution could be from a radioactive decay of long-lived radionuclides, including potassium, uranium, and thorium that are potentially present in the core (Labrosse et al., 2001; Nimmo et al., 2004).
Potassium, uranium, and thorium are lithophile elements and expected to be distributed in the silicate Earth. The concentrations of uranium and thorium in the bulk silicate Earth (BSE) are comparable with CI-chondrites, allowing only limited amounts to be incorporated into Earth's core (McDonough & Sun, 1995). In contrast to uranium and thorium, potassium is possibly favorable for being present in the core due to its chemical features. The concentration of potassium in the BSE is ~240 ppm, which is strongly depleted compared with CI-chondrites (McDonough & Sun, 1995). The substantial evaporation loss of potassium during the early stage of the Earth evolution was, however, denied by the evidence showing similar potassium isotopic ratios in samples from Earth, Moon, and chondrites (Humayun & Clayton, 1995). The existence of an incompatible-element-enriched layer hidden in deep mantle is also ruled out by geodynamics studies (Campbell, 2001; Helffrich & Wood, 2001). Another possibility is the incorporation of potassium into the core, which can compensate its depletion in the BSE compared to CI-chondrite.
The amount of potassium that can be incorporated into the core is determined by its partitioning between silicate and metallic systems. A number of laboratory experiments have therefore been performed to measure the metal/silicate partitioning coefficient of potassium (e. g., Blanchard et al., 2017; Bouhifd et al., 2007; Chabot & Drake, 1999; Corgne et al., 2007; Gessmann & Wood, 2002; Murthy et al., 2003; Watanabe et al., 2014). However, the results scatter significantly (by four orders of magnitude) potentially due to wide variations of the potassium partitioning depending on temperature, pressure, composition, and/or oxygen fugacity. In this study, free energy calculations have been conducted to clarify the metal/silicate partitioning behavior of potassium using the ab initio molecular dynamics simulations combined with thermodynamics integration technique. Potential controlling factors including temperature, pressure, and metallic composition with different types of light elements are investigated. Electronic and atomic structures are then analyzed to elucidate the mechanisms of potassium partitioning.
9. Ab initio chemical potentials of solid and liquid solutions and the chemistry of the Earth’s core
I. INTRODUCTION
We present here a set of techniques that allow the ab initio calculation of chemical potentials in solid and liquid solutions, and hence the ab initio treatment of chemical equilibrium between solid and liquid phases. There are many areas of chemical physics where such techniques might be important, but we believe they have a particular role to play in studying the partitioning of impurities between different phases under extreme conditions, where experiments are difficult or impossible. As an illustration of the power of the techniques, we will describe how we have applied them to study chemical equilibrium between the solid and liquid parts of the Earth’s core.
The techniques to be presented form a natural sequel to recent developments in the ab initio thermodynamics of condensed matter based on electronic density-functional theory (DFT). For many years, DFT has been used to calculate the phonon spectra of perfect crystals, and it is only a short step from that to the calculation of free energies and other thermodynamic quantities in the harmonic approximation. There have already been several reports of DFT calculations of high-temperature crystal thermodynamics, including solid–solid phase equilibria using this approach. For liquids, ab initio thermodynamics first became possible with the Car–Parrinello technique of DFT molecular dynamics simulation, which immediately gave a way to calculate such quantities as pressure, internal energy and temperature of a liquid in thermal equilibrium. The first DFT treatment of solid–liquid equilibrium was achieved by Sugino and Car, who used thermodynamic integration to compute the free energies of solid and liquid Si, and hence the melting properties of the material. Closely related is the work of de Wijs et al. on the melting of Al. We have recently reported DFT calculations of the free energies and melting curves of Fe (Refs. 9, 17, and 18) and Al (Ref. 19) over a wide range of pressures; Jesson and Madden have recently presented ab initio calculations of the zero-pressure melting properties of Al using their ‘‘orbital free’’ approach. The work of Smargiassi and Car and Smargiassi and Madden on the free energy of formation of defects in crystals is also relevant to the ideas to be presented here.
Thermodynamic integration has been the key to calculating the ab initio free energy of liquids and anharmonic solids, and hence to the treatment of solid–liquid equilibrium. It provides a means of computing the difference of free energy between the ab initio system and a reference model system whose free energy is known. We will show that it is also the key to calculating ab initio chemical potentials of liquids and anharmonic solids, but here it is used in a rather different way. The chemical potential of a species is the free energy change when an atom of that species is added to the system. The difference of chemical potentials of two species is therefore the free energy change when an atom of one species is replaced by an atom of the other, or equivalently when one atom is transmuted into the other. The role of thermodynamic integration here is to provide a way of calculating the free energy change associated with such transmutations, and we shall show how this can be accomplished in practical DFT simulations. This general approach is closely related to ideas that have been used for a long time in classical simulation. A recent example of classical thermodynamic integration with molecular transmutation to calculate solvation free energies in aqueous solution can be found in Ref. 23, which gives references to earlier literature.
Although the techniques we shall present are fairly general, we do impose two restrictions at present: First, the theoretical framework is developed for the case of a two-component mixture; second, one of the components is present at low, but not very low, concentration, in a sense to be clarified in Sec. II. The situation envisaged therefore consists of fairly dilute solid and liquid solutions in coexistence. There are vast numbers of problems both in the chemical industry and in the natural world that depend on the partitioning of chemical components between coexisting phases, and the ability to calculate chemical potentials ab initio should make it possible to address some of these problems in a new way. Our original incentive for developing these techniques was the desire to understand better the chemistry of the Earth’s core, and this is a good example of a problem where ab initio calculations can supply information that is difficult to obtain experimentally because of the extreme conditions of temperature (T; 6000 K!) and pressure (p; 330 GPa!). The core is composed mainly of Fe, and comprises an outer liquid part and an inner solid part. The density of the outer core is; 6% too low to be pure Fe, and cosmochemical and geochemical arguments show that the main light impurities are probably S, O, and Si. The inner core has grown over geological time by crystallization from the outer core, and the partitioning of impurities between liquid and solid is crucial for understanding the evolution and contemporary dynamics of the core. The size of the density discontinuity (~4.5%! Refs. 29 and 30!) at the inner-core/outer-core boundary can only be interpreted if one understands this partitioning, and also provides a constraint on possible chemical compositions. We shall show how our ab initio techniques for calculating chemical potentials shed completely new light on this problem. Brief reports of these calculations have already appeared.
In developing the theoretical basis of our techniques, we define our technical aims in Sec. II by summarizing the standard thermodynamic relations describing phase equilibrium. The difference of chemical potentials of solute and solvent atoms, and the free energy change associated with the transmutation of solvent into solute are discussed in Sec. III. In Sec. IV, we then develop the ab initio techniques themselves. We shall explain how thermodynamic integration can be used to perform the solvent–solute transmutation so as to obtain the difference of chemical potentials in the liquid; we also describe the techniques for calculating chemical potentials in solid solutions, both in the harmonic approximation and using thermodynamic integration to handle anharmonicity. Section V presents our results for the case of S, O, and Si dissolved in solid and liquid Fe under Earth’s core conditions, and summarizes the implications of the results for the partitioning of these impurities between the inner and outer core and the chemical composition of the core. Discussion and conclusions are given in Sec. VI.
我们在这里提出了一套技术，允许从头计算固体和液体溶液中的化学势，从而从头处理固液相之间的化学平衡。这些技术在许多化学物理领域可能都很重要，但我们认为它们在研究在极端条件下不同相之间杂质的分配的情况下尤为重要，而在这些条件下实验是困难或不可能的。作为这些技术强大性能的说明，我们将描述我们如何应用它们来研究地球核心的固体和液体部分之间的化学平衡。
这里将介绍的技术是从基于电子密度泛函理论（DFT）的凝聚态物质的从头热力学的最新发展自然而然的续篇。多年来，DFT已被用于计算完美晶体的声子谱，从而很容易从中计算在谐振近似中的自由能和其他热力学量。已经有几份关于使用这种方法计算高温晶体热力学，包括使用DFT进行的固-固相平衡的DFT计算的报告。对于液体，从头热力学首次成为可能是通过DFT分子动力学模拟的Car–Parrinello技术，它立即提供了计算液体处于热平衡状态时的压力、内能和温度等量的方法。Sugino和Car首次实现了对固-液平衡的DFT处理，他们使用热力学积分计算了固体和液体硅的自由能，从而计算出了该材料的熔化性质。de Wijs等人对铝熔化的工作与此密切相关。我们最近报告了在广泛压力范围内对Fe（参考文献9、17和18）和Al（参考文献19）的自由能和熔化曲线的DFT计算；Jesson和Madden最近使用他们的‘‘轨道自由’’方法对Al的零压熔化性质进行了从头计算。Smargiassi和Car以及Smargiassi和Madden关于晶体中缺陷的形成自由能的工作也与我们这里要介绍的思想相关。
热力学积分已成为计算从头自由能的液体和非谐固体的关键，因此也是处理固-液平衡的关键。它提供了计算从头系统与其自由能已知的参考模型系统之间的自由能差异的方法。我们将展示，它还是计算液体和非谐固体的从头化学势的关键，但在这里使用的方式略有不同。一种物种的化学势是当该物种的原子被添加到系统中时的自由能变化。两种物种的化学势差因此是当一种物种的原子被另一种物种的原子替代时的自由能变化，或者等效地，当一个原子被转变为另一个原子时的自由能变化。热力学积分在这里的作用是提供一种计算与这种转变相关的自由能变化的方法，我们将展示如何在实际的DFT模拟中实现这一点。这一通用方法与长期以来在经典模拟中使用的思想密切相关。可以在参考文献23中找到使用分子转变进行经典热力学积分计算水溶液中溶解自由能的最近例子，该文献提供了早期文献的参考文献。
虽然我们将要介绍的技术相当通用，但目前我们确实对其施加了两个限制：首先，理论框架是为两组分混合物的情况而开发的；其次，其中一种成分以低浓度存在，但并非非常低，这一点将在第二节中进行澄清。因此，所设想的情境包括相对稀薄的固体和液体溶液共存。在化工行业和自然界中，有大量问题依赖于化学成分在共存相之间的分配，而从头计算化学势的能力应该能够以一种新的方式解决其中一些问题。我们发展这些技术的最初动机是更好地理解地球核心的化学性质，这是一个很好的例子，说明了在实验条件因温度（T; 6000 K！）和压力（p; 330 GPa！）极端而难以获取实验数据的情况下，从头计算能够提供信息的问题。地核主要由Fe组成，包括一个外部液体部分和一个内部固体部分。外核的密度低约6%，不能是纯Fe，宇宙化学和地球化学的论据显示，主要的轻杂质可能是S、O和Si。内核已经通过从外核结晶而在地质时间内增长，而在液体和固体之间的杂质分配对于理解核的演变和当代动力学至关重要。只有在理解了这种分配的情况下，才能解释内核/外核边界处密度不连续性的大小（约4.5％！参考文献29和30！），并且还为可能的化学组成提供了一个约束。我们将展示，我们用于计算化学势的从头技术如何为这个问题带来全新的视角。已经出现了这些计算的简要报告。
在发展我们技术的理论基础时，我们在第二节中通过总结描述相平衡的标准热力学关系来定义我们的技术目标。在第三节中，讨论了溶质和溶剂原子的化学势差异以及将溶剂转变为溶质所伴随的自由能变化。在第四节中，我们进一步发展了从头技术本身。我们将解释如何使用热力学积分执行溶剂-溶质转变，以获得液体中化学势的差异；我们还描述了在谐振近似和使用热力学积分处理非谐振动的情况下计算固体溶液中化学势的技术。第五节呈现了在地球核心条件下S、O和Si溶解在Fe固体和液体中的情况的结果，并总结了这些结果对这些杂质在内外核之间分配以及核心化学组成的影响。讨论和结论在第六节中给出。