First Principles Calculation of the Stability of Iron Bearing Carbonates at High Pressure Conditions

Abstract

Carbonate minerals such as ferromagnesite (Mg,Fe)CO${}_3$ are suggested to be a possible major deep-carbon host in the lower mantle, because ferromagnesite is possibly stabilized by Fe spin crossover under pressure. However, the behavior of Fe-bearing carbonates under lower mantle pressure conditions has not been suitably examined thus far. Thus, in this study, we investigate the high-pressure stability of ferromagnesite and possible high-pressure structures with the chemical composition of (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ via first principles calculation using internally consistent local density approximation with Hubbard parameter (LDA+U) method, which can more accurately account for the electronic state of Fe than the LDA and generalized gradient approximation (GGA) approaches. The enthalpy values obtained via our calculations suggest that (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ undergoes phase transition from the R$\overline{3}$c structure (high spin) to the P$\overline{1}$ (high spin) at 50 GPa, and to C2/m (high-spin) structure above 80 GPa, under static 0 K conditions. Therefore, no spin transitions in these carbonate minerals is expected under the lower mantle pressure conditions.

1. Introduction

Carbon is one of the important volatile elements in the Earth. Accordingly, it is known to significantly affect the climate as well as dynamics processes on the planet, such as volcanic activity. Carbon is primarily transported into deep Earth in the form of carbonate minerals, such as magnesite (MgCO${}_3$), calcite and aragonite (CaCO${}_3$), dolomite (CaMg(CO${}_3$)${}_2$), and siderite (FeCO${}_3$) [1,2]. Among them, magnesite has attracted special attention from researchers because of its high stability under high pressure and temperature conditions in the Earth’s interior [3,4,5,6,7,8]. However, in general, the circulation mechanism of carbon in the deep interior remains largely unknown [9,10].

High pressure and temperature experiments conducted in some previous studies using a diamond-anvil cell (DAC) found that magnesite did not undergo phase transformation at 80 GPa between temperatures of 2000–2500 K [3,4,5]. In contrast, through an in situ X-ray diffraction experiment conducted in a later work, it was shown that magnesite transformed into an unknown form at ∼115 GPa and 2100–2200 K [7]. Recently, in a study wherein in situ X-ray diffraction measurement using a double-sided laser-heated DAC was conducted, a phase transformation from magnesite to magnesite-II (space group C2/m) was reported at 80 GPa [8], which had previously been predicted based on ab initio calculations [11,12].

The structural evolution of magnesite under ultrahigh pressure conditions has been explored theoretically [11,12]. Oganov et al. [11] reported that the magnesite structure transformed to space group C2/m with a lower enthalpy than that of magnesite above 82.4 GPa; this transformation was predicted based on calculations performed using ab initio evolutionary crystal structure searching algorithm. In addition, they also reported that the C2/m crystal structure transformed to one with space group P2${}_1$ above 138 GPa. Pickard and Needs [12] also conducted ab initio random structure searching of the MgCO${}_3$ system and found that its new crystal structure (space group P$\overline{1}$) has a lower enthalpy than that of the C2/m phase in the range of 80–101 GPa.

The presence of Fe in carbonate minerals has been shown to modify their phase stability under high pressure and temperature conditions [13,14,15]. Ferromagnesite (Mg,Fe)CO${}_3$ is a solid solution of magnesite and siderite. Several different results on the high pressure behaviors of ferromagnesite have been reported thus far. In particular, siderite-rich ferromagnesite (Fe${}_{0.75}$Mg${}_{0.25}$)CO${}_3$ has been reported to transform into a new Fe${}^{3+}$-bearing high pressure carbonate, Fe${}_3$O${}_4$ magnetite, and nano-diamonds above 80 GPa and 1850–2300 K [2]. In contrast, it was found that (Fe${}_{0.65}$Mg${}_{0.35}$)CO${}_3$ transforms into an orthorhombic phase (phase II with space group Pmm2) above ∼50 GPa and 1400 K; this was determined via a synchrotron X-ray diffraction experiment using a laser-heated DAC [14]. It should be noted that these previous experiments were strongly affected by the oxygen fugacity in the system at high pressure–temperature conditions [2,14,16].

Moreover, it has been suggested that spin crossover of Fe in high-pressure carbonate phases stabilizes them at relatively lower pressure conditions [14]. In particular, high-spin (HS) to low-spin (LS) transition of ferrous iron (Fe${}^{2+}$) in (Mg,Fe)CO${}_3$ ferromagnesite has been indicated via Raman spectroscopy, X-ray emission spectroscopy, and X-ray diffraction measurements at ∼50 GPa [14,17,18,19,20,21]. In addition, first principles calculation by Hsu and Huang [22] showed that a spin transition from HS to LS occurs in ferromagnesite at around 45–50 GPa, which is consistent with the previously conducted experiments. Furthermore, Hsu and Huang [22] also found that spin transition pressure does not depend on the Fe concentration in ferromagnesite.

In general, it is believed that the spin transition of Fe significantly affects the structure, elasticity, and thermodynamic properties of minerals owing to a large difference between the ionic radii of Fe in its HS and LS states. However, these effects of Fe spin transitions on the stability of high-pressure phases of ferromagnesite have not been investigated thus far. Thus, in this study, we investigated the structure and stability of ferromagnesite and its possible high-pressure phases via first principles calculation. Because the strongly correlated nature of Fe cannot be properly reproduced using local density approximation (LDA) and generalized gradient approximation (GGA) approaches, we used the internally consistent LDA+U method wherein the Hubbard parameter U was determined based on our first principles calculation [23].

2. Results

2.1. Calculation of MgCO${}_3$ Polymorphs under Pressure

First, we optimized the structural parameters of R$\overline{3}$c, P$\overline{1}$, and C2/m structures of MgCO${}_3$ polymorphs and confirmed that our calculations are consistent with previously reported results (

Table 1. Optimized structural parameters of P$\overline{1}$ and C2/m phases of MgCO${}_3$ polymorphs.

P$\overline{1}$ structureat 100 GPa
a = 5.178 Å, b = 5.206 Å, c = 7.242 Å, $\alpha =69.85°,\beta =81.64°,\gamma =78.04°$
(Pickard and Needs [12]: a = 5.211 Å, b = 5.238 Å, c = 7.268 Å, $\alpha =70.030°,\beta =81.904°,\gamma =78.272°$ at 100 GPa)
Atom	Wyckoff position		Atomic coordinates
		x	y	z
Mg1(M1)	2i	0.7414	0.2541	0.0022
Mg2(M2)	2i	0.5669	0.4731	0.3143
Mg3(M3)	2i	0.0636	0.1784	0.6501
C1	2i	0.2237	0.2078	0.0295
C2	2i	0.0812	0.3351	0.3092
C3	2i	0.4398	0.9790	0.3231
O1	2i	0.1277	0.1522	0.8965
O2	2i	0.0086	0.2930	0.1470
O3	2i	0.8794	0.3896	0.4232
O4	2i	0.3461	0.9572	0.1548
O5	2i	0.3752	0.3872	0.9623
O6	2i	0.2030	0.0682	0.4149
O7	2i	0.2249	0.5252	0.2622
O8	2i	0.5465	0.7397	0.4278
O9	2i	0.6110	0.1433	0.2692
C2/m structure at 120 GPa
a = 8.0417 Å, b = 6.4468 Å, c = 6.8273 Å, $\beta$ = 103.84°
(Oganov et al. [11]: a = 8.0945 Å, b = 6.4881 Å, c = 6.8795 Å, $\beta$ = 103.98° at 120 GPa)
Atom	Wyckoff position		Atomic coordinates
		x	y	z
Mg1(M1)	4g	0.0000	0.7536	0.0000
Mg2(M2)	4i	0.3229	0.5000	0.6975
Mg3(M3)	4i	0.9353	0.5000	0.6555
C1	8j	0.3701	0.3214	0.3273
C2	4i	0.2315	0.5000	0.0363
O1	4i	0.0926	0.5000	0.8961
O2	8j	0.3490	0.1665	0.4299
O3	8j	0.5081	0.3095	0.2675
O4	4i	0.3661	0.5000	0.9729
O5	4i	0.3522	0.5000	0.4276
O6	8j	0.2253	0.3320	0.1595

) [11,12]. In particular, the coordination number (CN) of all Mg sites is six in R$\overline{3}$c structure, whereas P$\overline{1}$ and C2/m structures contain three different crystallographic Mg sites (Figure 1). Some of the CNs of Mg in P$\overline{1}$ and C2/m structures are larger than six (P$\overline{1}$: CN${}_{M1}$ = 8, CN${}_{M2}$ = 6, CN${}_{M3}$ = 8; C2/m: CN${}_{M1}$ = 8, CN${}_{M2}$ = 8, CN${}_{M3}$ = 10), and thus their polyhedral volumes are significantly larger than that of R$\overline{3}$c (V${}_M$ = 9.0∼12.0 Å${}^3$) (Figure 1).

2.2. Structure Modeling of (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$

For the atomic configuration of Fe in (Mg,Fe)CO${}_3$, we selected the lowest enthalpy configuration of Fe in the cell. In particular, we used the 2 × 2 × 1 supercell of the conventional R$\overline{3}$c structure, and replaced four Mg with Fe atoms (Figure 2) to avoid the generation of an FeO${}_6$ octahedral layer. Consequently, there was almost no site preference for Fe configurations in the R$\overline{3}$c-type structure because all Mg sites in the R$\overline{3}$c structure have a CN of six with a Wyckoff 6b position. It should be noted that we retain the same crystal structure notations for (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ as those of MgCO${}_3$ though the symmetry (space group) of (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ polymorphs are changed to $P1$ (no symmetry) because of the incorporation of Fe.

In the P$\overline{1}$ structure, there are three independent Mg positions (Figure 1, Table 1); in addition, we confirmed that the structure with Fe in the M3 position shows the lowest enthalpy in the HS state structures, whereas Fe in the M2 position is the most stable in the LS state structures (Figure 3). Owing to the relatively larger ionic radius of HS-Fe, the structure with a large polyhedral volume (V${}_{M3}$ = 14.10 Å${}^3$ at 100 GPa) is most stable, whereas, for LS-Fe, the structure with the smallest polyhedral volume (V${}_{M2}$ = 7.87 Å${}^3$ at 100 GPa) is most stable (Figure 2).

A similar trend was observed for site preference of Fe in the C2/m crystal structure, wherein HS-Fe atom is the most stable at the M3 site, which has the largest polyhedral volume, whereas LS-Fe atom is most stable at the M2 site leading to the lowest enthalpy (Figure 4).

2.3. High Pressure Transitions of MgCO${}_3$ and (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$

We have calculated the enthalpy of MgCO${}_3$ and (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ with R$\overline{3}$c, P$\overline{1}$, and C2/m structures. Figure 5a shows the enthalpy of MgCO${}_3$ relative to that of the C2/m phase based on LDA. MgCO${}_3$-magnesite transforms into P$\overline{1}$ phase at 60 GPa, and the C2/m phase is stable above 80 GPa. In contrast, calculations in previous studies indicate the stability of the P$\overline{1}$ phase between 85–101 GPa, and that of the C2/m phase between 101–144 GPa based on GGA [12]. It is typically known that LDA/GGA leads to overestimation/underestimation of the bond strength; therefore, the volume calculated using LDA is usually smaller than that using GGA. Our calculated transition pressures for MgCO${}_3$ obtained using LDA are 10–15 GPa lower than those obtained using GGA; however, these GGA results are consistent with previously reported GGA results (

Table 2. Comparison between previously reported calculation results for high-pressure transition of MgCO${}_3$ and those obtained in this study.

	GGA [12]	GGA (This Study)	LDA (This Study)
R$\overline{3}$c → P$\overline{1}$	85	85	62
P$\overline{1}$ → C2/m	101	95	80
C2/m → P2${}_1$2${}_1$2${}_1$	144	145	130

) [12]. Hereafter, in our study, we use LDA results to investigate the stability of ferromagnesite. Nevertheless, it should be noted that the difference of about 10–15 GPa in the transition pressures between LDA and GGA results is typical, and thus, it provides lower and upper bounds for the transition pressures, respectively.

Figure 5b shows the relative enthalpies of ferromagnesite (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ and those of its high-pressure candidate phases with respect to the HS-C2/m phase obtained using the LDA+U method. The spin crossover of ferromagnesite (from HS-R$\overline{3}$c to LS-R$\overline{3}$c) occurs at ∼50 GPa which is consistent with previously reported calculations [22]. However, the HS-P$\overline{1}$ phase has a lower enthalpy than LS-ferromagnesite above ∼50 GPa. In addition, the HS-C2/m phase has a slightly lower enthalpy than that of the HS-P$\overline{1}$ phase above about 80 GPa, and the former is stable up to about 120 GPa, suggesting that no spin crossover is expected under lower mantle pressure conditions. This high stability of the HS state in P$\overline{1}$ and C2/m phases can be attributed to the large polyhedral volumes due to Fe${}^{2+}$ in these structures (Figure 1 and Figure 2). However, the energy differences between HS and LS states are relatively small (∼0.03 eV, equivalent to ∼350 K) in the P$\overline{1}$ and C2/m phases above 100 GPa. Therefore, there is still a possibility that the high-pressure transition sequence changes or spin crossover occurs around 100 GPa for above-ambient temperature conditions.

3. Discussion

Carbonate minerals are expected to exist in the subducting cold slabs. Because SiO${}_2$ is one of the most abundant components in a subducting slab, it has been previously indicated that a reaction between MgCO${}_3$ and SiO${}_2$ can occur in deeply subducted slabs [8]. We calculated enthalpies of the following chemical reactions to investigate the relative stability of carbonate phases compared with the mantle mineral phases:

$$MgC{O}_3+Si{O}_2\to MgSi{O}_3+C{O}_2,$$

(1)

wherein the high pressure phase of CO${}_2$ corresponds to the tetragonal I4$\overline{2}$d phase [11], SiO${}_2$ is stishovite, and the CaCl${}_2$-type structure below and above 50 GPa, respectively, and MgSiO${}_3$ has bridgmanite-type structure. The relative enthalpy curves shown in Figure 6a indicate that the transformation of MgCO${}_3$ and SiO${}_2$ into MgSiO${}_3$-bridgmanite along with CO${}_2$ formation is unfavorable, which is consistent with the results reported in previous studies [11,12].

Then, we calculated the reaction between (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ and SiO${}_2$,

$$\left(M{g}_{0.833}F{e}_{0.167}\right)C{O}_3+Si{O}_2\to \left(M{g}_{0.833}F{e}_{0.167}\right)Si{O}_3+C{O}_2.$$

(2)

For the calculation of enthalpy of (Mg${}_{0.833}$Fe${}_{0.167}$)SiO${}_3$, we used the supercell of bridgmanite containing 60 atoms with high spin Fe${}^{2+}$ at A sites [24]. Although both Fe${}^{2+}$ and Fe${}^{3+}$ have been observed in bridgmanite, ferric iron Fe${}^{3+}$ content in bridgmanite strongly correlates with the concentration of Al in bridgmanite [25]. In addition, the spin crossover of Fe${}^{2+}$ from HS to LS at an A site has not been theoretically shown at lower mantle pressure conditions (e.g., [26]). Moreover, the enthalpy of HS-(Mg${}_{0.833}$Fe${}_{0.167}$)SiO${}_3$-bridgmanite plus CO${}_2$ is higher than that of the Fe-bearing carbonate mineral plus SiO${}_2$ up to ∼120 GPa, which also suggests that decomposition of (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$-carbonate phases is unfavorable under lower mantle pressure conditions (Figure 6a).

Finally, we also investigated the stability of Fe${}^{2+}$ between bridgmanite and carbonate minerals by the following Fe exchange reaction,

$$\left(M{g}_{0.833}F{e}_{0.167}\right)C{O}_3+MgSi{O}_3\to MgC{O}_3+\left(M{g}_{0.833}F{e}_{0.167}\right)Si{O}_3.$$

(3)

In the case of R$\overline{3}$c structures, the enthalpy of the mixture of (Mg,Fe)CO${}_3$ and MgSiO${}_3$ is higher than that of MgCO${}_3$ and (Mg,Fe)SiO${}_3$ (Figure 6b), which indicates that Fe${}^{2+}$ is more stable in bridgmanite than in (Mg,Fe)CO${}_3$. HS-Fe${}^{2+}$ is more stable at the A site of bridgmanite because the polyhedral volume of the A site of the bridgmanite (16–21 Å${}^3$) is significantly larger than those of ferromagnesite (7–12 Å${}^3$). However, the HS-Fe${}^{2+}$ is more stable in P$\overline{1}$ and C2/m structures because HS-Fe${}^{2+}$ can be well accommodated at the sites with large polyhedral volumes in these phases.

4. Conclusions

In this study, we calculated the enthalpy of ferromagnesite (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ to investigate the effect of Fe on the stability of high-pressure phases of magnesite. Our enthalpy calculations suggest that (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ undergoes phase transition from the HS-R$\overline{3}$c structure to the HS-P$\overline{1}$ structure at 50 GPa, and to the HS-C2/m structure above 80 GPa at static 0 K conditions. Therefore, no spin transition in (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ is expected under the lower mantle pressure conditions. This result also indicates that the effect of Fe on the transition pressure of high-pressure phases of ferromagnesite is not significant.

Furthermore, the calculated enthalpy relationship suggests that the high-pressure phase of ferromagnesite plus SiO${}_2$ does not decompose even up to pressures of ∼120 GPa. In addition, the enthalpy of the Fe${}^{2+}$ exchange reaction between magnesite and bridgmanite shows a weaker preference in the high-pressure carbonate phases than in bridgmanite under deep lower mantle pressure conditions.

5. Methods

5.1. First Principles Calculation

We performed first principles calculation based on the density functional theory (DFT) [27,28] for the determination of the relative stability of (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ ferromagnesite and its possible high-pressure phases. We used the norm-conserving pseudopotentials [29] for Mg and Si, which have also been used extensively in previous studies e.g., [30]. In contrast, we used ultrasoft pseudopotentials for Fe, C, and O [31,32]. The electronic wave functions are expanded into plane-waves with a kinetic cutoff of 50 Ry. The irreducible Brillouin zone of R$\overline{3}$c (ferromagnesite), P$\overline{1}$, and C2/m phases were sampled on 2 × 2 × 2, 4 × 4 × 4, and 4 × 4 × 4 Monkhorst-Pack meshes with shifted grids, respectively [33]. The associated k-point sampling and kinetic energy cutoff are sufficiently converged within 0.001 eV/(Mg${}_{0.83}$Fe${}_{0.17})$CO${}_3$. All structure parameters of MgCO${}_3$ and (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ were fully relaxed at static 0 K and in the range of 0–120 GPa via the damped variable cell-shape molecular dynamics method implemented using a Quantum-espresso code until residual forces were less than 1.0 × 10${}^{-4}$ Ry/a.u.(a.u. = atomic unit) [34].

The number of atoms calculated in the MgCO${}_3$ cells are 30 for the R$\overline{3}$c (Wyckoff positions of Mg, C, and O are 6b, 6a, and 18e, respectively) and P$\overline{1}$ structures, whereas they are 60 for the C2/m structure (Wyckoff positions of P$\overline{1}$ and C2/m structures are listed in Table 1). For the (Mg${}_{0.833}$Fe${}_{0.167}$)CO${}_3$ composition, we constructed a supercell containing 120, 60, and 30 atoms for the R$\overline{3}$c, P$\overline{1}$, and C2/m structures, respectively (Figure 2).

5.2. LDA+U${}_{ic}$ Calculation

Because LDA and GGA typically fail to accurately determine the phase relationship of Fe-bearing mineral oxides owing to the large on-site Coulomb interactions of 3d Fe electrons, in this study, we performed calculations based on LDA with Hubbard U correction, where U values were determined using an internally consistent method (U${}_{ic}$, Figure 7) [23]. This method has been successfully applied for mantle mineral stability calculations [24,35,36]. The U${}_{ic}$ parameters for HS state of (Mg,Fe)CO${}_3$ are obtained in the range of 3.5–4.0 eV, except for ferromagnesite (HS-R$\overline{3}$c), for which they are higher than those reported in a previous study [22]. The large U${}_{ic}$ of the HS ferromagnesite could be attributed to the Fe atom occupying the small octahedral site in the R$\overline{3}$c structure. In contrast, the U${}_{ic}$ parameters of the LS state are in the range of 5.0–7.0 eV.