Diffusion of CO2 in Magnesite under High Pressure and High Temperature from Molecular Dynamics Simulations

United Laboratory of High-Pressure Physics and Earthquake Science, Institute of Earthquake Forecasting, Chinese Earthquake Administration, Beijing 100036, China State Key Laboratory of Geological Processes and Mineral Resources, and School of Earth Sciences and Resources, China University of Geosciences, Beijing 100083, China Institute of Microstructure and Property of Advanced Materials, Beijing Key Lab of Microstructure and Property of Advanced Materials, Beijing University of Technology, Beijing 100124, China Faculty of Land Resources Engineering, Kunming University of Science and Technology, Kunming, 650093 Yunnan, China


Introduction
Prediction of the changes in the Earth's climate due to the increase of carbon dioxide (CO 2 ) concentrations needs to deeply understand the chemical and physical processes that control the carbon cycle [1]. Generally, the carbon cycle involves the deep and shallow subcycles of carbon. So far, intensive attention has been paid to the shallow carbon cycle that is mainly related to the atmosphere, oceans, and shallow crustal. However, there is relatively little known about the deep carbon cycle. Due to the Earth's deep interior that contains more than 90% of the total Earth's carbon [2,3], the carbon transport and storage in the deep Earth play a key role in determining the global carbon cycle [1,[4][5][6][7][8][9]. Therefore, it is very necessary to understand the carbon cycle in the Earth's deep interior.
Carbon in the mantle exists in the different oxidized forms (e.g., CO, CO 2 , and CO 3 2− ) and the different states (solid, molten carbonate, fluid, etc.) [10][11][12][13]. Magnesite (MgCO 3 ) is the most common carbonate phase in the subducted slabs [14]. High-pressure experiments suggest that magnesite with a rhombohedral structure can be stable up to 115 GPa at 2000~3000 K which may be the major carbon host in the mid-lower mantle [14][15][16][17]. CO 2 is an important carbon-bearing form in the Earth's interior. About 40 megatons of CO 2 emits from the mantle to the surface every year by diffuse outgassing [18]. Therefore, CO 2 plays a fundamental role in determining the global carbon cycle. In situ Raman spectroscopy in laser-heated diamond anvil cells revealed that CO 2 molecular was stable up to 34 GPa at~1770 K [6,19]. Therefore, the transport properties of CO 2 in the deep Earth under high pressure and high temperature are a vital link to understand the deep carbon cycle.
Except for high-pressure and high-temperature (high P-T) experiments, molecular dynamics (MD) simulations are powerful tools to investigate the molecular-scale dynamic properties of materials ( [20][21][22][23][24][25][26]; Hu et al. 018). There are some achievements on researching diffusion of CO 2 in different materials by MD simulation [27]. For example, the CO 2 diffusion coefficients of crystalline polystyrene at 1 atm. are 2:3 × 10 −7 cm 2 s −1 (at 333.15 K), 4:5 × 10 −7 cm 2 s −1 (at 343.15 K), and 7:7 × 10 −6 cm 2 s −1 (at 353.15 K) [28]. Diffusion coefficients of supercritical CO 2 and its mixtures monotonously decreased as its increasing molar density at 317.5 K [29]. The diffusion coefficients of H 2 , CH 4 , CO, O 2 , and CO 2 in water near the critical point (600~670 K, 250 atm) were investigated by MD simulation [30], and CO 2 can greatly influence the diffusion of mobile species in clay minerals [31]. However, most of those works focused on model structures, which are different from the potential minerals in the mantle. In addition, previous works were mainly performed under the static conditions without concerning high temperature and high pressure. Accordingly, the transport properties of CO 2 under high temperature-pressure in the Earth's interior conditions have not yet been identified. In this work, the diffusion of CO 2 in magnesite was investigated under the high temperature-pressure conditions mimicking the extreme environment in the mantle.

Crystalline Models.
Magnesite has a space group of R3c. Its structure is composed of the alternative layers of Mg and CO 3 groups. Among them, one C atom and three neighbor O atoms lie in the same plane that is parallel to the a-b plane. Under the ambient condition, the cell parameters of magnesite are as follows: a = b = 4:637 Å, c = 15:023 Å, α = β = 90°, γ = 120° [32]. A 6 × 6 × 2 supercell of magnesite was established. A vacuum layer was added between two CO 3 planes to build slit-like pores. The pore size (H, Å) was determined by the height of a vacuum layer. The twolayer slit-like pore structures of magnesite with electrically neutral were built and simulated ( Figure 1).

Equilibrium Adsorption Configurations: The Grand
Canonical Monte Carlo Method (GCMC) Simulation. The GCMC employed the μVT ensemble. The chemical potential μ, the volume (V), and the temperature (T) were kept constant, while the number of molecules N was allowed to fluctuate [33]. GCMC simulations were carried out by sorption code as implemented in Materials Studio 2017. The structures of CO 2 in magnesite pores with energy minimized adsorption were obtained under the different conditions. During the GCMC simulation, CO 2 molecules within the framework were randomly rotated and translated. Only were the positions and orientations of CO 2 molecules sampled. Each configuration was treated as a rigid body. The COM-PASS force field was used to describe the hypersurface of the potential energy on which the atomic nuclei moved [34]. The Ewald method was adopted to calculate the Coulomb force interaction. The atom interaction-based method was used to calculate van der Waals force with the Lennard-Jones potential. The cut-off distance was 0.15 Å. The number of load step in each calculated run is 6:0 × 10 6 , in which the balance step number is 3:0 × 10 6 and the process step number is 3:0 × 10 6 . This GCMC method was successfully used to calculate the absorption of CO 2 in magnesite in our previous works [35], and its validity also was proved by different works [31,[36][37][38][39][40].

Molecular Dynamics (MD)
Simulations. MD simulations were carried out by the Forcite code as implemented in Materials Studio 2017. The equilibrium configurations from GCMC calculations were further used in MD simulations. The MD simulations were performed with a time step of 1 fs (up to 1000 ps) using an NPT ensemble at 350~2500 K and 3~50 GPa. Interlayer molecular configurations were sampled every 1000 steps for further analysis. The temperature was kept constant using a velocity scale thermostat algorithm, and the convergence criterion is 1 K. The Ewald summation method was employed with an Ewald accuracy of 1 × 10 −5 kcal/mol. Like the GCMC simulations, the COM-PASS force field was used for describing the potential energy hypersurface [34].

Geofluids
where r i ! ðtÞ is the position of the atom I at time t, and the angular bracket h i denotes that the quantity is an ensemble average property. In an equilibrium ensemble, the average squared displacement is independent of the time t, which can be averaged out, leaving the mean square displacement over an interval Δt: If there are N equivalent particles, the MSD can be further averaged as If the particle is diffusing, the MSD becomes linear as a function of time (diffusive regime). The self-diffusion coefficient (D) can be described by the slope according to Einstein equation: where N is the number of diffusion molecules, t is the time, and r i ! ðtÞ is the displacement vector of the ith molecule at time t. The angular bracket h i denotes that the quantity is an ensemble average property.

The Diffusivities at High Pressure and High Temperature.
Diffusions of CO 2 in magnesite pore were calculated with the pore size of 25 Å and the CO 2 pressure of 50 MPa at high pressures (3~50 GPa) and high temperatures (350~2500 K); the results are shown in Figure 2. Diffusion of CO 2 in magnesite increased as temperature increases but decreased as pressure increases. Specifically, the lg D (D in 10 -12 m 2 /s) increased from 1.67 at 350 K to 4.45 at 2500 K (3 GPa), while it decreased from 4.45 at 3 GPa to 3.29 at 50 GPa (2500 K). However, the change trend of the diffusion with temperature is different under different pressures. At low pressure (<9 GPa), the diffusion coefficient consistently increased with increasing the temperature at the whole range explored in this work. However, as the pressure increased, the diffusion coefficients showed the different trends when the temperature is higher or lower than 900 K. For the temperature higher than 900 K, the lg D values showed a good linear relation with the 1/T, which corresponds well with the Arrhenius law: where D 0 is the preexponential factor, R is the gas constant, and E A is the activation energy for diffusion. Fitting parameters at 900~2500 K and 3~50 GPa are shown in Figure 3 and The different diffusion properties of CO 2 in magnesite under different pressure and temperature (9 GPa and 900 K) might be related to the phase transition of CO 2 and stability of magnesite under different pressure and temperature. The CO 2 might transform from fluid phase to Pa3 molecular structure, Pbcn, P42/mnm, and I42d with increasing pressure [41][42][43][44]; however, the stability of CO 2 under high temperature is not yet well understood. For example, at ambient temperature, Pa3 CO 2 undergoes a pressureinduced transformation to the Cmca symmetry at about 9-12 GPa; when pressure is larger than 12 GPa, there are two molecular phases: P42/mnm phase for T > 470 K and Pbcn phase for T > 500 K [43,44]. At the same time, the presence of CO 2 might affect the stability of magnesitite [45,46]. Figure 4 shows an anisotropic feature of CO 2 diffusion in magnesite. For the magnesite structure, it has the same a-axis and b-axis; therefore, the diffusion in the a -b plane of magnesite is the same. Accordingly, we only calculated diffusion coefficient along a-axis and c-axis direction under high T-P. The diffusion coefficients (D) of CO 2 along the a-axis and c-axis varied differently with temperature and pressure (Figure 4). At temperature of 350~500 K, the D a (D along a-axis direction) was larger than the D c (D along c-axis direction) at a low pressure (<9 GPa) and then the D a became smaller than the D c with increasing pressure. When temperature was higher than 700 K, the D a was larger than the D c , but the differences became very small with increasing pressure. Anisotropy of the CO 2 diffusion in magnesite decreased with increasing pressure and temperature. At 25 GPa, anisotropy increased from 3% at 2500 K to 60% at 350 K. For comparison, at 3 GPa, the anisotropy increased from 3.5% at 2500 K to 280% at 350 K. As the  4 Geofluids pressure was higher than 25 GPa, the difference between D c and D a was very small. The results indicated diffusion of CO 2 in magnesite is anisotropic at a low pressure and is isotropic under high pressure. Anisotropy of CO 2 diffusion in magnesite is important to understand the anisotropic properties of deep carbon cycle. A study on olivine-hosted melt inclusions from the Mid-Atlantic Ridge indicated that carbon content in the upper mantle is highly heterogeneous, varying by almost two orders of magnitude globally [47]. Volatiles in the mantle can be hosted either in independent phases or as dissolved components in minerals. The content of magnesite in the upper mantle is 1000~8000 ppm [17,48]. Under the upper mantle (temperature <~1900 K and pressure <~15 GPa), diffusion of CO 2 in magnesite shows an obvious anisotropy. Therefore, this diffusion difference of CO 2 probably leads to the heterogeneous transport and storage of carbon in the mantle.

The Diffusion under Different CO 2 Pressures.
The diffusion of CO 2 in magnesite is closely related to its content in the mantle. Therefore, D values of CO 2 in magnesite under different CO 2 pressures were calculated in order to understand the effect of CO 2 content on its diffusion behavior. Figure 5 shows the D values of CO 2 in magnesite with the CO 2 pressures of 10, 30, 50, 70, and 100 MPa under the 900~2500 K and 3~50 GPa. The diffusibility increased as the CO 2 pressure increased. The averages of D values increased by 1.5%-4.8% as the CO 2 pressure increased from 10 MPa to 100 MPa under the different pressures. As detected by Raman spectroscopy, the diffusion coefficient of CO 2 in aqueous solutions decreases significantly at high CO 2 concentrations [49]. Our results here show different diffusion behavior from the diffusivity of CO 2 in aqueous solutions under geologic carbon sequestration conditions. In addition, the effects of CO 2 pressure on its diffusion show a slight decrease with increasing the temperature. Compared with the change caused by temperature or pressure (orders of magnitude), the effect of CO 2 pressure on its diffusion is limited.  Figure 6 shows diffusions of CO 2 in these models. It is clear that D firstly increases and then decreases with increasing pore sizes. Under the same simulated condition of temperature and pressure, the diffusion has the highest ability in magnesite with a pore size of 20~25 Å. The differences of the diffusion caused by the change of pore size become smaller with increasing pressure.

Discussion
Plenty of CO 2 is emitted from the deep Earth to the surface by diffusion [18,47] and made significant influence on the

Geofluids
Earth's carbon cycle. Therefore, it is critical to understand the diffusion process of CO 2 in the Earth's interior. The diffusion coefficient of CO 2 in magnesite is 9~28000 × 10 −12 m 2 s −1 at 350~2500 K and 3~50 GPa. CO 2 diffusion in crystalline polystyrene was 2:3~77 × 10 −11 m 2 s −1 at 1 atm and 333.15~353.15 K [28]. The diffusion coefficients of CO 2 in   6 Geofluids silicate melt (Jadeite+albite+quartz) are 2:8~87 × 10 −12 m 2 s −1 at 1323~1673 K and 500 MPa [6]. Under the similar condition of temperature and pressure, our calculated diffusion coefficient of CO 2 is smaller than that in crystalline polystyrene, however is larger than that of CO 2 in silicate melt. The results indicate the faster diffusion ability of CO 2 in magnesite crystal than in silicate melt. Another critical issue on deep carbon cycle is the distribution of the diffusion ability of CO 2 along the Earth's geothermal profile. From Figure 7, we can estimate diffusion distance of CO 2 in magnesite is 0.017~0.933 m each year. Accordingly, CO 2 diffusion from the mantle to the surface takes around dozens of Ma. Carbon hosted in the mantle ultimately affected the surface carbon budgets over billionyear time scales [50,51]. Similarly, isotopic results indicated three time-scale carbon cycles in the entire upper mantle: 1: >120 km, 5~10 Ma; 2: 120~410 km, 48~53 Ma; and 3: 410~660 km, >60 Ma [52][53][54][55][56]. The diffusion of CO 3 2would be three to four orders of magnitude slower than that of molecular CO 2 [6]. Therefore, CO 2 diffusion plays the important role in determining the deep carbon cycles. For the first time, we present the CO 2 diffusion behavior in real minerals (magnesite) under the Earth's interior conditions. The results here provide new data to understand the deep carbon cycle.
The above results indicate that temperature, pressure, and pore size of magnesite have important effect on CO 2 diffusion. The excess adsorption amounts of CO 2 in the pores of magnesite decrease with increasing temperature but increase with increasing pressure and pore size [35]. So interior circumstance changes from deep to shallow of the Earth might drive the transport of CO 2 in its interior, such as the changes of stress and porosity caused by the movement of fault or tectonic slab. However, the grain or mineral boundaries are always the important channel for fluid migration in the Earth's interior, which also needs to be paid more attention.

Conclusion
For understanding the Earth's deep carbon cycle, diffusion of CO 2 in magnesite pore with different CO 2 pressure and pore sizes was calculated under 350~2500 K and 3~50 GPa. The diffusion coefficient of CO 2 in magnesite increases with increasing temperature but decreases with increasing pressure. The E A increases with increasing pressure, indicating the diffusion of CO 2 in magnesite becomes more difficult under high pressure. CO 2 has the strongest diffusion in the magnesite when the pore size is 20~25 Å, and the diffusion slightly increases with increasing CO 2 pressure. Diffusion of CO 2 in magnesite shows obvious anisotropy and might be responsible for the inhomogeneous distribution of carbon in the upper mantle. The diffusion coefficient of CO 2 in magnesite is in the range of 9 × 10 −12 m 2~2 8000 × 10 −12 m 2 s −1 under the pressure and temperature investigated here. Our calculated results indicate that diffusion of CO 2 is faster in magnesite than in silicate melt. The results of CO 2 diffusion

Data Availability
All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest
The authors declare that they have no conflicts of interest.