Carbon Storage in Earth's Deep Interior Implied by Carbonate‐Silicate‐Iron Melt Miscibility

Carbonate melts have been proposed to exist in the lower mantle, but their interaction with other lower mantle melt compositions is poorly understood. To understand miscibility in the carbonate‐silicate‐metal melt system, we simulate endmember, binary, and ternary melt mixtures and study how their Gibbs free energies of mixing evolve with pressure. We find that carbonate‐metal and carbonate‐silicate melts have miscibility gaps that close with increasing pressure, while silicate‐metal melts are immiscible at all lower‐mantle pressures. Extending this analysis to the core‐mantle boundary, we suggest three miscible melt fields near the endmember carbonate, silicate, and iron melt compositions. Analysis of the densities of these miscible melt compositions indicates that some carbonate‐rich and some silicate‐rich melt compositions are gravitationally stable at the core‐mantle boundary and could be candidate compositions to explain ultra‐low velocity zones. Additionally, we evaluate the speciation of an example immiscible melt composition at various pressures throughout the mantle and identify reduced carbon species that would be expected to form in the melt. Our analysis reveals that a majority of Earth's carbon could have been transported to the core during core‐mantle differentiation and that much of Earth's carbon may be stored in the deep interior today.


Methods
Ab initio molecular dynamics simulations using the projector-augmented wave method (Kresse & Furthmuller, 1996) of density functional theory were performed with the Vienna ab initio simulation package (Blochl, 1994). We used the generalized gradient approximation in the Perdew-Burke-Ernzerhof form (Perdew et al., 1996) to treat electron exchange and correlation. The kinetic energy cutoffs for the plane-wave expansion of the wavefunctions were set to 600 eV. We used the canonical ensemble (NVT) with a Nosé-Hoover thermostat (Hoover, 1985;Nosé, 1984) with a time step of 1-2 fs for 18-80 ps, depending on the density. The Brillouin zone was sampled at the gamma point. The calculations were spin-polarized at all pressures. A Hubbard U eff (U − J) parameter of 4 eV was applied, which enhances the magnetic moment of the Fe atoms and corrects for their volume and coordination environment. The mean-square displacement as a function of time shows a ballistic regime below approximately 1,000 fs, after which the atoms reach a diffusive regime. For the carbonate-silicate-metal melt composition, calculations were run with a minimum of two starting configurations, and the results were averaged. We employ the Universal Molecular Dynamics package for the analysis of the results .
We work with seven melts representing endmember (MgCO 3 , MgSiO 3 , and Fe), binary (Mg(C,Si)O 3 , MgCO 3 + Fe, MgSiO 3 + Fe), and ternary (Mg(C,Si)O 3 + Fe) melt compositions (Table 1), with supercells ranging from 108 to 133 atoms. Simulations span a pressure range of 0-200 GPa, and all calculations are performed at 4000 K. Bond distances were determined from the pair distribution functions. The first peak in the pair distribution function marks the radius of the first coordination sphere for the reference atom, and the first minimum translates to the maximum acceptable bond distance for a bonding pair. The fitted minimum values were used in the speciation analysis to determine carbon clusters.

Melt Miscibilities
The Gibbs free energy of mixing, ΔG mix , determines whether a given solution of melt components will mix or unmix. Negative ΔG mix values indicate that a mixture is energetically favorable and therefore miscible. Positive ΔG mix values indicate that a mixture is energetically unfavorable and therefore immiscible. The Gibbs free energy of mixing was estimated using the following equation: Note. The numbers refer to the number of atoms included in the simulation. where ΔH mix is the enthalpy of mixing, T is the temperature, ΔS mix is the entropy of mixing, P is the pressure, and ΔV mix is the mixing volume. An example of the contribution of each term to ΔG mix on the carbonate-silicate binary is shown in Figure S1 of the Supporting Information S1. We describe the calculation of each term in the equation in the following sections.

Enthalpy of Mixing
To calculate the enthalpy of mixing, we use the equation: where β ij represents the binary parameter along the i − j binary, and X i and X j represent the mole fractions of the i and j components, respectively. Along each of the three binaries calculated in this study, X i and X j are known, but β is unknown and is dependent on the chemistry of the mixture. ΔH mix values for these compositions are hard to find, making it difficult to select the correct value of β along each binary. Solubility numbers, however, are easier to find in the literature, and can be directly compared to numbers extracted from our calculated ΔG mix curves (see Section 3.4). Thus, we employ an iterative guess and check method to determine which value of β best matches the solubility along a given binary under specific pressure and temperature conditions that can be compared to data from the literature or from our simulations. The selected β value is then used to calculate ΔH mix (Equation 2), and this ΔH mix value is used in all further calculations of ΔG mix (Equation 1) under different pressure and temperature conditions. To find appropriate values for β, we plot ΔG mix at 0 GPa and 4000 K along each binary and select β values that match the expected degree of mixing ( Figure S2 in Supporting Information S1). At 0 GPa, the mixing volume term of the ΔG mix equation goes to zero, and the entropy term at 4000 K is the same along all three compositional binaries studied (see Section 3.2). Thus, at 0 GPa and 4000 K, ΔG mix across all three binaries is influenced exclusively by the enthalpy term, and by changing the β value of ΔH mix , we directly change the degree of mixing along a given binary. Reasonable degrees of mixing are determined by examining solubility experiments on binary systems in addition to our own simulation results at 0 GPa. A study of orthopyroxene solubility in carbonate melts reports that carbonate melts contain 4 atomic percent silicate at 2 GPa and 1273 K (Kamenetsky & Yaxley, 2015). We select a value of 95 kJ for β on the carbonate-silicate join, which leads to silicate solubility of 8 atomic percent in carbonate melts at 4000 K. Silicate-metal melts are immiscible at 0 GPa (Fichtner et al., 2021), and our simulations show groupings of silicon and iron atoms that is suggestive of immiscibility. Thus, we select a β parameter of 135 kJ, which leads to limited miscibility (2 atomic percent Fe in silicate melt) at 0 GPa and 4000 K. Experimental reports of carbonate solubility in iron melt are lacking, but our carbonate-metal simulation indicates less miscibility than the carbonate-silicate simulation and more miscibility than the silicate-metal simulation at 0 GPa. We select a value of 115 kJ for the β parameter, which leads to 4 atomic percent Fe in the carbonate melt. ΔG mix values calculated with the chosen β parameters at 0 GPa are plotted in Figure S4 of the Supporting Information S1.

Entropy of Mixing
To calculate the entropy of mixing, we use the ideal entropy of mixing where R is the gas constant and X i , X j , and X k are the mole fractions of the i, j, and k components, respectively.

Mixing Volumes
Mixing volumes were calculated by taking weighted averages of molar volumes of individual melt components.
To calculate molar volumes, we first fit either second-or third-order Birch-Murnaghan equations of state for our simulated melt compositions ( Figure 1a). We also included the pyrolite and pyrolite + 8CO compositions from Solomatova et al. (2019) for comparison. The equations of state fit parameters are reported in Table S8 of the Supporting Information S1. The fit parameters reveal that melts with a carbonate component are highly compressible, which is in agreement with previous studies of carbon-bearing melts Ghosh et al., 2007;Sakamaki et al., 2011). There is a significant covariance between K 0 ', K 0 , and V 0 values for all melts, which is common in finite strain equations of state. One continuous equation of state was fit across multiple structural transitions, which stem from gradual coordination changes in the melt. However, the majority of coordination changes occur between 0 and 20 GPa, below the pressure regime of interest. In the pressure regime of the lower mantle, the fits closely match the data. Using the Birch-Murnaghan equation of state fits, we calculate molar volumes of the melts at pressures from 0 to 200 GPa ( Figure 1b). Due to the non-stoichiometric nature of the melt mixtures, we calculate volumes per mole of atoms instead of per formula unit, allowing the molar volumes to be directly compared. Iron and iron-bearing melts have the largest molar volumes, while pyrolite melts have the smallest. The densities and molar volumes of each melt composition are reported in Tables S1-S7 of the Supporting Information S1.
To calculate the mixing volumes, we compared the molar volume of our simulated mixture with the weighted average of the molar volumes of the mixture components. The magnitude of the mixing volume indicates the nonideality of a melt mixture and the degree of interaction between melt components. As the pressure derivative of ΔG mix , the mixing volume is the tendency of the mixture to become more or less energetically favorable with changing pressure. Thus, the sign of the mixing volume is suggestive of miscibility in multicomponent mixtures. This is especially true at high pressures, where the mixing volume term dominates the contribution to ΔG mix . For instance, a composition that is 50% MgSiO 3 and 50% Fe has a ΔG mix value of 58 kJ at 136 GPa and 4000 K (see Figure 3). Of the 58 kJ, 47 kJ is from the mixing volume component, accounting for 82% of the contribution to ΔG mix . Mixtures with positive mixing volumes become larger upon mixing and become less stable with increasing pressure, enforcing immiscibility. The mixing volumes for the four multicomponent melts in this study were plotted as a function of pressure in Figure 1c. To expand our analysis to any composition in the MgCO 3 −MgSiO 3 −Fe ternary system, the mixing volume data are fit to the following power series multicomponent mixing model (Ganguly, 2001;Wohl, 1946Wohl, , 1953: where ΔV mix is the mixing volume, the W G 's are the binary interaction parameters, and C ijk is the ternary interaction term. X i , X j , and X k are the mole fractions of the i, j, and k components, and X ji and X ij are the projected mole fractions of the i and j components in the binary join i − j. X ij is given analytically by  all at the temperature of 4000 K, calculated using the fit parameters reported in Table S9 of the Supporting Information S1. Along the carbonate-silicate join, mixing volumes are negative at all lower mantle conditions and become more negative with increasing pressure, suggesting continuous solubility. Along the silicate-metal join, mixing volumes are positive across all lower mantle pressure conditions, and decrease with pressure. At pressures beyond those of the Earth's mantle, we would expect silicate and metal melt to become miscible. Along the carbonate-metal join, mixing volumes are always negative but become less negative with increasing pressure, suggesting potential immiscibility beyond the core-mantle boundary pressure. Additionally, the magnitude of the mixing volumes represents the degree of interaction between the melts. Generally, silicate-metal melts have the most interaction, followed by carbonate-metal melts, and carbonate-silicate melts. Carbonate-silicate melt interaction terms are very small, even at their most negative point at 136 GPa, indicating that this mixture is close to ideal. Figures 2d-2f show the calculated mixing volumes for ternary compositions. Note that the mixing model was fit using four data points (three binary compositions and one ternary composition), which are labeled in Figures 2d-2f. At 24 GPa, melts with greater than 50% carbonate have negative mixing volumes. For melts less than 50% carbonate, mixing volumes are more negative with increasing iron percentage and more positive with increasing silicate percentage. With increasing pressure, the negative mixing volume regime shrinks and the positive mixing volume regime grows to cover more iron-and carbonate-rich parts of the ternary plot. By 136 GPa, only compositions that are greater than 70% carbonate and compositions close to the carbonate-metal and carbonate-silicate binaries have negative mixing volumes. Additionally, mixing volume magnitudes decrease with increasing pressure, indicating that these melts tend to become more ideal with increasing pressure. This conclusion is supported by the trends in the binary and ternary interaction parameters ( Figure S3 and Table S9 in Supporting Information S1). With increasing pressure, the interaction parameters trend toward 0, indicating that the interaction between melt components becomes increasingly less important with depth.

Gibbs Free Energy of Mixing
With the equations and approximations describing ΔH mix , ΔS mix , and ΔV mix , we determine how ΔG mix evolves along binary and ternary joins. ΔG mix is plotted along the carbonate-silicate, carbonate-metal, and silicate-metal binaries in Figures 3a-3c. MgCO 3 and MgSiO 3 melts demonstrate limited miscibility at all lower mantle pressures, with a large miscibility gap. Note that although ΔG mix is negative for many mixture compositions plotted along the MgCO 3 − MgSiO 3 binary, not all mixture compositions are miscible (see Figure 4 for the translation of the Gibbs free energy curve to a binary phase diagram). In the immiscible region, two melt compositions coexist, and these compositions are determined by the common tangent of the ΔG mix curves. These tangents are quasi-horizontal and have support points that are very close to the minima of the free energy. The shared tangents of the curves are plotted in Figure 4 and are at ∼9 and ∼91 mol percent MgCO 3 at all pressures examined for the carbonate-silicate binary. ΔG mix decreases with pressure, suggesting eventual closing of the miscibility gap at higher pressures than those reached by the Earth's mantle. MgSiO 3 and Fe are immiscible at all lower mantle pressures. The metallic character of the pure Fe melt makes it incompatible with the insulating character of the molten silicate melts, and any iron that is dissolved in the silicate is always incorporated as an ionic phase, FeO or Fe 2 O 3 . Similar to the MgSiO 3 and MgCO 3 binary, the MgCO 3 and Fe binary also has a miscibility gap that begins to close with increasing pressure. At 24 GPa, the two coexisting melt compositions are at 6 and 94 mol percent MgCO 3 , but by 136 GPa, the two coexisting melt compositions are at 38 and 62 MgCO 3 mole percentage (Figure 4). The ternary diagrams (Figures 3d-3f) show a range of miscibilities that expand with increasing pressure. As in the case of the binary phase diagrams, the miscibility fields for the ternary mixture compositions are set by the shared tangent in the ΔG mix curves, and not solely by where ΔG mix is calculated to be negative. These miscible melt fields, outlined by solid gold lines, are extended from ΔG mix values along each of the binaries. There are three miscibility fields, and each is located near an endmember composition. As the carbonate-metal and the carbonate-silicate miscibility gaps close with pressure, the miscibility fields grow to accommodate more mixing. In between the miscibility fields are two-melt regions, and compositions that fall in these regions will exsolve two immiscible melts. Dashed gold lines represent example tie lines in these regions. The central triangle is the three-melt region, and compositions that fall in this region will exsolve three immiscible melt compositions. With increasing pressure, the two-melt regions grow and the three-melt region shrinks, indicating the overall increase of miscibility in this system at high pressure. It is important to note that the chosen value for ΔH mix affects the miscibilities of the melt mixtures. We estimate values for ΔH mix based on chemical speciation and literature experiments, but without additional constraints, there is some ambiguity in the selected value. The blue regions of the plot indicate compositions with negative ΔG mix values, and thus, show a possible range of miscible compositions that are available under smaller ΔH mix values. Nonetheless, the results based on these reasonable estimates of ΔH mix illustrate the plausibility of reduced immiscibility with increasing pressure in the carbonate-silicate-metal system, such that an Fe-rich carbonate melt and a carbon-rich Fe melt would be expected to segregate from other phases at the base of the mantle.

Carbon-Bearing Clusters
As evidenced by Figures 3d-3f, many carbonate-silicate-iron melt compositions are immiscible at lower mantle conditions, even at the high-pressure conditions of the core-mantle boundary. Although ab initio molecular dynamics cannot model phase separation because of size effects, the clustering of species observed in our simulations is suggestive of the process of melt separation. In this section, we identify the species that segregate in an example ternary melt mixture through speciation analysis and determine individual cluster densities to understand how elements distribute through the lower mantle.
The example melt composition is the ternary melt composition simulated for this study (10% Fe, 45% MgCO 3 , 45% MgSiO 3 ). A more complete speciation analysis of this melt composition is reported by Davis et al. (2022). In general, large increases in C-Fe and C-C bonding with pressure at the expense of C-O bonding are observed. No evidence of Fe-Si bonding is found at any pressure. These tendencies are in good agreement with the predicted miscibilities for binary solutions (Figures 3a-3c). In the melt, large carbon-carbon clusters and, at higher pressures, carbon-iron clusters form and have limited interaction with the silicate melt network, indicating the types of carbon-bearing melt species we might expect to segregate from a silicate melt. However, the extent of carbon-iron interaction is difficult to quantify. Carbon-iron clusters often consist of iron atoms surrounding a polymerized carbon core, which could be classified as either a diamond seed nucleus or an iron carbide cluster. Previous simulations (Davis et al., 2022;Karki et al., 2020;Solomatova et al., 2019) show that we would expect carbon to bond to O, Fe, C, and Si. For each carbon atom, we would expect 12% (13/108) of the bonds to be C-Fe bonds as there are 13 iron atoms from available 108 coordinating anions for carbon (12 silicon, 72 oxygen, 11 out of the 12 carbon, and 13 iron). Additionally, we anticipate that 67% (72/108) of carbon bonds are to oxygen. Starting from these estimates based on the statistical sampling, we classify the carbon-based clusters in the melt. Bond abundances greater than the abundances predicted from statistical sampling indicate that there is a chemical preference for the bonding element. In our cluster analysis, we classify carbon clusters with C-O bond abundances greater than 67% as carbonates and carbon clusters with C-Fe bond abundances greater than 12% as carbides.
Clusters with less C-O and C-Fe abundances than expected from statistical sampling are classified as carbon polymers. The abundances of the different types of clusters at 1, 74, and 148 GPa are plotted in Figure 5d. At all three pressures, clusters of each type are formed, but the relative abundances of the cluster types evolve with pressure. Carbonates are the most abundant cluster type at 1 GPa and account for 49% of the total clusters, but that number drops to 31% at 74 GPa and 30% at 148 GPa. Carbides almost match the number of carbonate clusters at 1 GPa, at 48% of the total, and are the most abundant cluster type at 74 and 148 GPa, at 55% and 64% of the total, respectively. Polymers are always the least abundant cluster type. They increase in abundance from 3% to 14% from 1 to 74 GPa, and then decrease in abundance to 5% at 148 GPa. In Davis et al. (2022), we noticed that the majority of changes in C-O and C-C bond abundances occur in the first 25 GPa. Thus, we expect diamond formation to peak around 25 GPa, and this expectation is reflected in the relative increase in polymer formation between 1 and 74 GPa. Similarly, we expect carbonate cluster abundance to decrease rapidly in the first 25 GPa before plateauing, and this result is also observed. Finally, the large and linearly increasing number of carbide clusters matches the speciation results in both Davis et al. (2022) and Solomatova et al. (2019), which report linear increases in C-Fe bond abundances with increasing pressure.
The composition and the volume of the carbon clusters determine their relative density within the mantle. Using the Bader charge analysis algorithm (Henkelman et al., 2006;Sanville et al., 2007;Tang et al., 2009;Yu & Trinkle, 2011), we calculate the volumes of individual atoms within carbon clusters to determine cluster densities. The densities of example carbon clusters isolated at 74 and 148 GPa are plotted in Figure 6. The selected clusters are grouped according to their classification as a carbonate, carbide, or polymer. We directly compare the density of the cluster to the calculated density of MgSiO 3 melt under the same conditions. Carbide clusters are much denser than MgSiO 3 melt, and with enough time and aggregation, we expect these clusters to segregate from the multicomponent melt and sink to the core. Similarly, polymers are slightly denser than MgSiO 3 melt. Carbonate clusters are lighter than MgSiO 3 melt, and we expect these clusters to be buoyant within the mantle.

Implications
From the miscibility analysis, we found three miscible melt compositional fields: carbonate-rich, silicate-rich, and iron-rich melts (Figures 3d-3f). We consider the densities of these melt compositions to determine their buoyancy in the lower mantle and to evaluate their implications for carbon distribution and sequestration in the lower mantle and core. Densities of liquids calculated from molecular dynamics methods have been shown to systematically deviate from experimental values depending on the approximation used for the exchange correlation functional (Zhang et al., 2013;Zhao et al., 2014). However, relative comparisons of density between calculated melts are useful, provided the same approximations are made. As an example, density differences from MgSiO 3 melt at the core-mantle boundary are plotted in Figure 7. Here, miscible melt compositions could be formed from a deep Earth carbonatite melt interacting with iron melt at the core-mantle boundary. Compositions are shaded in red, white, and blue to represent buoyant, neutrally buoyant, and dense compositions, respectively, as compared to the density of MgSiO 3 melt, which is used as a proxy for the lower mantle composition. Three groups of miscible compositions emerge. Compositions outlined in red are iron-rich and denser than MgSiO 3 . We expect these compositions to sink into the outer core, dragging carbon and silicon out of the mantle and enriching the outer core with light elements over time. Compositions outlined in blue are carbonate-rich and buoyant. These compositions are anticipated to rise through the mantle and return carbon to shallower depths. Finally, the green compositions are neutrally buoyant and thus gravitationally stable at the core-mantle boundary. These compositions could serve as possible contributors to ULVZs. ULVZs have multiple proposed explanations, including FeSi formed through core-mantle reactions (Mergner et al., 2021), hydrous phases such as (Al,Fe)OOH (Thompson et al., 2021), Fe-rich post perovskite (Garnero & McNamara, 2008), Fe-rich (Mg,Fe) O (Solomatova et al., 2016;Wicks et al., 2010), and patches of partial melt (Williams & Garnero, 1996). Partial melt is a likely explanation for ULVZs due to the 3:1 ratio of S-to-P wave velocity reduction (Garnero & McNamara, 2008;Williams & Garnero, 1996). Therefore, buoyantly neutral melt compositions, such as the carbonate-silicate-metal melt compositions calculated in this study, could serve as one possible explanation and contribute to the ULVZs.
Within the immiscible melt compositions, carbon, carbon-iron, and carbon-oxygen clusters form ( Figure 5). Given enough time and aggregation, we expect the carbon and carbon-iron clusters to exsolve from the melt, as has been previously suggested (Dasgupta & Hirschmann, 2010;Karki et al., 2020;Mysen et al., 2011;Stagno et al., 2013). In our example ternary melt composition, the majority of the clusters formed at 148 GPa are carbide (64%), and the propensity for carbon to bond with iron indicates the high siderophility of carbon under these thermodynamic conditions. Carbide clusters are denser than the surrounding mantle ( Figure 6). Thus, a significant amount of Earth's carbon contained in the lower mantle may bond with iron and sink to the core, matching previous ab initio predictions of carbon's fate under reduced conditions in the lower mantle (Karki et al., 2020;Rohrbach & Schmidt, 2011). This not only prevents carbon from being recycled back to the Earth's surface but also changes the evolution of the core composition. An increasingly carbon-rich core composition would evolve to have density, sound velocities, and electrical and thermal conductivity more similar to the carbon-rich alloys Fe 3 C and Fe 7 C 3 (Fiquet et al., 2009;Wood et al., 2013). Moreover, given the chemical preference of carbon to be bonded to iron, we propose that during core formation, iron droplets that segregate from the magma ocean and fall downwards would constitute strong attraction basins for carbon. In this way, the magma ocean would be leached of its carbon. After the Moon-forming impact, metal and silicate melt would be well-mixed, and siderophile elements, such as carbon, would be segregated with iron into the core, supporting the idea that carbon is a candidate element to explain the density deficit in the core (Prescher et al., 2015;Solomatova et al., 2019).
In addition to carbide clusters, we observed the formation of carbon polymers in our simulated ternary melt composition, which could be precursors for diamonds. Our simulations reveal a possible mechanism for diamond formation, where carbon polymers exsolve from a silicate melt. Previously, this formation mechanism was observed in oxygen-deficient carbon-bearing silicate melts , and the addition of iron in our simulations may actually increase carbon polymerization (Belonoshko et al., 2015). Our previous speciation analysis of this melt composition (Davis et al., 2022) indicates that carbon-carbon bond formation reaches a peak around 25 GPa, beyond which it plateaus. C-Fe bonding, however, increases linearly and with increasing depth, the percentage of polymers decreases as carbide clusters are preferentially formed. Thus, our analysis indicates a diamond formation zone around 25 GPa. This depth in the Earth matches reports of diamonds with a deep Earth origin, which are returned from either the transition zone or the top of the lower mantle (Smith et al., 2016;Stachel et al., 2005). Ultradeep diamonds are thought to crystallize from a carbon-rich melt composition (Haggerty, 1999), which matches the formation mechanism outlined by our speciation analysis. Additionally, this formation mechanism is thought to play a role in eclogitic (E-type) diamonds (Gao et al., 2017;Haggerty, 1999). E-type diamonds crystallize within an eclogitic source rock, and have been reported to form at depths ranging from 150 to >660 km (Kirkley et al., 1991;Meyer, 1985;Moore & Gurney, 1985;Shirey & Shigley, 2013). They tend to form at higher pressure and temperature conditions than peridotitic (P-type) diamonds (Deines et al., 1993(Deines et al., , 2000Shirey & Shigley, 2013). E-type diamonds tend to have a wide range of δ 13 C values and are a good match for the carbon isotopic values of diamonds returned from the transition zone (Shirey & Shigley, 2013). Thus, the crystallization of a melt composition similar to our studied carbonate-silicate-metal melt, which forms abundant carbon polymers around the depth of the transition zone, may be a formation mechanism for ultradeep E-type diamonds. Our cluster analysis ( Figure 5) also indicates that polymers are formed even at the core-mantle boundary, and if these polymers aggregate to form diamonds, these diamonds may be brought to the surface by deep mantle plumes. Diamonds with a lowermost mantle origin may be identified through compositional analysis of fluid inclusions. Diamonds containing silicate-poor metal-rich carbonate melt compositions that fall into the miscible melt regions indicated in Figure 3f would indicate a core-mantle boundary origin and would provide evidence for carbonate-silicate-metal melt reactions in the lowermost mantle.
Finally, we examine the carbon distribution at pressure and temperature conditions of the core-mantle boundary to provide some insight into possible carbon distributions between core and mantle phases. We examine an equimolar composition (i.e., 1/3 Fe, 1/3 MgCO 3 , and 1/3 MgSiO 3 ), which falls into the three-melt region at the center of the ternary plot ( Figure 3f). Of the three melts that exsolve from this composition, 31% is a carbonaterich melt, 37% is an iron-rich melt, and 32% is a silicate-rich melt, where the melt compositions that exsolve are given by the corners of the miscible melt fields in Figure 3f. From mass balance calculations, we determine that for this case, 90% of the carbon is contained in outer-core compositions (carbonate-rich and iron-rich melts) and 10% is contained in a potential ULVZ composition (silicate-rich melt). In fact, carbon is distributed to varying degrees between outer-core and ULVZ-type compositions for the majority of compositions in this ternary system, and buoyant carbonate-rich melt phases only form when there is less than ∼5% Fe in the system. At the core-mantle boundary where iron melt is abundant, we anticipate that carbon is preferentially stored in the lower-mantle and core phases, indicating that the ultimate fate of Earth's carbon may be storage in the deep interior.

Conclusions
Carbonates are important compounds in the crust and upper mantle and may play a role in the lower mantle as well. Carbonate melts in the deep Earth may react with silicates and metals, especially at the core-mantle boundary where these phases are abundant. The chemical and physical properties of the melts that form from these reactions have important consequences for the distribution and storage of carbon in the deep Earth. Ab initio molecular dynamics simulations of carbonate-silicate-iron melt compositions allow for the examination of melt miscibilities, densities, and speciation. We find that carbonate-silicate and carbonate-iron melts have miscibility gaps that close with increasing pressure, and that carbonate-iron melts have the highest affinity for mixing. Silicate-iron melts are immiscible at all lower mantle pressures. By expanding our analysis to the ternary carbonate-silicate-iron system, we find that three miscible melt fields exist near each of the endmember compositions (Fe-rich, carbonate-rich, and silicate-rich melts). Iron-rich melts are dense and sink into the core, providing a mechanism to enrich the outer core with light elements, such as carbon, oxygen, and silicon. Silicate-rich melts are neutrally buoyant and sit at the core-mantle boundary, providing one possible explanation for the existence of ULVZs. Carbonate-rich melts, depending on their iron content, may sink into the core, remain at the core-mantle boundary, or rise through the mantle. Thus, depending on the composition that forms through the reaction of carbonate, silicate, and iron phases, carbon may be stored in the deep Earth in core-or ULVZ-type compositions or may return to shallower depths. The majority of melt compositions have densities that classify them as core-or ULVZ-type compositions, indicating that the fate of carbon may be to be stored in the Earth's deep interior. Finally, the speciation of carbonate-silicate-iron melts indicates that carbon polymers, iron carbides, and carbonate clusters are formed in the melt, and the relative proportions of these clusters at various pressures indicate carbon's changing affinity for the other elements. Iron carbides, which are favorably formed at higher pressures, indicate carbon's increasingly siderophile nature with depth. Carbon polymers, when aggregated, could form diamonds and are abundant at transition zone pressures, indicating a propensity for diamond formation in and around the transition zone. The distribution of carbon throughout the Earth's interior is a complicated topic, affected by many thermodynamic variables, including pressure, temperature, composition, and oxygen fugacity. More experimental and computational studies of carbonate melts and their interactions with other phases at lower-mantle and especially core-mantle boundary conditions will help elucidate the role of carbon in the Earth's deep interior.