Carbon-bearing silicate melt at deep mantle conditions

Knowledge about the incorporation and role of carbon in silicate magmas is crucial for our understanding of the deep mantle processes. CO2 bearing silicate melting and its relevance in the upper mantle regime have been extensively explored. Here we report first-principles molecular dynamics simulations of MgSiO3 melt containing carbon in three distinct oxidation states - CO2, CO, and C at conditions relevant for the whole mantle. Our results show that at low pressures up to 15 GPa, the carbon dioxide speciation is dominated by molecular form and carbonate ions. At higher pressures, the dominant species are silicon-polyhedral bound carbonates, tetrahedral coordination, and polymerized di-carbonates. Our results also indicate that CO2 component remains soluble in the melt at high pressures and the solution is nearly ideal. However, the elemental carbon and CO components show clustering of carbon atoms in the melt at high pressures, hinting towards possible exsolution of carbon from silicate melt at reduced oxygen contents. Although carbon lowers the melt density, the effect is modest at high pressures. Hence, it is likely that silicate melt above and below the mantle transition zone, and atop the core-mantle boundary could efficiently sequester significant amounts of carbon without being gravitationally unstable.

Volatiles such as carbon and hydrogen are important elements that are often present in shallow mantle as CO 2 and H 2 O fluids, which help in driving volcanic eruptions 1,2 . In the mantle, carbon also occurs as carbonates [3][4][5] , and as diamond and metal carbides under reducing mantle conditions [6][7][8] . Unlike hydrogen, carbon has very low solubility in major mantle silicate minerals 9,10 . It is likely that most of the carbon gets partitioned into molten silicates [11][12][13][14][15] . Although the present-day mantle is mostly solid, partial melts as evidenced from anomalous geophysical observations exist in various mantle conditions such as the mantle transition zone and ultra-low velocity zones at the core mantle boundary [16][17][18] . In addition, in the early history of the Earth, most of the silicate mantle might have been molten, resulting in a large scale magma ocean 19 . The remnant of such magma ocean is likely to have been preserved in the lowermost part of the present-day mantle [17][18][19][20] . It is possible that silicate magmas could have actually served as one of the largest repositories of mantle volatiles including carbon throughout the earth's history. However, we know very little about the incorporation of carbon in silicate melts and how it influences the atomistic scale structure and thus physical properties of the melt at pressure-temperature conditions pertaining to the whole mantle.
Extensive studies on the speciation of carbon in silicate melts [21][22][23] , the solubility of carbon 24,25 , and the effect of carbon on the melt density 13,14,26 have been mostly confined to conditions relevant for the deep crust and upper mantle. Recent numerical simulations on carbon-bearing basaltic and kimberlitic melts are also confined to pressures up to 15 GPa 27,28 . Here we investigate the high-pressure behavior of carbon bearing silicate melt using first-principles molecular dynamics (FPMD) simulation method 29 . In order to explore the effects of oxygen content on the speciation of carbon, we consider three distinct scenarios of dissolved carbon in MgSiO 3 system as carbon dioxide (CO 2 ), carbon monoxide (CO), and elemental carbon (C). We also explore a wide range of carbon concentration in terms of 5.2, 16.1, and 30.5 wt.% CO 2 (and equivalent CO and elemental carbon proportions). The calculated results on the structure and thermodynamic properties of carbon-bearing silicate melts under conditions of the whole mantle allow us for the first time to check whether carbon in its oxidized and reduced forms can be effectively be sequestered in high-pressure silicate melts and whether these carbon-bearing melts are sufficiently dense with respect to the surrounding solid mantle to be gravitationally stable.
We characterize the speciation of carbon dioxide in MgSiO 3 liquid in terms of coordination environments of carbon and oxygen (Fig. 1). The C-O coordination mainly consists of two-fold species (a molecular CO 2 ) and three-fold species (a carbonate ion, − CO 3 2 ). Their respective abundances at zero pressure and 2200 K are 5.7 and 94%. Both molecular and carbonate species have been experimentally detected in silicate melts and glasses below 10 GPa 21-23 . At higher pressures, the speciation changes: At the expense of CO 2 molecules, the abundance of carbonate ions increases, reaching as high as 98% at 5 GPa and 2200 K. On further compression, the − CO 3 2 population declines and is replaced by tetrahedrally coordinated species (CO 4 ) also found in solid carbonate minerals 4,5 . The CO 4 species becomes dominant at pressures above 60 GPa. Strong C-O short-range order is manifested in the corresponding radial distribution function (RDF) with a well-defined peak ( Supplementary Fig. S1). The CO 2 molecules and carbonate ions exist mostly in isolation at low pressures (Fig. 2, top) because oxygen rarely is coordinated with two carbons. However, our simulations for highly compressed system indicate the presence of polymerized di-carbonate complexes, C 2 O n (Fig. 2, bottom). Such di-carbonates have been also predicted previously in basaltic and kimberlitic melts 28 . The polymerized complexes are known to form when carbon dioxide dissolves in a carbonate melt via 30 Simulations show other types of species including C 2 O 6 and tetrahedrally coordinated di-carbonates, C 2 O 7 . The presence of -C-C-dimers at high pressure is manifested in the C-C RDF with new peak appearance at ~1.5 Å ( Supplementary Fig. S1).
The Si-O RDF and mean coordination remain largely unaffected by the carbon (Supplementary Figs S1 and S2). However, the Mg-O RDF and mean coordination are somewhat enhanced as a function of carbon-dioxide content in the silicate melt. Our results show that the conversion of molecular CO 2 → − CO 3 2 → CO 4 → polymerized di-carbonates upon compression involves increased sharing of oxygen with both magnesium and silicon. Almost all CO 2 molecules and majority carbonate ions (about 65% at 2200 K and 0 GPa) exist as free units at low pressure, i.e., these species are majorly bound to magnesium via Mg-O-C connections (Fig. 2, top). This would imply that as the mantle melt upwells through the upper mantle, the dominant − CO 3 2 species might serve as nucleating seeds for the formation of MgCO 3 from a carbonated silicate melt. Upon segregation of volatiles in the form of MgCO 3 , the remaining melt composed largely of networked SiO 2 units would freeze. Hence, faster removal of carbonated species and precipitation as MgCO 3 may account for elevating the solidus temperature to a volatile free silicate melting of >2000 K compared to that of 600 K for magnesite. At high pressure, carbon dioxide is much more integrated in the melt with mean C-O coordination exceeding 3.5 ( Supplementary Fig. S2). Pressure enhances the formation of Si-O-C connections as shown by the presence of the Si-C RDF peak ( Supplementary Fig. S1) and over 90% carbons form such polyhedral connections at pressures above 50 GPa. Upon compression, the nature of bonding in carbon species appears to evolve from ionic to covalent. The polyhedral association occurs mostly via the formation of bonds between carbon and non-bridging oxygen. Besides, there exist bridging CO n groups at elevated pressures arising from the formation of bonds between carbons and bridging oxygens.
Our results on elemental carbon (C) and carbon monoxide (CO) in MgSiO 3 liquid help us understand how the carbon speciation in silicate melt is affected under oxygen-poor, i.e., reducing conditions. We find that carbon incorporation occurs via both C-O and C-C bonding. In reduced silicate melts, the C-O coordination is severely suppressed whereas the C-C coordination is enhanced relative to the more oxidized CO 2 bearing silicate melt (Fig. 3). For the most reduced condition with elemental carbon in the amount corresponding to 16.1 wt.% CO 2 , about 50% carbons are bonded with oxygen. More than 80% carbon atoms are coordinated with each other forming carbon-carbon clusters, that is, the C-C bonds prevail in the melt. This clustering is less prominent at lower carbon content but it is enhanced with pressure. In contrast, 90% or more C atoms are not connected with each other in the oxidized case at pressures below 20 GPa so the presence of polymerized species such as di-carbonates is severely suppressed. The CO component shows an intermediate level of C-C bonding and clustering. The high C-C bonding activity of elemental carbon-bearing MgSiO 3 was previously predicted in the C-O-H-Fe system under conditions of reduced activity of oxygen 31 . Hence, our results demonstrate that under reducing conditions expected in the present day lower mantle, the carbon atoms may form clusters in silicate melt or may tend to exsolve out from the silicate melt and possibly form separate metal carbide melt at the core-mantle boundary conditions as recently suggested 32 .
The equation of state of the carbon bearing MgSiO 3 liquid as in the case of pure liquid 33 can be accurately described with P(ρ, T) = P(ρ, T 0 ) + B TH (ρ)(T − T 0 ), where P(ρ, T) is the total pressure at temperature T and density ρ. The reference isotherm P(ρ, T 0 ) at T 0 = 3000 K is the third order Birch-Murnaghan equation for each composition, with the fit parameters being somewhat sensitive to CO 2 concentration. The derived equation of state agrees well with experimental results 13,14 (Table 1). The thermal pressure coefficient B TH defined as (dP/dT) V increases by a factor of ten over the compression range explored in this study and is largely insensitive to composition ( Supplementary Fig. S3). All cases of carbon-bearing MgSiO 3 liquids are highly compressible but always remain less dense than the pure (carbon free) liquid (Fig. 4). In the low-pressure regime, the CO 2 component decreases the melt density by large amounts (Fig. 4, Inset). However, the density contrasts between the carbon bearing and pure melts decreases considerably at higher pressures. The density effects of carbon in the reduced forms of C and CO are similar to that of CO 2 for the same amount of carbon (Fig. 4, Supplementary Fig. S4).
The calculated density reductions (Δρ) of 0.015 to 0.005 g cm −3 per wt.% CO 2 between 10 and 140 GPa are smaller than those due to the water component (∼0.03 g cm −3 per wt.% H 2 O 34 ). Carbon-induced partial melting has been proposed atop mantle transition zone and in the uppermost parts of the lower mantle either owing to the changes in the oxidation state 6, 7, 15 and/or due to the lowering of solidus temperatures (in the presence of minor amounts of alkalies, for instance) so that the solidus intersects the mantle adiabat 35,36 . Based on the carbon-associated density differences found by this study, carbon-bearing silicate melts provide us with a viable mechanism to explain geophysical anomalies of the deep mantle 16,18 . Considering a peridotitic melt density of 3.66 gcm −3 at 13.6 GPa and 1873 K 37 and the mantle density of 3.54 gcm −3 at 410 km depth 38 , we find that with Δρ = 0.015 gcm −3 per wt.% CO 2 , the addition of up to 8 wt.% of CO 2 can produce a neutrally buoyant silicate melt at the top of the mantle transition zone. To assess the buoyancy situation at greater depths, we estimate the density of carbon-bearing MgSiO 3 liquid by including iron contribution (Supplementary text S1) for which we use the iron partitioning coefficient of ≤0.4 between the bulk solid mantle and partial melt 39 , and the bulk iron content of of 0.1. The estimated density of the silicate melt containing 5 wt.% CO 2 at 23.5 GPa and 2200 K (Fig. 4) is greater than the mantle density of 4.0 gcm −3 at 660 km depth 38 . Similarly, the melt density at 4000 K exceeds the mantle density near the CMB (Fig. 4).
The calculated thermodynamic results allow us to explore the behavior of the silicate melt-CO 2 solution. The calculated value of partial molar volume V CO 2 for the melt containing 5.2 wt.% CO 2 is 33 ± 4 cm 3 /mol at zero pressure and 3000 K. As pressure increases, V CO 2 decreases rapidly initially and then gradually (Fig. 5). The partial molar volume V CO 2 is in good agreement with the experimental estimates for basaltic 13 , peridotitic 14 , and komatitic 26 melts at low pressures. Pressure systematically suppresses both compositional and thermal effects so the Similarly, the C-C speciation is divided into two groups: free C atoms (not coordinated/bonded with any other carbon), and the other C atoms which are bonded/coordinated with one or more carbon atoms. Also shown are the results for 5.2 and 30.5 wt.% equivalent cases at high pressures. values of V CO 2 for all cases converge to 10 cm 3 /mol above 100 GPa. Compared to the partial molar volume of the CO 2 component, the estimated partial molar volumes of the CO and C components in MgSiO 3 liquid are systematically smaller, converging to 7 and 3.5 cm 3 /mol, respectively, at high pressure (Fig. 5). We have also calculated the the molar volume (V CO 2 ) of pure carbon dioxide fluid as a function of pressure. Our results show that the partial molar volume of melt CO 2 is significantly smaller than the molar volume (V CO 2 ) of pure fluid at zero pressure. Two volumes approach each other rapidly as pressure increases (Fig. 5). The volume of the melt-CO 2 solution, ΔV = V CO 2 − V CO 2 , is large and negative at low pressure (Fig. 5, Inset). For all isotherms and concentrations, its magnitude decreases rapidly with pressure and becomes nearly zero, implying an ideal solution within the computational uncertainty. Predicted negative ΔV (=d(ΔG mix )/dP < 0) means that the Gibbs free energy of mixing decreases with pressure. The calculated ΔH indeed decreases from positive values at low pressure to negative values at high pressure (Fig. 5, Inset) for all temperatures and concentrations. The melt-CO 2 solution is thus energetically favorable over a wide pressure range for carbon contents.
Our study based on first-principles molecular dynamics simulations provides valuable insights into the energetics of CO 2 miscibility in the silicate melts. The CO 2 solubility remains high over most of the mantle pressure regime. However, under reducing conditions, our results indicate carbon-carbon clusters implying that carbon solubility in silicate melts may reach a limit. It is generally known that the oxygen fugacity decreases significantly with increasing depth and the lower mantle is likely to be strongly reducing 8,40 . The oxidation state in deeper parts probably is inhomogeneous, consisting of locally or regionally oxidized domains associated with deeply subducting slabs where CO 2 -bearing silicate melts could be stabilized 15,36 . An upwelling silicate mantle could undergo redox melting across the lower mantle and mantle transition zone discontinuities (i.e., 660 km depth) so the mantle melting itself could be oxidizing 7,8 .
Although subducted carbonate-rich oceanic crust could transfer carbon up to a depth of 700 km, deeper subduction of slabs containing reduced carbon as metal carbide and subsequent melting to form carbon-bearing metallic melts has been recently proposed 32,36 . This carbon-rich melt is also likely to replenish the carbon in the deep mantle and eventually interact with any deep-seated silicate melt. If the present day carbon budget of the mantle has its origin dating back to the early Earth core formation period when a metallic-core segregated from a silicate mantle, the process might have sequestered significant amounts of carbon in a residual dense basal magma ocean at the bottom of the mantle and in possibly melt pockets associated with local anomalies detected in other parts of the mantle. Irrespective of the origin, such carbon may remain dissolved in silicate melt near the CMB perhaps in the oxidized forms because more of disproportionated metallic iron from the melt could have descended into the core thereby implying increased oxidation level. However, if the deep mantle is indeed strongly reduced as generally assumed, the reduced forms of carbon may still prefer the silicate melt (perhaps  38 . The inset shows the changes in the silicate melt density due to the CO 2 component at 4000 K for three concentrations calculated as ∆ρ = (ρ mix − ρ pure )/c. Here, ρ mix and ρ pure are the densities of the carbonated and pure MgSiO 3 liquids, respectively, and c is the wt.% of CO 2 content.
being incorporated as carbon-carbon clusters and/or metal carbides 32 ) over the solid mantle. Once incorporated, the carbon in silicate melts is likely to be sequestered owing to greater thermodynamic and gravitational stability. That mantle contained in deep-seated dense melts and perhaps in early magma ocean should represent an important part of the global carbon cycle 11 .

Methods
First principles molecular dynamics simulations within the local density approximation and projector augmented wave potentials were performed using the VASP program 29 . We explored three distinct scenarios: (a) MgSiO 3 melts with dissolved carbon dioxide (with 5. wt.% CO 2 . The pure CO 2 system used 72 atoms. For silicate melts with dissolved C and CO, the same numbers of corresponding formula units (C or CO) as in the case of CO 2 were included in the selected simulations.
Many canonical NVT ensembles were generated to cover a volume range of V = 1.5V X to 0.5V X (where V X = 2067.1 Å 3 ) corresponding to the entire mantle pressure regime (0 to 140 GPa), and a temperature range of 2200 to 6000 K. A cutoff energy of 400 eV (with Gamma point Brillouin zone sampling) was used, requiring the Pulay stress correction of 3 to 8 GPa for the volume range considered. At each volume, the system was melted and thermalized at 8000 K, and then subsequently quenched down to desired lower temperatures. The run durations with a time step of 1 femtosecond ranged from 50 to 300 picoseconds, depending on V-T conditions. Unlike most previous simulations of liquids, no empirical corrections were applied. The liquid state at each condition was confirmed by examining the radial distribution functions (RDF's) and mean square displacements (MSD's). The radial distribution function for each atomic pair species show a well-defined first peak in most cases and even a second peak sometimes, then converges to the unity value implying strong short-range order but no long-range order ( Supplementary Fig. S1). The MSD plots show that all atomic species reached diffusive regime within the simulation durations ( Supplementary Fig. S5). The calculated results on pressure, energies and various structural parameters were well converged with respect to the system size and runtime. Further details can be found in previous publications 41,42 . The calculated partial molar volume (V CO 2 ) of the carbon dioxide component in MgSiO 3 liquid at 3000 K (green circles), 4000 K (black circles), and 5000 K (red circles) for 16.1 wt.% CO 2 and at 4000 K also for 5.2 and 30.5 wt.%. The partial molar volume is defined as V CO 2 = (V mix − V pure )/n, where V mix and V pure are the volumes of the carbonated and pure silicate liquids, respectively, and n is the number of CO 2 molecules added to the system. The calculated molar volume of pure CO 2 fluid is shown along the corresponding isotherms (curves). The experimental data (asterisks) are for basalt at 19.5 GPa and 2573 K 13 , peridotite at pressures up to 3.8 GPa and 2050 K 14 , and komatite at 5.5 GPa and 2123 K 26 , containing 2.3 to 5.5 wt.% CO 2 . Also calculated are the partial molar volumes of the CO (open triangles) and C (filled triangles) components in the melt at 3000 K for the amount of carbon corresponding to 16.1 wt.% CO 2 . The partial molar volume of H 2 O in the hydrous melt (filled grey circles) from previous calculations and the molar volume of pure water (dotted curve) at 3000 K are shown for comparison 34,42 . The inset shows the calculated volume of mixing: ∆V = V CO 2 − V CO 2 and enthalpy of mixing:ΔH=H CO 2 − H CO 2 at the corresponding conditions. Here, H CO 2 = (H mix − H pure ) is the enthalpy per CO 2 formula unit in the melt and H CO 2 is that of pure CO 2 fluid.