Article Cite This: ACS Earth Space Chem. 2018, 2, 1−8
Structures and Transport Properties of CaCO3 Melts under Earth’s Mantle Conditions Xiangpo Du,† Min Wu,*,‡ John S. Tse,*,†,§ and Yuanming Pan*,∥ †
Key State Laboratory of Superhard Materials, Jilin University, Changchun, 130012, P.R. China College of Materials Science and Engineering, Zhejiang University of Technology, Hangzhou, 310014, P. R. China § Department of Physics and Engineering Physics, University of Saskatchewan, Saskatoon, Saskatchewan S7N 5E2 Canada ∥ Department of Geological Sciences, University of Saskatchewan, Saskatoon, Saskatchewan S7N 5E2 Canada ‡
ABSTRACT: Carbonatitic and carbonated silicate melts generated from melting in the mantle are the chief agents for liberating carbon from the solid Earth and exert important controls on Earth’s deep carbon cycle. However, significant gaps in our knowledge about the carbonatitic melts under conditions pertinent to Earth’s deep interior remain due to experimental challenges. Here, we report on a first-principles molecular dynamics (FPMD) calculation of calcium carbonate (CaCO3) melts at pressures up to 52.5 GPa. Our FPMD calculations reproduce the ultralow viscosity measured by experiments and confirm the ideal liquid behavior of calcium carbonate melts at pressures below 11.2 GPa. However, calcium carbonate melts are characterized by a pronounced nonideal liquid behavior at pressures above 11.2 GPa, arising from (1) the temporal formation of small carbonate clusters and (2) increased interactions between Ca2+ and CO32− ions. It is found that the Stokes−Einstein equation relating the viscosity with the diffusion coefficients still holds at high pressure provided that a suitable effective particle size can be chosen. KEYWORDS: carbonate, Earth’s mantle, viscosity, melt, molecular dynamics
■
INTRODUCTION Calcium carbonate is the main component of skeletons of living organisms and also occurs as the dominant minerals in limestones and marbles, as well as a common phase in many rocks, constituting an enormous carbon reservoir.1 At ambient conditions, the two common crystalline phases of calcium carbonate are calcite and aragonite. Aragonite is metastable with respect to calcite and is the stable form above ∼650 K.2 Under ambient pressure, calcite decomposes into CaO and CO2 above 1170 K.3,4 High pressure experiments in laserheated diamond anvil cells (DAC) have shown that calcium carbonate melts congruently at least up to 3400 K when the pressure is lower than 43 GPa. Moreover, compressed calcium carbonate melts decompose into CaO and CO2 at temperatures higher than 3900 K.29,30 In addition, the CO2 produced from the decomposition of carbonate melts disassociates further into molecular oxygen and perhaps to form diamond at 30−80 GPa.5,6 Carbonated melting that generates carbonatitic and carbonated silicate melts in the mantle has been proposed to be the chief agents for liberating carbon from the solid Earth and thus exerting important controls on the Earth’s deep carbon cycle.29,30 Multiple lines of evidence from inclusions in superdeep diamond, high P-T melting experiments and theoretical calculations suggest that carbonatitic melts are present down to the lower mantle.7−10 Therefore, knowledge about the structures and properties of carbonatitic melts under pressure not only has significant impact on the geological processes research from mantle geodynamics to volcanic activities but © 2017 American Chemical Society
also important for the understanding of the related deep carbon cycle and global climate change.11−15 Compared to the much studied silicate melts, the physical properties of carbonate melts, in particular under extreme conditions (high pressure and high temperature), are scarcely investigated.12,15−18 Recently, an ultralow viscosity of carbonate melts (i.e., about 1 order of magnitude less than calcium silicate melts) up to 6.2 GPa has been measured experimentally.19 The higher viscosity of silicate melts is attributable to the polymerized SiO4 units that impede diffusion. The observed ultralow viscosity of calcium carbonate melts was supported by a molecular dynamics (MD) calculation using an empirical potential fitted to first-principles calculations.17 However, because of the difficulty of in situ high temperature and high pressure experiments, the viscosity of calcium carbonate melts under conditions relevant to Earth’s transition zone and lower mantle has not been studied. To better understand the transport properties of calcium carbonate melts in Earth’s deep interior, experimental investigations at higher pressures are needed but remain technically challenging. In this Article, we extended the investigation of the structures and transport properties of calcium carbonate melts to Earth’s lower mantle conditions. Received: Revised: Accepted: Published: 1
September 6, 2017 November 8, 2017 November 21, 2017 November 21, 2017 DOI: 10.1021/acsearthspacechem.7b00100 ACS Earth Space Chem. 2018, 2, 1−8
ACS Earth and Space Chemistry
■
RESULTS AND DISCUSSION Equation of States of Crystalline and Molten Calcium Carbonates. In calcite (Figure 1a), triangular CO3 anions
It is well-known that high pressure can alter the electronic structures of a material dramatically, and consequently, changes its geometrical structures and related physical properties. For instance, elemental calcium is known to transform from face centered cubic (FCC) to body centered cubic (BCC) structures at around 8 GPa.20 Theoretical work has also shown that orbital hybridization between Ca and Al starts to emerge at ∼12 GPa in compressed Ca−Al metallic glasses.21 Therefore, it is often not possible to extrapolate the property information on minerals and melts obtained at low pressures to higher pressures. To circumvent the uncertainty of MD calculations employing force fields parametrized at ambient conditions, first-principles molecular dynamic (FPMD) calculation has been proven to be able to correctly describe the interatomic interaction at high temperature and high pressure. The goals of the contribution are to unravel the mechanisms for the viscosity evolution and accompanied structural changes of calcium carbonate melts at high pressures.
■
Article
Figure 1. Structure of calcium carbonate polymorphs: (a) calcite. (b) A snapshot of molten calcium carbonate at 2000 K and 0.6 GPa. The gray, brown, and red balls represent Ca, C, and O atoms, respectively.
METHODS formed from sp2 hybridized C with three oxygen atoms were arranged to form planes where a layer of Ca atoms were intercalated between these planes. To evaluate the accuracy of the parameters used in the first-principles calculations, the theoretical densities of calcite are compared to those from experiments and previous calculations in Table 1. For both calcite and aragonite, the present theoretical calculations gave the best densities as compared to the experimental values. The density of calcite calculated with the Ca_sv potential at ambient pressure is 2.75 g/cm3 and agrees well with the experimental observation of 2.72 g/cm3.32 In comparison, the calculated density of 2.78 g/cm3 from the Ca_pv potential is slightly overestimated. The density of aragonite calculated with the Ca_sv potential of 3.00 g/cm3 agrees well with the experiment observation of 2.93 g/cm3 and, once again, the density of 3.04 g/cm3 obtained from the Ca_pv potential is slightly overestimated.33 These results show the importance of including the inner-core 3s2 electrons in the calculations even at ambient pressure. It is noteworthy that all previous theoretical calculations employing different GGA functionals underestimated the densities of both calcite and aragonite (Table 1). For instance, the densities of calcite and aragonite predicted by Vuilleumier et al.17 using GGA (BLYP) together with Van der Waal dispersion at 0 K are ∼2% smaller than those from experimental results at ambient condition and ∼3% smaller than our calculated result using the Ca_sv potential at 0 K, respectively. The results demonstrate that the structures predicted by the present choice of computational parameters are more reliable than earlier works. From Bader atomic charge analysis,39 each cationic Ca loses 1.7 electrons to the CO3 anion. The electronic density of states shows calcite is an insulator with an energy gap of 5 eV. The highest occupied crystal orbitals (HOCO) were composed mainly of O lone pair electrons. This description is similar to that from previous theoretical calculation.40 In the melts (Figure 1b), the arrangement of CO3 group is more disordered and the bandgap energy decreases dramatically to 1.2 eV. There is no spectroscopic study on the of CaCO3 melt. The evaluation on the accuracy of the computed density of molten calcium carbonate at high temperature and high pressure is more difficult due to lack of experimental data. A measurement on the CaCO3−Li2CO3−Na2CO3−K2CO3 quaternary system
FPMD calculations and electronic structure calculations were performed using the Vienna Ab Initio Simulation Package (VASP) based on the density functional theory.22,23 Effects of the core electrons were replaced by atomic projected augmented wave (PAW) potentials.24 Valence electronic configurations of 2s22p2 and 2s22p4 were used for C and O atomic potentials, respectively. To better describe the electronic structure of Ca atoms under pressure, a hard Ca potential with valence electronic configurations of 3s23p64s2 (Ca_sv) was used in the present study. Calculation using a softer Ca potential neglecting the 3s orbital (Ca_pv) was also performed for comparison. The Perdew−Burke−Ernzerhof (PBE) exchangecorrelation functional was used.25 The wave functions were represented by a plane-wave basis set with an energy cutoff of 400 eV. A molten structure of CaCO3 at ambient pressure was generated by melting a 2 × 2 × 1 supercell constructed from the calcite unit cell containing 120 atoms (24 CaCO3 formula units) by heating to 3000 K. The liquid structure at ambient pressure and 2000 K was obtained by cooling slowly from 3000 K in several temperature steps. Higher pressure liquid structures were obtained by isotropically reducing the volume of the liquid model obtained at ambient pressure. The pressure was then determined from the stress tensor. Shear viscosities were computed from the Green−Kubo autocorrelation function of the off-diagonal components of the stress tensor (SACF) in the canonical ensemble (NVT).26−28 The NoseHoover thermostat was used to control the temperature.27 A single Γ point was used for Brillouin zone sampling. A very long trajectory is required to achieve convergence in the calculated viscosity. Therefore, at each pressure, MD calculation ran for at least 100 ps. Since viscosity is usually difficult to measure at high pressure. Most often the experimental viscosity was estimated from the use of the Stokes−Einstein relationship from the diffusion constant, which is easier to determine using X-ray tomography.29,30 This simple relationship, however, has not been validated theoretically for high pressure. Therefore, we compare the viscosities calculated from the Green−Kubo correlation function31 with those estimated from the Stokes− Einstein equation from the corresponding diffusion coefficient computed from the same MD trajectory. 2
DOI: 10.1021/acsearthspacechem.7b00100 ACS Earth Space Chem. 2018, 2, 1−8
Article
ACS Earth and Space Chemistry
Table 1. Comparison of the Densities (g/cm3) of Calcite and Aragonite Crystalline Structures by First-Principles Calculations and Experiments calcite aragonite a
exp.
Ca_sv
Ca_pv
BLYP+D2
B3LYP
PBE
PBE
2.72a 2.93b
2.75c 3.00c
2.78c 3.04c
2.67d 2.83d
2.604e 2.821f
2.603g 2.851h
2.650i
Ref 32. bRef 33. cThis work. dRef 17. eRef 34. fRef 35. gRef 36. hRef 37. iRef 38.
Figure 4. (a) Calculated pair distribution functions at 3000 K, 0.7 GPa. (b) Effective particle diameters at 2000 K (red) and 3000 K (black) calculated from the equation (dCC + 2dCCa + dCaCa)/4.
Figure 2. (a) Pressure-dependent viscosities of molten calcium carbonates at 2000 and 3000 K. The experiment result at 2063 K and 6 GPa was marked as a star. The solid red and black circles represent the calculated viscosity from the Green−Kubo relationship at 2000 and 3000 K, respectively. The open pink and blue circles represent the calculated viscosity from the Stokes−Einstein (SE) equation using the theoretical diffusion coefficients (see text) at 2000 and 3000 K, respectively. The dashed lines are guides for eyes. (b) The SACF (black solid line) and the integral (blue solid line) giving the viscosity calculated at 3000 K and 32.6 GPa. The pink dashed line is a guide for eyes.
0.1 GPa. These values are close to the present calculated density of ∼2.52 g/cm3 at 2000 K and 0.7 GPa. On the other hand, the FPMD calculation by Vuilleumier et al.17 predicted a density of 2.94 g/cm3 at 2073 K and 12.2 GPa, which is close to the present value of ∼3.1 g/cm3 at 12 GPa and 2000 K. The agreement is even more encouraging taking into the consideration that the densities of calcite and aragonite were underestimated from previous theoretical studies as discussed above (vide infra). Moreover, the different exchange-correlation functionals used in the calculations and the procedures to generate the liquid structure and the length of the MD run may also contribute to the discrepancies. We conclude that the calculated densities of the melts in the present study are reliable. Viscosity of Molten Calcium Carbonate under Pressure. Viscosity Calculated from Green−Kubo Relation. The goals of this study are to investigate the CaCO3 melts and to test the validity of the Stokes−Einstein relationship29,30 at high pressure. The viscosities can be calculated from firstprinciples from the time correlation function of the temporal stress tensors from the MD calculation via the Green−Kubo relation,31
Figure 3. Variation of calculated diffusion coefficient D with pressure at 2000 and 3000 K.
η=
by Liu et al.41 suggested the density of CaCO3 melts to be ∼2.6 g/cm3 at 1073 K and 0.7 GPa. From thermodynamic data, Treiman et al.42 estimated a density of 2.4 g/cm3 at 1596 K and
V kBT
∫0
∞
⟨Pxy(t )Pxy(0)⟩dt
(1)
where η is the viscosity, V is the volume of the system, T is the temperature, kB is the Boltzmann constant, ⟨⟩ denotes the 3
DOI: 10.1021/acsearthspacechem.7b00100 ACS Earth Space Chem. 2018, 2, 1−8
Article
ACS Earth and Space Chemistry
Figure 5. Comparison of the pressure dependent pair distribution function (PDF) computed at 2000 and 3000 K.
ensemble average, and Pxy is the off-diagonal component of the instantaneous stress tensor. In order to maintain comparable viscosity at high pressure in the experiment, it was necessary to increase the temperature. For example, the experiment of Kono et al.19 showed that a measured viscosity of ∼6 mPa·s at 0.9 GPa and 1380 °C (1653 K) required an increase of temperature to 1790 °C (2063 K) at 6.2 GPa. Compared to the experiment, the calculated viscosity at 0.7 GPa and 2000 K is ∼3.8 mPa·s. To maintain a similar viscosity at a higher pressure of 7 GPa, it was necessary to increase the temperature to ∼3000 K. The predicted trend is consistent with experimental observations although the absolute temperature difference is higher. The calculated pressure dependent viscosities at 2000 and 3000 K are shown in Figure 2. To illustrate the convergence of the viscosity calculation, the SACF and the integral (viscosity) as a function of the simulation time at 32.6 GPa were presented in Figure 2b. The calculated viscosity at 4.1 GPa and 2000 K of 5.8 mPa·s agrees with the experimental value of ∼6 mPa·s measured at 2063 K and 6.2 GPa.19 However, at 2000 K and above 8 GPa, calcium carbonate starts to solidify and becomes a viscous liquid with very low self-diffusion coefficients and very long simulation time is needed. Therefore, we can only predict reliably the viscosity of the melt below this pressure. Since the
experimental melting temperature of CaCO3 is close to 2000 K, the very slow relaxation of the CaCO3 at 2000 K is indicative the melting temperature has increased with pressure. In comparison, the calcium carbonate prepared at 3000 K is more fluid like and accurate viscosities can be computed up to 32.6 GPa. The viscosity at 3000 K and 0.6 GPa is 2 mPa-s and incidentally, close to the water viscosity of 0.89 mPa·s at ambient conditions. The result confirms the experimental finding. Furthermore, the calculated viscosity at 3000 K increases slowly at low pressure (11.2 GPa) and may signal a structural change (vide supra). Self-Diffusion Coefficient. The self-diffusion coefficients D for the ions in calcium carbonate melts were calculated from the Einstein mean square displacement formula N
D(t , t → ∞) =
⟨ri(t ) − ri(0)⟩2 1 ∑ N i=1 6t
(2)
where ri is the position of the ion i and the bracket represents the ensemble average over many time origins. The results at 2000 and 3000 K are shown in Figure 3a and 3b with the 4
DOI: 10.1021/acsearthspacechem.7b00100 ACS Earth Space Chem. 2018, 2, 1−8
Article
ACS Earth and Space Chemistry
should be noted that an alternative relationship Deff = [DCa×DCO3]/[0.5DCa + 0.5DCO3] has been suggested to be suitable for mixed charged species. In this case, since DCa is close to DCO3, therefore both expressions gave similar effective diffusion coefficients. As shown in Figure 3, the calculated diffusion coefficient of the CO32− anion of ∼10−9 m2/s at 2000 K at low pressure is similar to previous theoretical studies.17,43 The D values at both 2000 and 3000 K decrease significantly upon compression at low pressure, indicating that at higher pressure (>4 GPa at 2000 K, >15 GPa at 3000 K), the ions are much more restrictive with compression. Viscosity Calculated from the Stokes−Einstein Equation. In high-pressure experimental studies, the liquid viscosity is often estimated from the diffusion coefficient using the Stokes− Einstein equation ηSE = (kBT)/(2πDd), where ηSE is the viscosity, D is the self-diffusion coefficient, d is the particle diameter, T is the temperature, and kB is the Boltzmann constant. However, the Stokes−Einstein equation is valid only in an ideal fluid composed of spherical particles that have no slippage against the layer of the fluid in contact with it.29 It is well-known that the Stokes−Einstein equation has to be corrected for fluids with nonspherical particles because the friction coefficient is expected to be different. The viscosity calculated from the Stokes−Einstein equation also depends on how the “particle” diameter is determined. In the calcium carbonate melts, the diffusing particles are charged Ca2+ and CO32− ions with different shapes and sizes. The condition is very different from an ideal fluid. For instance, the Stokes− Einstein relation has been shown to fail in describing the viscosity of silicate liquids.44 Therefore, it is a somewhat surprising that the Stokes−Einstein relation seems to work quite well for the calcium carbonate melts below 15 GPa in the study of Vuilleumier et al.17 To examine the Stokes−Einstein relation in a broader pressure range, we employed the same procedure described earlier to determine the particle diameters: viz. the effective particle diameter d was calculated from the equation (dCC + 2dCCa + dCaCa)/4, where dCC, dCCa, and dCaCa are the first maximum of the respective pair distribution functions (PDFs) obtained from MD calculations. For this case, the results at 0.7 GPa at 3000 K are shown in Figure 4a. Instead of using a constant particle diameter for all the viscosity calculations as assumed in the previous work, we employed the pressure dependent particle diameters derived from the respective pair distribution functions calculated from the MD trajectories. The results are shown in Figure 4b. A comparison of the estimated ηSE from the Stokes−Einstein equation and ηGK computed directly from the Green−Kubo relationship is shown in Figure 2. The calculated ηSE values at 2000 K agree well with ηGK up to 7 GPa. This observation is consistent with a previous theoretical study.17 Intuitively, it is not expected that the ηSE values at 3000 K would also agree well with ηGK at pressures up to 32.6 GPa, even though there is an obvious abrupt change in the viscosity at ∼11.2 GPa which may indicate a structural transformation (vide supra) has occurred. The present results show that the Stokes−Einstein relationship is valid for the calcium carbonate melts within the pressure range studied here as long as an appropriate method is used to estimate the particle size. In retrospect, the success of the Stokes−Einstein equation with an appropriate variable “particle size” at different pressures is not too surprisingly. The FPMD calculations show that below
Figure 6. Pressure-dependent partial coordination numbers of atoms in calcium carbonate melts at 3000 K. (a) Number of O atoms surrounding a selected C atom and (b) number of C atoms surrounding a selected O atom.
Figure 7. Instantaneous formation of small CO3 clusters in the calcium carbonate melt at 3000 K and 32.6 GPa.
estimated errors. The diffusion coefficients of Ca2+ and CO32− are very similar. The results suggest both Ca2+ and CO32− diffused in unison. The diffusion coefficients doubled from ∼0 GPa when the temperature was raised from 2000 to 3000 K. However, at both temperatures, the diffusion coefficients decrease to ∼0.1 × 108 m2/s at high pressure (8 GPa at 2000 K and 25 GPa at 3000 K). We have followed the suggestion from a previous study17 that the effective diffusion coefficients D of the system were computed from the arithmetic mean of DCa and DCO3 as the CO3 diffuse as a unit and no decomposition was observed in the P-T range studied here. It 5
DOI: 10.1021/acsearthspacechem.7b00100 ACS Earth Space Chem. 2018, 2, 1−8
Article
ACS Earth and Space Chemistry
Figure 8. Pressure-dependent velocity cross correlation functions of (a) Ca−C and (b) Ca−O in calcium carbonate melts, together with (c) the velocity density of states at 3000 K.
K and 32.6 GPa (Figure 7), although there is no extended network structure connected by the C−O covalent bonds, a few clusters connected instantaneously to two or even three CO3 units are formed to give rise to a few 4 coordinated carbons. The observation of these temporal poly carbonate clusters at high pressures support the previous report of Corradini et al.18 The change from the isolated and independent CO3 to correlated instantaneous clusters is responsible for the abrupt increase in the viscosity above 11.2 GPa. Ca and CO3 Interactions. Increasing interaction between Ca2+ and CO32− ions at high pressures also has an important impact on the viscosity. This interaction, as a function of pressure, can be analyzed from the corresponding velocity cross-correlation between Ca and O and between Ca and C atoms. The velocity cross-correlation functions (VCC) at 3000 K are depicted in Figure 8. It can be seen that the VCC profiles are liquid like at pressures below 7.1 GPa. Weak coupling between the Ca and the CO3 motions is reflected in the lack of variation in the viscosity of liquid CaCO3 in this pressure range. Above 11.2 GPa, however, the interactions between Ca and C as well as Ca with O increase noticeably. This is related to the temporal formation of small clusters as described above (vide infra). A band attributed to the Ca···CO3 vibrations starts to emerge in the plot of the Ca projected density of states (Figure 8c). This vibrational band corresponds to the three oscillation peaks observed at 11.2 GPa in VCC with a period of ∼71 fs, or a frequency of 470 cm−1. This frequency matches well the emergence of the new band in the Ca PDOS between 400 and 500 cm−1. We can conclude unambiguously from the analysis of the atomic dynamics that at high pressure (>11.2 GPa at 3000 K), the Ca2+ and CO32− ions are compressed closer to each other and their motions are correlated. As a consequence of the decrease in the diffusion rates, the viscosity of the melt increases by almost 1 order of magnitude. This observation may have important impact to the transport of melts under the mantle conditions.10,16 Moreover, the increased viscosity may affect the heat transport (thermal conductivity) in the mantle as well.
11.2 GPa, the structure of the carbonate melts, in particular at higher temperature (3000 K), does not changed fundamentally. As will be discussed in greater detail below, at 3000 K and below 11.2 GPa, there is no noticeable change in the CO distances of the carbonate anions indicating a more or less rigid group. Furthermore, the diffusions of the Ca2+ and CO32− ions are almost independent of each other. Therefore, the assumption of a constant particle size is appropriate in this pressure region. Above 11.2 GPa, an abrupt increase in viscosity is predicted. As will be shown below, the higher viscosity is due to increase in the correlation of the Ca2+ and CO32− diffusive motions and appearance of instantaneous clustering of CO32−. In this case, pressure-dependent corrections to the “particle” size are needed. We must caution that the Stokes−Einstein equation should not be regarded a universal relationship. Structures of the CaCO3 Melts. A comparison of the pair distribution functions (PDF) for 2000 and 3000 K MD calculations are shown in Figure 5. The Ca PDF at 2000 K and 0.7 GPa shows three principle features at 2.5, 4.2, and 6.6 Å. In comparison, the first peak of the Ca PDF at 3000 K is much broader indicating a more disordered nearest environment at higher temperature. The first peak of C PDF and O PDF at ∼1.3 Å can be assigned unambiguously to the C−O covalent bonds of the CO3 units.32,33 Interestingly, the positions of the peak are not sensitive to pressure, indicating that the CO32− anions are rigid because of the strong C−O covalent bonds. Compression only affects the distributions of the second or third nearest neighbors. A shoulder at ∼2.7 Å starts to emerge at high pressure in the O PDF for the structures at both 2000 and 3000 K. Inspection of the structure revealed that the distance between neighboring CO3 units decreases with increasing pressure. Upon closer examination, it is found that although the first peak of the O PDF is still separated from the second peak at 2000 K, the two peaks start to merge at higher pressures at 3000 K. At 2000 K and 40 GPa, the calculated coordination numbers of carbon and oxygen are 3 and 1, respectively. At 3000 K, the coordination numbers of both carbon and oxygen hardly change at pressure below 11.2 GPa (Figure 6). However, above 11.2 GPa at 3000 K, the coordination numbers of both carbon and oxygen atoms increase steadily. The increase in the CNs is best illustrated by the snapshot on the formation of instantaneous clusters at 3000
■
CONCLUSIONS In summary, the structures and viscosities of calcium carbonate melts under pressures up to 52.5 GPa have been investigated 6
DOI: 10.1021/acsearthspacechem.7b00100 ACS Earth Space Chem. 2018, 2, 1−8
ACS Earth and Space Chemistry
■
with the FPMD method. The present study reproduced the ultralow viscosity of calcium carbonate melts measured at low pressures (