# Chemical speciation and structure of kimberlite melts from ab initio molecular dynamics simulations

Ana Anzulović, Anne H Davis, Razvan Caracas

► To cite this version:

Ana Anzulović, Anne H Davis, Razvan Caracas. Chemical speciation and structure of kimberlite melts from ab initio molecular dynamics simulations. Geochimica et Cosmochimica Acta, 2025, 410, pp.203-217. 10.1016/j.gca.2025.10.016. hal-05377421

HAL Id: hal-05377421

https://hal.science/hal-05377421v1

Submitted on 22 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Chemical speciation and structure of kimberlite melts from *ab initio* molecular dynamics simulations

Ana Anzulović<sup>a,†</sup>, *, Anne H. Davis<sup>a</sup>, Razvan Caracas<sup>a,b,c †</sup>

## ARTICLE INFO

Associate editor: Xiandong Liu

Dataset link: https://doi.org/10.5281/zenodo.16400603

**Keywords:**

Kimberlite melt

Carbon

Volatiles

Water

Polymerization

Vaporization

## ABSTRACT

Kimberlite magmas, known for their rapid ascent from the upper mantle, transport volatiles, deep-sourced xenocrysts, and xenoliths to the surface. Their ascent rates suggest highly mobile, low-viscosity fluids, with their velocities further accelerated by the exsolution of volatile-rich fluid and vapor phases. Despite their role in driving kimberlite ascent, the atomic-scale mechanisms by which melt composition and volatile content influence physical behavior remain elusive, in part due to the scarcity of direct experimental constraints under relevant conditions. Here, we use *ab initio* molecular dynamics to address this gap. We model kimberlite melts with varying amounts of CO₂ (0-11 wt%) and H₂O (0-9 wt%). Our simulations provide a quantitative basis for the distinct structural roles of the major cations. Silicon acts as a rigid network former, with the SiO₄ tetrahedron being exceptionally stable (average lifetime ~0.6 ps), while network-modifying cations (Mg, Ca, Fe) form highly transient coordination environments (lifetimes <0.1 ps), quantifying their role in promoting the melt's characteristic mobility. The dominant speciation change with pressure is the conversion of molecular CO₂ to carbonate ions (CO₃²⁻), a transition largely completed below 3 GPa. The carbonate ion is the most stable polyanionic complex in the entire system, with a lifetime exceeding 1.2 ps—more than double that of SiO₄. This stability is more strongly controlled by the melt's redox state than by its total carbon or water content. Hydrous kimberlite melts have the capacity to start degassing at significant mantle depths (>90 km). This early and sustained degassing generates a buoyant, carbon-bearing vapor phase that provides a continuous driving force for the rapid ascent of kimberlite magmas.

## 1. Introduction

Kimberlites are volatile-rich magmas that transport mantle-derived material to the Earth's surface. They are derived from low-degree partial melting of carbonated peridotite where the peridotite solidus meets the cratonic geotherm (Wyllie, 1980). Characterized by highly depolymerized, silica-undersaturated nature, they exhibit high CO₂ solubility (Bailey, 1980). This positions them as transitional melts between silicic and carbonatitic melts (Giuliani et al., 2023; Sparks, 2013). Their source region is characterized by temperatures around 1500 K and depths ranging from 150 km to 450 km, corresponding to pressures between 5 and 15 GPa (Priestley et al., 2006; Sparks, 2013).

Carbon and hydrogen are the main volatile species of interest when considering kimberlite speciation and structure. Carbon dioxide, as the main driver of bubble nucleation (Holloway, 1976), is responsible for rapid devolatilization and eruption of kimberlite melts. Water, however, is often an overlooked volatile species. xperimental studies show

that H₂O-rich systems outperform peridotite-CO₂ systems in producing primitive kimberlite magmas at high pressures (Brey et al., 2008), which demonstrates the importance of water in the genesis of parental kimberlite melts. The minimum content of primary water to begin producing the parental magma necessary is estimated to be ~0.5 wt% (Hao et al., 2024). Measured water contents in natural kimberlites range from about 3 to 12 wt%, indicating the importance of water as a volatile species in kimberlite melts (Giuliani et al., 2023).

Despite extensive literature on kimberlites (Mitchell, 2013; Sparks, 1978; Sparks et al., 2006; Sparks, 2013; Giuliani et al., 2023; Russell et al., 2012; Kopylova et al., 2007), many questions regarding the detailed structure of kimberlite melts remain. In particular, the mechanisms through which carbon and water alter the structure of the silicate network is unclear. Previous work has indicated that the speciation of carbon in kimberlite melts strongly depends on pressure,

<footnote><sup>a</sup> Centre for Planetary Habitability, University of Oslo, Sem Sælands vei 2A, Oslo, 0371, Norway</footnote>

<footnote><sup>b</sup> Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, Paris, F-75005, France</footnote>

<footnote><sup>c</sup> The Research Institute of the University of Bucharest (ICUB), University of Bucharest, Bucharest, Romania</footnote>

<footnote>* Corresponding author.
E-mail address: anaanz@uio.no (A. Anzulović).</footnote>
carbon content, and melt composition (Guillot and Sator, 2011; Moussallam et al., 2015). The balance between carbonate ions ($CO_3^{2-}$) and molecular $CO_2$ is affected by the degree of melt polymerization, the concentration of network-forming cations (Si and Al), and the abundance of non-bridging oxygen (NBO) (King and Holloway, 2002; Guillot and Sator, 2011). A higher concentration of network-forming cations enhances polymerization by forming interconnected tetrahedral $SiO_4$ units. As the silicon content rises, the number of bridging oxygens (BO) increases, leading to a more polymerized melt structure. Carbonate ions are preferentially associated with NBOs, meaning they are more common in mafic and ultramafic melts, whereas molecular $CO_2$ is positively correlated with silicon content and the degree of polymerization (Blank and Brooker, 1994; Guillot and Sator, 2011; Mitchell, 2013; Ni and Keppler, 2013; Morizet et al., 2015; Moussallam et al., 2016b). As pressure decreases during magma ascent, $CO_2$ becomes less soluble, leading to its exsolution. However, this process is highly sensitive to melt composition. Water content further modifies these dynamics, reducing $CO_2$ solubility and promoting $CO_2$ vapor formation at higher pressures, enabling rapid ascent from greater depths (Field and Scott Smith, 1999; Moussallam et al., 2016b). The timing and conditions under which bubble nucleation occurs are not well constrained.

Previous studies have advanced our understanding of how carbon and water affect mafic and ultramafic melts, using both experiments and simulations (Konschak and Keppler, 2014; Moussallam et al., 2015, 2016b,a; Vuilleumier et al., 2015; Morizet et al., 2017; Solomatova et al., 2020; Solomatova and Caracas, 2021; Karki et al., 2021; Davis et al., 2022, 2023). But none of these studies consistently addresses kimberlite melts across a wide pressure-temperature-volatile space. Hence, our study is the first investigation based on Ab Initio Molecular Dynamics (AIMD) simulations that explores not only the behavior of dry kimberlite melts within a comprehensive thermodynamic and thermochemical space, but also how carbon and water influence their properties. We explore the structural evolution of kimberlites at the atomic level, under pressure-temperature conditions characteristic of the lithosphere and upper mantle. In addition to determining how carbon and water are incorporated into and break up silicate melt networks, we also gain insight into carbon solubility and degassing behavior. We determine the depth of nano-bubble nucleation in both water-free and water-bearing melts and describe the effect of water on carbon exsolution.

## 2. Methods

### 2.1. Ab Initio Molecular Dynamics

We performed AIMD simulations using the Vienna Ab Initio Simulation Package (VASP) (Kresse and Furthmüller, 1996). The electronic density was calculated using the projector augmented wave method (Blöchl, 1994) within Density Functional Theory (Kohn and Sham, 1965). We employed the generalized gradient approximation in the Perdew-Burke-Ernzerhof (Perdew et al., 1996) formulation (PBE) for the exchange-correlation function. While PBE might induce a shift in density/pressure compared to experiment, it remains the standard functional for studying silicate melts, and its efficiency is vastly superior to that of other, more exact functionals, such as hybrids, whose huge computational costs cannot balance their gain in accuracy. The kinetic energy cutoff was set to 600 eV for the plane-wave expansion of the wave functions and to 800 eV for the augmentation charges. The electronic wavefunctions were sampled in reciprocal space at the Gamma point. These parameters ensured a convergence of the internal energy to better than 15 meV/atom.

The simulations were performed within the NVT canonical ensemble with constant volume, V, number of atoms, N, and temperature, T. The temperature was controlled using a Nosé-Hoover thermostat (Hoover, 1985; Nosé, 1984). Simulations were performed at temperatures of 2500 and 3000 K. Each simulation was thermalized for 1 ps before

beginning production runs, which lasted at least 7 picoseconds (ps), with most extending to 10 ps. Additionally, a few longer simulations (15–60 ps) were used as support for the analysis of the structural properties, for each composition along each isotherm. The time step of the simulations was 1 femtosecond (fs). To ensure that the systems were fully molten at all temperatures, we computed the mean squared displacements of the atoms as a function of time. The positive value of the slope at 2500 K and 3000 K as well as the distance traveled by the atoms confirms that our systems are entirely molten at these temperatures at all pressures below ~35 GPa (Figure S4).

Our compositions were modeled off the hypabyssal Jericho kimberlite (Table 1), which is considered to be similar to those of primitive kimberlite melts (Price et al., 2000; Kopylova et al., 2007). We used the reconstructed kimberlite melt composition prior to the assimilation of orthopyroxene as reported in Kopylova et al. (2007). To build a ‘dry’ kimberlite melt, devoid of volatiles, we removed the volatile and trace elements from the measured composition and renormalized the values to 100%. Our chemical composition was based on this renormalized measurement. We generated a 265-atom cubic supercell with stoichiometry $Mg_{52}Si_{33}Ca_{17}Fe_8Ti_2Al_1K_1O_{151}$. Our theoretical composition agrees with the measured composition to within 1 wt% or better. The chemical compositions of the Jericho kimberlites and our theoretical dry kimberlite are presented in Table 1, with all values expressed in weight percent.

Using our dry composition as a base, we added varying amounts of $H_2O$ and $CO_2$ to create four new compositions:

* Carbon-poor kimberlite melt (CPK) = $DK + 8 CO_2$ (5.6 wt% $CO_2$)

* Carbon-rich kimberlite melt (CRK) = $DK + 17 CO_2$ (11.3 wt% $CO_2$)

* Water-poor kimberlite melt (WPK) = $CRK + 8 H_2O$ (11.3 wt% $CO_2$, 2.2 wt% $H_2O$)

* Water-rich kimberlite melt (WRK) = $CRK + 38 H_2O$ (11.3 wt% $CO_2$, 9.3 wt% $H_2O$)

The carbon and water contents were based on the range of volatile contents exhibited by the Jericho kimberlites (Kopylova et al., 2007). Although 9 wt% water may be unstable in natural kimberlite melts at high pressures, we analyze these simulations with the intent to represent conditions where kimberlites may interact with groundwater.

### 2.2. Analysis of the liquid

We use the UMD/Magmatomix package (Caracas et al., 2021) to process the AIMD results. The atomic trajectories and the thermodynamic parameters from each snapshot are extracted from the OUTCAR files generated by VASP.

We infer the presence of chemical bonds from the radial distribution functions (RDFs), which describe the probability of encountering a specific atomic type as a function of distance from a central atom. The first maximum of the RDF marks the most common interatomic distance. The first minimum of the RDF curve represents the limit of the first coordination sphere. We use this limit as a geometrical criterion for bond determination. If two atoms lie at a distance smaller than the radius of their corresponding first coordination sphere, they are considered bonded. From the identification of bonding pairs, we construct a connectivity matrix, which identifies different clusters of atoms and coordination polyhedra. In systems with coexisting liquid-vapor, the connectivity matrix reveals a bimodal distribution of atomic cluster sizes, where large clusters are identified as belonging to the liquid phase and small clusters (less than 10 atoms) represent gas molecules.

The pressure as a function of density and temperature is directly computed by VASP. At all conditions, we add the thermal contribution stemming from the kinetic energy of the atoms to the value of the pressure. At low densities, the liquid may enter an extension regime, where calculated pressures are negative. This indicates the coexistence
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Average bond lengths as a function of pressure, derived from the first peak of the radial distribution function. The figure compares results at 2500 K (dashed lines) and 3000 K (solid lines) for four compositions: (a) dry kimberlite melt (DK), (b) water-rich kimberlite melt (WRK), (c) carbon-rich kimberlite melt (CRK), and (d) water-poor kimberlite melt (WPK). Bond lengths remain largely constant across the 45 GPa pressure range, highlighting their incompressible nature.</figcaption></figure>

**Table 1**

Oxide compositions based on measurements on real rocks (Kopylova et al., 2007) are in column 'Rock'. For each composition, in the first column, we list the renormalized oxide values from Kopylova et al. (2007), followed by the compositions from the kimberlite model and the calculated differences ('Diff') between them. The final column in each composition ('Cat') indicates the number of cations in the kimberlite composition.

<table><thead><tr><th colspan="2">Oxides</th><th colspan="4">Dry kimberlite</th><th colspan="4">Carbon-poor kimberlite</th></tr><tr><th>wt%</th><th>Rock</th><th>Renorm</th><th>Model</th><th>Diff</th><th>Cat</th><th>Renorm</th><th>Model</th><th>Diff</th><th>Cat</th></tr></thead><tbody><tr><td>SiO<sub>2</sub></td><td>26.70</td><td>33.38</td><td>33.25</td><td>-0.12</td><td>33</td><td>31.42</td><td>31.40</td><td>-0.02</td><td>33</td></tr><tr><td>TiO<sub>2</sub></td><td>1.73</td><td>2.16</td><td>2.68</td><td>0.52</td><td>2</td><td>2.04</td><td>2.53</td><td>0.50</td><td>2</td></tr><tr><td>Al<sub>2</sub>O<sub>3</sub></td><td>1.57</td><td>1.96</td><td>1.71</td><td>-0.25</td><td>1</td><td>1.85</td><td>1.61</td><td>-0.23</td><td>1</td></tr><tr><td>FeO</td><td>7.58</td><td>9.48</td><td>9.64</td><td>0.16</td><td>8</td><td>8.92</td><td>9.10</td><td>0.18</td><td>8</td></tr><tr><td>MgO</td><td>28.25</td><td>35.32</td><td>35.15</td><td>-0.17</td><td>52</td><td>33.24</td><td>33.19</td><td>-0.05</td><td>52</td></tr><tr><td>CaO</td><td>12.90</td><td>16.13</td><td>15.99</td><td>-0.14</td><td>17</td><td>15.18</td><td>15.10</td><td>-0.08</td><td>17</td></tr><tr><td>K<sub>2</sub>O</td><td>1.26</td><td>1.58</td><td>1.58</td><td>0.00</td><td>1</td><td>1.48</td><td>1.49</td><td>0.01</td><td>1</td></tr><tr><td>CO<sub>2</sub></td><td>9.88</td><td></td><td></td><td></td><td></td><td>5.88</td><td>5.58</td><td>-0.31</td><td>8</td></tr></tbody></table>

<table><thead><tr><th colspan="2">Oxides</th><th colspan="4">Carbon-rich kimberlite</th><th colspan="4">Water-poor kimberlite</th><th colspan="4">Water-rich kimberlite</th></tr><tr><th>wt%</th><th>Rock</th><th>Renorm</th><th>Model</th><th>Diff</th><th>Cat</th><th>Renorm</th><th>Model</th><th>Diff</th><th>Cat</th><th>Renorm</th><th>Model</th><th>Diff</th><th>Cat</th></tr></thead><tbody><tr><td>SiO<sub>2</sub></td><td>26.70</td><td>29.71</td><td>29.55</td><td>-0.16</td><td>33</td><td>29.06</td><td>28.93</td><td>-0.14</td><td>33</td><td>26.99</td><td>26.81</td><td>-0.17</td><td>33</td></tr><tr><td>TiO<sub>2</sub></td><td>1.73</td><td>1.93</td><td>2.38</td><td>0.46</td><td>2</td><td>1.88</td><td>2.33</td><td>0.45</td><td>2</td><td>1.75</td><td>2.16</td><td>0.41</td><td>2</td></tr><tr><td>Al<sub>2</sub>O<sub>3</sub></td><td>1.57</td><td>1.75</td><td>1.52</td><td>-0.23</td><td>1</td><td>1.71</td><td>1.49</td><td>-0.22</td><td>1</td><td>1.59</td><td>1.38</td><td>-0.21</td><td>1</td></tr><tr><td>FeO</td><td>7.58</td><td>8.43</td><td>8.56</td><td>0.13</td><td>8</td><td>8.25</td><td>8.38</td><td>0.13</td><td>8</td><td>7.66</td><td>7.77</td><td>0.11</td><td>8</td></tr><tr><td>MgO</td><td>28.25</td><td>31.43</td><td>31.23</td><td>-0.20</td><td>52</td><td>30.75</td><td>30.57</td><td>-0.18</td><td>52</td><td>28.55</td><td>28.34</td><td>-0.21</td><td>52</td></tr><tr><td>CaO</td><td>12.90</td><td>14.35</td><td>14.21</td><td>-0.15</td><td>17</td><td>14.04</td><td>13.91</td><td>-0.14</td><td>17</td><td>13.04</td><td>12.89</td><td>-0.15</td><td>17</td></tr><tr><td>K<sub>2</sub>O</td><td>1.26</td><td>1.40</td><td>1.40</td><td>0.00</td><td>1</td><td>1.37</td><td>1.37</td><td>0.00</td><td>1</td><td>1.27</td><td>1.27</td><td>0.00</td><td>1</td></tr><tr><td>CO<sub>2</sub></td><td>9.88</td><td>10.99</td><td>11.15</td><td>0.16</td><td>17</td><td>10.75</td><td>10.92</td><td>0.16</td><td>17</td><td>9.99</td><td>10.12</td><td>0.13</td><td>17</td></tr><tr><td>H<sub>2</sub>O</td><td>9.07</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0</td><td>2.18</td><td>2.10</td><td>-0.07</td><td>8</td><td>9.17</td><td>9.26</td><td>0.09</td><td>38</td></tr></tbody></table>

of liquid and vapor phases (Caracas et al., 2021; Caracas and Stewart, 2023; Solomatova et al., 2020; Davis and Caracas, 2024). In these special cases, we followed the procedure of Davis and Caracas (2024) and computed the pressure of the vapor phase using a modified Van der Waals equation of state (Van der Waals, 1873) of the form:

$$P = \left( \frac{N k_B T}{V} - \frac{a_{vdW} N^2}{N_A^2 V^2} \right)$$

where P is the pressure, V is the volume, N is the number of molecules, k<sub>B</sub> is Boltzmann's constant, T is the temperature, N<sub>A</sub> is Avogadro's number, and a<sub>vdW</sub> is a constant that corrects for intermolecular forces. The a<sub>vdW</sub> constants for the most common gas species in our simulations (CO, CO<sub>2</sub>, and O<sub>2</sub>) were taken from the Handbook of Chemistry and

Physics (Weast, 1972). The volume of the vapor phase was obtained by subtracting the volume of the liquid from the volume of the simulation box. To find the volume of the liquid phase, we identified the volumes of individual atoms using the Bader charge analysis algorithm (Tang et al., 2009; Sanville et al., 2007; Henkelman et al., 2006) with the default vacuum flag. Volumes were summed for atoms we identified as part of the liquid.

Our simulations, involving 300 to 400 atoms, are larger than those in most previous studies, which typically include around 100–150 atoms. However, this size might still introduce finite-size effects that would influence the onset and morphology of vapor–liquid phase separation. In particular, nanobubble nuclei may approach the box size, and periodic boundaries can distort interface dynamics. While this limits
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Average cation-oxygen coordination numbers as a function of pressure. Each row shows a different cation coordinated by oxygen: (a) Mg-O<sub>x</sub>, (b) Si-O<sub>x</sub>, (c) Ca-O<sub>x</sub> and (d) Fe-O<sub>x</sub> as a function of pressure. Each column represents a distinct kimberlite composition (abbreviations defined in the text). The figure highlights the stability of the Si-O network-forming backbone, whose coordination remains between 4 and 5 across the pressure range. In contrast, the network-modifying cations show a clear hierarchy of adaptability, with the coordination of the large Ca<sup>2+</sup> ion increasing most steeply, followed by the smaller, less flexible Mg<sup>2+</sup> and Fe<sup>2+</sup> ions.</figcaption></figure>

quantitative interpretation of nucleation pressures, the consistent emer-
gence of nanocavities across compositions and conditions supports the
trends we report. These results provide insight into the early degassing
behavior of volatile-rich melts, despite the inherent scale limitations.

3. Results

3.1. *Bond lengths*

The first maximum of the RDF curve is indicative of the average
bond length, and these are plotted in Fig. 1 for each bonding pair.
Notably, pressure does not have a large effect on the observed bond
distance; rather, increasing pressure is accommodated through more
efficient packing of atoms, resulting in increasing coordination numbers
for the cations (Solomatova and Caracas, 2019; Davis et al., 2022) (Fig.
2).

All cations form bonds to oxygen. K-O, Ca-O, and Mg-O bond lengths decrease as pressure increases from 0 to 40 GPa. The variability in K-O bond lengths reflects fluctuations in potassium coordination caused by limited sampling—due to only one K atom per simulation cell— and by transient changes in the immediate oxygen environment around potassium Ti-O and Fe-O bonds, however, slightly increase in length across the same pressure range due to changes in coordination state. Bond lengths for the network-forming cations, aluminum and silicon, as well as bond lengths for volatile-oxygen pairs (C-O and H-O) show minimal variation under pressure.

Beyond cation-oxygen bonding, we see evidence of C and H assum-
ing the role of oxygen and forming bonds with cations. In the hydrated

compositions, H bonds with C and H. We see H-H bonds forming only
at low pressures and at 3000 K (Fig. 1b), indicating a hydrogen gas
phase (Table T1). In some of the simulations, carbon displays bonding
behavior with Fe, C, and H. Many of these minor bonding events, occur
only intermittently and may result from transient atomic proximity due
to thermal motion rather than stable bonding. Given their inconsistency
and the high-temperature conditions, these are likely random collisions.

Previous work on carbon-bearing melt systems (Solomatova et al., 2020; Ghosh et al., 2017; Davis et al., 2022) indicates that C–C bonding is more abundant at higher pressures, lower oxygen fugacities, and larger carbon concentrations. While kimberlite melts tend to be present at lower pressures (0–10 GPa) and comparably oxidizing environments, their high carbon contents seemingly indicate an inclination towards C–C bonding. However, even in our longest simulation in a CPK composition at 3000 K and 1.1 GPa (60 ps), we do not observe consistent C–C bonding or clustering of carbon atoms. This can be explained by the polymerization of the kimberlite melt (see Section 3.3), which allows carbon to play a more important structural role compared to more depolymerized silicate melts. Carbon is incorporated as a cation into the kimberlite melt structure, forming part of the polymerized network that restricts its movement and inhibits its ability to find and bond with other carbon atoms. In the DK and CRK, we observe no significant difference in bond lengths when comparing them at 2500 K and 3000 K.

3.2. *Atomic coordination*

3.2.1. *Magnesium-oxygen*

At 3000 K, in all simulated kimberlite compositions, the average
coordination state of magnesium increases from 4.2 to 7 (Fig. 2a). Our
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Relative abundance of specific cation-oxygen coordination states as a function of pressure at 3000 K. This figure details the distribution of coordination states for (a) Mg, (b) Si, (c) Ca, and (d) Fe. Each row corresponds to a different kimberlite composition. This figure illustrates the distinct pressure-response mechanisms of the cations. The network-former Si exhibits a discrete structural transition from an overwhelmingly dominated (> 80%) four-fold (SiO₄) to a five-fold (SiO₅) state above 20 GPa. In stark contrast, the network-modifiers (Mg, Ca, Fe) show a broad, continuously shifting distribution of multiple coexisting coordination states, highlighting their flexible role in accommodating melt densification.</figcaption></figure>

results closely align with other studies (Sun et al., 2011; Solomatova and Caracas, 2019; Davis et al., 2022), despite differences in melt composition. Mid-Ocean-Ridge-Basalt (MORB) melts (Bajgain et al., 2015) tend to have higher average magnesium CNs above 5 GPa when compared to kimberlites, pyrolite (Solomatova and Caracas, 2019), CaMgSi₂SiO₆ (Sun et al., 2011), and Mg(C,Si)O₃ + Fe melts (Davis et al., 2022) (Fig. 2a). Comparing hydrous basalt (Bajgain et al., 2015) with hydrated kimberlite melts (WPK and WRK), we see that the average magnesium CN for kimberlite is slightly lower than the hydrous basalt (Fig. 2a).

Magnesium displays a smooth transition from lower to higher coordination states for all kimberlite compositions (Fig. 3a). For example, in all compositions at 3000 K, fourfold coordination is initially dominant, before fivefold coordination becomes the most abundant coordination state between 0 and 9 GPa. Magnesium transitions to primarily sixfold and sevenfold coordination at 9 and at 25 GPa, respectively. At no point is the most dominant coordination state more than 50% of all magnesium species for all compositions studied.

The lifetimes of different Mg coordination states do not show significant variation. In Fig. 4, we plot the lifetimes of MgOₓ coordination states for the WPK composition at 1.4 GPa and 3000 K. We find that fivefold-coordinated Mg species are the most abundant, exceeding 20 000 individual instances, and have an average lifetime of 43 fs. Fourfold-coordinated species are the next most abundant, with more than 14 000 occurrences and an average lifetime of 47 fs.

### 3.2.2. Silicon-oxygen

At 3000 K, in all simulated kimberlite compositions, the average CN of Si to O changes from 4 to 5.2 (Fig. 2b). The average CN of Si to O is higher in basalts (Bajgain et al., 2015; Sanloup et al., 2013) and CaMgSi₂SiO₆ (Sun et al., 2011) compared to kimberlite melts, both in the anhydrous and the hydrous case. The consistently lower Si-O CN in kimberlites distinguishes these ultramafic, carbonated melts from more silica-rich compositions, such as basalt. This structural difference can be attributed to the high concentration of network-modifying cations (e.g., Mg²⁺, Fe²⁺), which depolymerize the silicate framework and the presence of carbonate groups. CO₃ groups bond with the NBO of the silicate structure (Fig. 9). The planar CO₃ group, associated with the silicate framework, is sterically bulky. This creates sterically hindered, relatively rigid local structures. As pressure increases, these bulky carbonate groups act as ‘spacers’ within the melt, physically impeding the approach of additional oxygen atoms to silicon and thus slowing down the structural transformation that would result in higher Si-O coordination states.

The experimentally derived Si-O CNs reported by Sanloup et al. (2013) exhibit a sharper increase with pressure than observed in simulations, rising from approximately 4.5 at moderate pressure to 6.0 at higher pressures. This contrasts with first-principles molecular dynamics simulations on related compositions, such as molten anorthite and diopside (de Koker et al., 2008; Sun et al., 2011), which predict a more gradual coordination increase and yield values more consistent with our results. The divergence likely stems from methodological
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Lifetimes of individual cation-oxygen coordination polyhedra in the WPK composition at ~1.4 GPa and 3000 K. Each column represents a cation bonded to oxygen, and each row corresponds to a specific oxygen coordination number (indicated on the left). The x-axis shows the number of occurrences of each cation-oxygen cluster, sorted in decreasing order of lifetimes, while the y-axis represents species lifetimes in femtoseconds. Average lifetime (Avg) expressed in femtoseconds is shown in the top right corner of each panel. This figure provides a quantitative basis for the distinct roles of different ions in the melt. The network modifiers (Mg, Ca, Fe) form highly transient coordination environments with average lifetimes consistently below 100 fs. The network-forming SiO₄ tetrahedron is an order of magnitude more stable with average lifetime of 570 fs. The most significant finding, however, is the stability of the carbonate ion (CO₃), which, with an average lifetime over 1.2 ps, is the most robust structural unit in the system. The long lifetimes of CO₃ and other volatile species like OH (231 fs) confirm their role as integral, chemically-bonded components of the melt.</figcaption></figure>

differences: Sanloup et al. (2013) derive coordination from total structure factor fitting, which incorporates all atomic pair contributions and depends on challenging background subtractions and density estimations, potentially introducing uncertainty. Additionally, compositional complexity, temperature, and melt stiffening at high pressure may further contribute to the discrepancy. Interestingly, their results align more closely with simulations on pure silica, highlighting the influence of composition and experimental resolution on coordination measurements.

Pyrolite (Solomatova and Caracas, 2019) and Mg(C,Si)O₃ + Fe melts (Davis et al., 2022) have very similar average Si-O coordination states compared to the kimberlites.

Silicon adopts the fewest unique coordination states of the cations in this study. Over the studied pressure range, silicon primarily forms SiO₄ tetrahedra (Fig. 3b), which are the fundamental building blocks of silicate melts at crustal and subcrustal conditions. SiO₄ polyhedra exhibit the longest average lifetimes among all species, measuring 658 fs in the WPK composition (Fig. 4). In all of the analyzed compositions, more than 80% of Si is present as SiO₄ polyhedra from 0 to 9 GPa,
<figure><img src="image_7.png" /></figure>

**Fig. 5.** Pressure-dependent carbon speciation in kimberlites and related ultramafic melts. The main panels show the relative abundance of (a) carbonate (CO₃), (c) molecular CO₂, and (d) CO and CO₃ species. Panel (b) provides a detailed view of the CO coordination number distribution at low pressures (0–8 GPa), illustrating the transition from two-fold (CO₂) to three-fold (CO₃) coordination. Data from this study's kimberlite models (filled/unfilled symbols) are compared with literature data for other mantle melts: pyrolite+CO, pyrolite+CO₂ (Solomatova and Caracas, 2019) and Mg(C,Si)O₃ + Fe melt (Davis et al., 2022) (dashed lines). A comparison of the literature data in Panel (a) shows that as melts become more reduced, the abundance of CO₃ decreases significantly. This is compensated by an increase in molecular CO (d), as a greater proportion of carbon engages in C-C, C-Fe, and C-Si bonding rather than C-O bonding (non C-O bonds not shown). The abundance of CO₂ (c) and CO (d) decreases sharply with pressure, mirrored by a rapid increase in carbonate (CO₃) abundance (a). Panel (b) details this structural shift at low pressures, showing the crossover from two-fold to three-fold C-O coordination. This conversion is the dominant process across our kimberlite compositions. As melts get more reduced less CO₃ is present and aligns with trends from literature data (dashed lines). At pressures above 20 GPa, a secondary transition towards four-fold coordinated carbon (CO₄) begins.

and fourfold coordination remains the most common coordination state up to about 30 GPa. At low pressures, SiO₃ polyhedra are present, but their abundance decreases linearly to almost 0% above 12 GPa. SiO₅ polyhedra also exist in the melt in minute amounts at surface conditions. Their abundance increases linearly with pressure, and they become the dominant Si-O species above ~30 GPa. This transition from a SiO₄⁻ to a SiO₅-dominated system occurs at 27 GPa for both the DK composition and the CPK composition. When more carbon is present in the system (CRK), the transition occurs at 35 GPa. Across all compositions, SiO₆ octahedra start to appear in the melt around 10 GPa. Their abundance increases with pressure, and they should ultimately become the dominant species at lower mantle pressures, as is the case for other silicate melts (Solomatova and Caracas, 2019; Karki et al., 2021; de Koker et al., 2008).

We counted over 2000 SiO₄ tetrahedra in the WPK composition at 1.4 GPa and 3000 K (Fig. 4). They are long-lived, with average lifetime of 570 fs. There are over 1000 SiO₅ species, with an average lifetime of 53 fs, and fewer than 1000 SiO₃ clusters, with an average lifetime of 95 fs.

### 3.2.3. Calcium-oxygen

The average CN for Ca-O species is 5.35 under ambient pressure conditions and increases to 8 between 10 and 13 GPa, and between 9 and 10 at the highest pressure across all analyzed compositions (Fig. 2c). The average Ca-O coordination state is weakly sensitive to composition as its values are within 10% for the different liquids compared here (Bajgain et al., 2015; Solomatova and Caracas, 2019; Sun et al., 2011). As the second-largest cation in this study (after

potassium), calcium displays a particularly wide range of coordination states. Its large ionic radius and flexible bonding environment enable a continuous increase in coordination with pressure, accommodating melt densification. Consequently, Ca-O exhibits the steepest coordination-pressure slope among all species considered (Fig. 2c).

At 0 GPa, most Ca-O species are sixfold-coordinated, transitioning to sevenfold-coordinated at 2 GPa for all studied compositions (Fig. 3c). Interestingly, the inclusion of volatile species causes the transition to higher coordination states to occur at lower pressures. For instance, eightfold coordination becomes the most common coordination state at 9 GPa for the DK composition, 7.8 GPa for the CPK composition, 7.5 GPa for the CRK composition, 6 GPa for the WPK composition, and 7 GPa for the WRK composition.

There are over 5000 CaO₆ and CaO₅ species in the WPK composition at 1.4 GPa and 3000 K (Fig. 4). Ca-O species have short average lifetimes, of around 30 femtoseconds.

### 3.2.4. Iron-oxygen

The average CN of Fe to O varies from 3.5 under ambient pressure conditions to 4.5–5.2 at 9 GPa depending on the volatile content of the kimberlite melt (Fig. 2d). Our results align well with pyrolite (Solomatova and Caracas, 2019) and fall between MORB (Bajgain et al., 2015) and Mg(C,Si)O₃ + Fe melts (Davis et al., 2022). Note that the Mg(C,Si)O₃ + Fe numbers include non-oxygen bonding in the calculation of average CN and the numbers are elevated due to the presence of Fe-C bonding in the system, especially at higher pressures. Additionally, our calculations are spin-polarized, while the MORB calculations (Bajgain et al., 2015) are not, which may explain the discrepancy between values. In the DK composition, iron adopts tetrahedral coordination
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Carbon speciation and coordination in kimberlite melts. The pressure-driven transition from molecular CO₂ to carbonate is demonstrated by (a) the steep increase in the average C-O coordination number towards 3 and (c) the rapid decrease in the $\frac{N_{CO_2}}{N_{CO_2}+N_{CO_3}}$ ratio below 3 GPa. (b) Carbon-oxygen coordination state abundances across different temperatures for the carbon-rich kimberlite melt (CRK). Our results for this transition (c) are consistent with literature data for a range of ultramafic melts. Results from this study are shown with shaded regions, where empty symbols denote the CPK composition and full symbols denote the CRK composition. Values from the literature with carbon contents comparable to the CPK composition are shown with empty symbols: carbonated MORB at 2272 K (Guillot and Sator, 2011) (cyan), pyrolite + 5 wt% CO₂ (Solomatova and Caracas, 2019) (peach), and forsterite + 7 wt% CO (Solomatova et al., 2020) (orange). Enstatite + 16 wt% CO₂ (Ghosh et al., 2017) (teal) has carbon contents comparable to the CRK composition. The temperature dependence is shown in (b), where we see minor variation in carbon-oxygen abundances at 2500 K and 3000 K.</figcaption></figure>

from 0 to 9 GPa, transitions to fivefold coordination between 9 and 20 GPa, and adopts octahedral coordination above 20 GPa (Fig. 3d).

In the WPK composition at 1.4 GPa and 3000 K, there are more than 2000 FeO₄ species, more than a 1000 FeO₅ and FeO₃ species (Fig. 4). The average lifetimes of Fe-O species are between 26 and 62 femtoseconds.

### 3.2.5. Carbon-oxygen

The chemical environment of carbon in kimberlite melts is dominated by oxygen bonding. C-O bonds make up more than 96 percent of all carbon bonds, and this number stays consistent across all thermodynamic conditions and chemical compositions explored in this study.

The average C-O CN in kimberlite melts exhibits a significant dependence on the pressure (Fig. 6a). Between 0 and 3 GPa, the coordination increases sharply from 2.25 to 3, suggesting a rapid structural reorganization of the melt. This is likely due to the incorporation of gaseous CO₂ species into the melt as CO₃ groups. As pressures continue to rise, the CN increases further, ultimately reaching values between 3 and 3.3 at pressures exceeding 30 GPa. This indicates a gradual transition from sp²- to sp³-hybridized carbon, as previously witnessed in carbonate minerals and carbon-bearing melts at high pressure (Boulard et al., 2011; Yao et al., 2018; Sanloup et al., 2019; Davis et al., 2022).

In all compositions, C-O speciation is dominated by CO₃, starting between 1–3 GPa and continuing to at least ~45 GPa, the highest pressure in this study (Fig. 5a). This behavior is consistent with previous findings for oxidized carbon-bearing silicate melts, such as pyrolite+CO₂ (Solomatova and Caracas, 2019), where carbonate is the dominant species across upper mantle pressures. Carbon speciation is governed to a greater extent by the redox state than the CO₂ amount (CPK and CRK) or even than the presence of water (WPK and WRK). The abundance of the CO₄ species in the kimberlite melt is between that of the reduced and oxidized pyrolite models. It has less CO₄ than the fully oxidized pyrolite. In contrast, more reduced systems like pyrolite+CO and carbonatitic-type melts (Solomatova and Caracas, 2019; Davis et al., 2022) have less CO₃ and CO₂, and more CO (Fig. 5). Compared to kimberlite where 96% of bonds are C-O bonds, a larger part of the C-O bonding in more reduced systems is replaced with C–C, C-Fe, and C-Si bonding, which eliminates those species from being expressed in bonding plots that only show C-O species. These results indicate that the redox state is the primary control of carbon speciation, rather than carbon content or the presence of water. The abundance of CO₄ species in kimberlite compositions lies between those of reduced and oxidized pyrolite melts, remaining lower than in fully oxidized systems.

At pressures above 3 GPa, carbonate groups account for more than 70% of all C-O species across all analyzed compositions (Fig. 5). The rise in CO₃ groups below 1 GPa is compensated by a decrease in CO₂ and CO abundances. In the WRK composition, CO₂ species decrease from 54% of all C-O species at 1.4 GPa to less than 10% at 3 GPa. Additionally, the presence of water in the kimberlite melt composition enhances CO₂ formation (Fig. 5), matching previous observations (Dixon and Stolper, 1995). For example, CO₂ abundance at 0.8 GPa is 19% higher in the WRK composition compared to the CRK composition. The difference between the anhydrous and hydrous compositions is also reflected in the decay of CO₂ abundances with pressure. In anhydrous compositions, CO₂ abundances fall off to less than 10% by 1 GPa, but only reach the same abundance at 1.8 and 3 GPa for the WPK and WRK compositions, respectively. Finally, CO₄ groups only begin to exist at higher pressures. At 5 GPa, they are the least common C-O species, but become the second-most prevalent above ~20 GPa.

In Fig. 6b, we compare C-O coordination state abundances across different isotherms for a CRK composition. Species abundances show similar trends regardless of temperature. Notably, at the hottest simulated temperature, 3000 K, there is a spike in CO₂ abundance, indicating preferential formation of the gas phase at higher temperatures. The increase in CO₂ abundance is compensated by a drop in CO₃ abundance.
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Relative abundances of hydrogen-oxygen species as a function of pressure in the water-poor kimberlite melt (WPK) and water-rich kimberlite melt (WRK) at 3000 K. CN stands for coordination number of hydrogen. The plots show the dominance of hydroxyl groups (OH).</figcaption></figure>

The effect of temperature on C-O species is also apparent in their lifetimes (Figure S5). For instance, CO₃ species have average lifetimes of 2829 fs at 2.09 GPa and 2500 K. At 2.3 GPa and 3000 K, the average lifetime falls to 1629 fs. Similarly, the lifetimes of the CO₂ species decrease at higher temperatures, with an average lifetime of 335 fs at 2.09 GPa and 2500 K and 225 fs at 2.3 GPa and 3000 K.

Our lifetime analysis for the water-poor kimberlite at 1.4 GPa and 3000 K reveals the exceptional stability of the carbonate ion (CO₃), which has an average lifetime of 1279 fs (1.28 ps) (Fig. 4). This value is more than double the average lifetime of the SiO₄ tetrahedron (570 fs), the primary network-forming unit in the melt. This establishes the carbonate ion as the most robust and durable polyanionic complex in the kimberlite melt. Other carbon species show varied stability: molecular CO is also a relatively stable component with a lifetime of 505 fs. In contrast, molecular CO₂ (332 fs) is less persistent, while the four-fold coordinated CO₄ species is transient (38 fs). This hierarchy of lifetimes definitively demonstrates that dissolved carbon exists as a suite of distinct chemical species, with the carbonate ion in particular acting as a long-lived, integral component of the melt's structure.

Our lifetime analysis of the anhydrous kimberlite (carbon-rich kimberlite melt (CRK)) at 2500 K and 3000 K (Figure S5) reveals that carbonate ion (CO₃) consistently emerges as the most stable species, with lifetimes often reaching several picoseconds (e.g., 3968 fs at 1.64 GPa, 2500 K). In contrast, molecular CO and CO₂ exhibit much shorter lifetimes, typically in the hundreds of femtoseconds, while CO₄ remains a transient, high-pressure species with limited stability.

Temperature exerts a systematic destabilizing effect: raising the temperature from 2500 K to 3000 K reduces the lifetimes of all carbon species. For instance, at ~2 GPa, the CO₃ lifetime decreases from 2829 fs to 1917 fs, consistent with enhanced bond-breaking dynamics at elevated thermal energy.

Pressure effects are more complex and species-dependent. CO₃ lifetimes show a non-monotonic trend, peaking at low to moderate pressures (~1–4 GPa) before declining at higher pressures, suggesting an optimal stability window in the shallow upper mantle. In contrast, the lifetimes of CO and CO₂ generally decrease with increasing pressure, reflecting their compressible molecular nature. CO₄, although short-lived, becomes more stable with pressure, consistent with its high-pressure character.

Because simulation durations vary across pressure-temperature conditions, these observations reflect qualitative rather than strictly quantitative trends. Nevertheless, the results demonstrate that both the abundance and stability of carbon species in kimberlite melts are strongly modulated by pressure and temperature, with the carbonate

ion acting as the dominant structural component under upper mantle conditions.

The $\frac{N_{CO_2}}{N_{CO_2} + N_{CO_3^{2-}}}$ ratio is an indication of gas formation due to decomposition of CO₃ groups. We compare the $\frac{N_{CO_2}}{N_{CO_2} + N_{CO_3^{2-}}}$ ratio to literature values for carbon-bearing pyrolite (Solomatova and Caracas, 2019), MORB (Guillot and Sator, 2011), enstatite (Ghosh et al., 2017), and forsterite (Solomatova et al., 2020) (Fig. 6c). At 2.5 GPa, we obtain a ratio of 0.07 for both the CPK and CRK compositions. This matches the value for forsterite + 7 wt% CO (Solomatova et al., 2020), and is very close to the ratio for enstatite + 16 wt% CO₂ (Ghosh et al., 2017) at 0.11 GPa. Pyrolite + 5 wt% CO₂ (Solomatova and Caracas, 2019) follows a similar trend with pressure but tends to have higher values, reaching 0.26 at the same pressure conditions. These compositions, which do not integrate carbon into the silicate melt structure as efficiently as the kimberlites (see Section 3.3), force carbon out of the melt and into the vapor phase, resulting in higher values. Carbonated MORB melts (Guillot and Sator, 2011) show the highest values, reaching 0.67 at 2.5 GPa. These results need further validation, as they were obtained with empirical potentials, which may have inconsistencies at high pressures.

### 3.2.6. *Hydrogen-oxygen*

Hydrogen bonding is dominated by H-O bonds. As expected for the pressure range studied here, hydrogen predominantly exhibits one-fold coordination with oxygen. A more interesting parameter to instead examine is the ratio of OH to H₂O in the melt, as the incorporation method of H affects the polymerization and the density of the melt. We expect that the OH to H₂O ratio decreases with increasing pressure at higher-pressure conditions (Mookherjee et al., 2008; Karki et al., 2021).

When water molecules enter a silicate melt, most of them dissociate into OH + H. In this way, they break the Si-O-Si sequence of bonds. Two silicate tetrahedra that share one oxygen can be separated, both ending up with one OH group. The breakup of the water molecule removes one bridging oxygen and creates two non-bridging oxygens (NBOs). This depolymerizes the melt, reducing viscosity and increasing the NBO/T ratio (Mysen, 2014). In contrast, undissociated molecular water (H₂O) interacts primarily through hydrogen bonding with other water molecules or oxygen atoms in the melt structure. They do not directly form covalent bonds with the silicate network, resulting in less dramatic changes in melt properties (Schoize, 1960). Throughout the pressure-temperature range investigated here, more than 80% of hydrogen-oxygen species are OH, with the remainder present as H₂O (Fig. 7), regardless of water content. Other transient hydrogen-oxygen species such as H₃O⁺ and HO₂, appear only in trace amounts and with very short lifetimes. The abundance of molecular water increases with pressure. The predominance of OH groups suggests that introducing water primarily alters the kimberlite melt structure by breaking Si-O-Si bonds and promoting depolymerization.

Hydrogen-oxygen lifetimes provide quantitative support for these distinct structural roles. At 2.3 GPa and 3000 K in the WPK melt, hydroxyl groups have an average lifetime of 231 fs, more than twice that of H₂O species with average lifetime of 85 fs. In terms of frequency, OH groups occur approximately ten times more often than H₂O. This difference in both stability and number of occurences highlights that the dominant and structurally significant role of water in kimberlite melts is through its dissociation into hydroxyl groups, which actively depolymerize the silicate network.

### 3.2.7. *Minor elements*

**Titanium-oxygen.** The CN for titanium bonded to oxygen is 4 at 0 GPa, and changes to 5 at 2 GPa for all compositions (Figure S6). CNs increase more rapidly for titanium than for aluminum. Most titanium transitions to sixfold coordination between 9 and 12 GPa. This shift occurs at higher pressures for the DK composition and lower pressures for volatile-bearing compositions.
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** (a) The ratio of non-bridging oxygens (NBO) per tetrahedrally coordinated cation (where T represents Si and Al) is an indicator of silica network polymerization. Lower values indicate higher polymerization. Our results show that carbon-bearing compositions are more polymerized (lower NBO/T), while hydrogen-bearing compositions are depolymerized (higher NBO/T values). These trends are compared with NBO/T values of kimberlite rocks, where the green shaded region represents a range of Jericho kimberlites (Kopylova et al., 2007; Price et al., 2000), and the beige shaded region represents different localities from Southern Africa (Becker and Roex, 2006; Le Roex et al., 2003; Harris et al., 2004). Note that the Jericho and Southern African kimberlite values are recorded only at ambient conditions; the pressure effect is unknown. The values have been extended to higher pressures for illustrative purposes. Snapshots of (b) the dry kimberlite melt structure at 2 GPa and 3000 K, and (c) the water-rich kimberlite melt at 2.1 GPa and 3000 K. The largest silica clusters seen in the snapshots are highlighted in the melt and shown separately. Note that in the more depolymerized water-bearing composition, the chains of interconnected silica tetrahedra tend to be shorter.</figcaption></figure>

**Aluminum-oxygen.** The CN of aluminum to oxygen increases from 4 to 6 between 0 and 40 GPa (Figure S6). The trend is similar across all compositions. Similarly to silicon, a large proportion of aluminum forms tetrahedra at low pressure. However, $AlO_5$ becomes the most abundant species by 15 GPa for all compositions, while $SiO_4$ tetrahedra are the most abundant species well over 20 GPa.

**Potassium-oxygen.** The CN of potassium to oxygen displays more variability and a higher CN compared to the other atomic pairs examined (Figure S6). This variability is primarily due to statistical limitations: each simulation contains only one potassium atom, which amplifies the effect of local structural fluctuations and limits the robustness of the averaged CN value. Consequently, even minor thermal perturbations or transient structural rearrangements can lead to large variation in the reported CN. The calculated values should therefore be taken as qualitative estimates rather than a quantitative determination. Across all compositions, there are variation in the CN between 0 and 3 GPa, after which the CN varies between 8 and 12.

### 3.3. Silica polymerization

The degree of polymerization in silicate melts directly affects their physical and transport properties. More depolymerized melts have

lower viscosities (Mysen et al., 1982), denser and more compressible structures (Wang et al., 2014), and higher diffusivities and electrical conductivities (Wolf and McMillan, 1995; Wang et al., 2014). Kimberlite melts are characterized by very high Mg:Si ratios; as a result, they should be highly depolymerized.

To quantify the effect of pressure and composition on silicon polymerization, we calculated the number of non-bridging oxygens per tetrahedrally coordinated cation NBO/T (Fig. 8a). NBO/T, as defined by Mysen and Virgo (Mysen and Virgo, 1980; Mysen et al., 1982; Mysen, 1990), indicates degree of polymerization in silicate melts. Although this is a highly empirical parameter, it remains widely used in the literature to describe the degree of polymerization in melts, glasses, and minerals. For each step in our simulation, we count the number of NBOs and divide by the total number of silicon and aluminum atoms. Although titanium is also fourfold coordinated at low pressures, its rapidly changing coordination under pressure increases the scattering of the results. As only two titanium atoms are present in our compositions, we choose to exclude them from the NBO/T calculation. We then calculate the average NBO/T value for each simulation. A completely depolymerized silicate structure has an NBO/T value of 4, and this value decreases for more polymerized phases.
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Size distribution of silicate and silico-carbonate linked clusters in various kimberlite melts as a function of pressure. The y-axis represents the number of linked Si and C cations in a given cluster. The color scale indicates the cluster's lifetime. **\*\*This figure directly visualizes the opposing effects of volatiles: carbon promotes the formation of larger, more stable polymeric chains, while water leads to smaller, shorter-lived clusters.**</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** The percentage of carbon exsolved into a vapor phase as a function of pressure and depth during magma ascent. The onset of bubble nucleation, indicated by the start of the arrows in the inset, is strongly dependent on the melt's volatile content. Compositions containing water (WRK and glisten-shortWPK) begin to degas at higher pressures (deeper in the mantle) than anhydrous compositions (CPK and CRK). Once nucleated, the vapor fraction increases exponentially with decreasing pressure. The data also indicate that higher temperatures (3000 K vs. 2500 K) promote more extensive vaporization.</figcaption></figure>

We assume that all Si and Al atoms act as network formers, regardless of their coordination number. However, at highest plotted pressures (e.g., 11 GPa), 15% of the Si atoms exist in fivefold coordination ($SiO_5$). These species have average lifetimes comparable to those of network-modifying cations such as Mg, Ca, and Fe, suggesting that $SiO_5$ units may not function as stable network formers in the melt. Instead, they may behave more like transient, flux-like species, facilitating structural rearrangements rather than contributing to long-lived polymerized networks.

Therefore, our current method of calculating *T* as the total number of Si + Al atoms, irrespective of coordination, likely leads to a systematic overestimation of the number of network-forming units, and hence an underestimation of NBO/T values at high pressures. While this limitation does not affect the qualitative trends we observe with

pressure and composition, it highlights the need for caution when interpreting absolute NBO/T values at highest pressure.

In the DK composition, the NBO/T values vary from 2.85 to 3.1 between 0 and 11 GPa, indicating a highly depolymerized melt even without volatiles. The addition of carbon seems to have a polymerizing effect on the silicate network, particularly at pressures greater than 2 GPa, where $CO_3$ groups are more prevalent (Fig. 5). This aligns with previous experiments, where carbonate groups were observed to be directly incorporated into the silicate network of depolymerized melts, adding to their polymerization (Brey, 1976). However, the addition of water increases the NBO/T, indicating depolymerization of the melt. This effect is particularly pronounced in the WRK. The water molecules split into H and OH groups that separate the chains of silica tetrahedra (Fig. 7).

We compare our results with the estimated NBO/T ratios of kimberlites from different localities; these values are obtained from the bulk chemical analysis. For five measured kimberlite samples from the Jericho mine (Kopylova et al., 2007; Price et al., 2000), NBO/T values vary between 2.77 and 3.11 at 0 GPa. This range of values completely encompasses the NBO/T values for all our simulated compositions at similar conditions. Estimates of southern African kimberlite melts have a higher degree of polymerization than Canadian kimberlites, with NBO/T values ranging from 2.57 to 2.75 at 0 GPa (Becker and Roex, 2006; Le Roex et al., 2003; Harris et al., 2004). Southern African kimberlites with similar water content as Jericho kimberlites (~9 wt%) contain more $CO_2$ (14 to 29 wt%) and may be more polymerized due to a higher percentage of carbon.

An example of the depolymerizing effect of water is depicted in Fig. 8b-c. In these snapshots of the WRK and DK compositions at 2–2.1 GPa and 3000 K, we display the connectivity of the silicate network. Silica polyhedra, represented as dark blue tetrahedra, form shorter interconnected chains in the WRK composition, highlighting the disrupting effect of water on polymerization. In these snapshots, the largest cluster formed in the DK composition consists of four tetrahedra, while the largest cluster formed in the WRK composition contains only three tetrahedra. In CRK, the NBO/T is higher than in CPK at pressures above 6 GPa, which appears counterintuitive given carbon's typical polymerizing role. However, CRK contains approximately 15% more $CO_3$ groups than CPK. This elevated carbon content may lead to saturation of polymerization sites within the silicate network. Once these
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Atomistic snapshots of carbon-rich kimberlite (CRK, left) and water-rich kimberlite (WRK, right) compositions at 0.6 GPa and 3000 K. Each panel displays a representative simulation cell, with atomic species colored as indicated in the legend below. The snapshots capture phase separation between the silicate melt and the vapor phase. In the CRK composition, CO and CO₂ molecules are clearly visible in the vapor, whereas in the WRK composition CO₂ and H₂O species are present in the vapor phase.</figcaption></figure>

sites are saturated, additional CO₃ groups may no longer contribute to network connectivity. Instead, they could begin to interfere with existing linkages, resulting in a net depolymerization effect. This could explain the observed increase in NBO/T despite the higher carbon content.

Finally, to further quantify the effect of pressure and composition on silicon polymerization, we identify SiₓO_z and SiₓC_yO_z clusters in the different melts, which represent interconnected chains of polyhedral units that constitute the backbone of the melt structure. We define the size of these clusters by the number of cations, either x silicon atoms (for the DK composition) or the sum x+y of silicon and carbon atoms (for the rest of the melts). Comparisons of these numbers between melt compositions yield the effect of water and carbon on melt polymerization. We plot the results in Fig. 9. We find that the size of the polymers increases with increasing pressure. In the carbon-bearing compositions, polymers consisting of silicon and carbon polyhedra tend to be larger than the silicon-only polymers observed in the DK composition. As CO₃ groups are incorporated into the melt network, they extend and glue together clusters of SiO₄ groups, resulting in a more polymerized structure. In contrast, water has a depolymerizing effect. Polymers consisting of silicon and carbon polyhedra, when compared directly to the CRK composition, are smaller in both the WPK and WRK compositions, indicating that the addition of even a small amount of water has a disruptive effect on the polymerized melt network.

### 3.4. Vaporization

Upon decompression, at sufficiently low densities, nanocavities spontaneously separate in our simulation cells. Small atomic species populate these cavities. The cavities represent nanobubbles, and the atomic species represent gas phases (Fig. 11). Gas species are identified through examination of the connectivity matrix (see Methods). For the DK composition, the temperatures employed here are insufficient to promote the evaporation of atoms at the timescale of our simulations. This is largely due to the refractory nature of the elements in the DK composition, but may also be partly due to the strong surface tension associated with the small sizes of the cavities, a limitation of our simulations. Atoms can overcome these energy barriers at the scale of our simulations only if they have enough kinetic energy, i.e. if the temperature is sufficiently high.

For the carbonated and the hydrated melts, however, we verify the presence of gas species (Table T1). In the anhydrous carbonated

melts, we see CO, CO₂ and free oxygen in the vapor phase. In hydrous kimberlites, isolated molecules consist primarily of hydrogen-bearing species (H₂O and H), alongside CO and CO₂. Hydroxyl groups (OH) and H₂ are minor components of the vapor phase in both hydrous compositions at less than 1%. The emergence of small amounts of more complex species, such as COH and CH₂O₂, suggests ongoing chemical interactions between C, O, and H in the vapor phase. At higher pressures, complex species involving Fe, Mg, and Al also appear in trace amounts.

Following decompression, the pressure at which gas species first appear depends on the volatile content. Compositions with greater volatile contents first nucleate vapor-filled nanobubbles at higher pressures. For instance, at 2500 K, in the CRK composition the gas species are first observed at 2.1 GPa. At 3000 K, gas species are first observed at 2.6 GPa for the CPK composition, 5.3 GPa for the CRK composition, 2.7 GPa for the WPK composition, and 6.3 GPa for the WRK composition (Fig. 10). For the WPK composition, there is a gap in the data between 2.7 and 7.3 GPa; bubbles may actually nucleate at pressures greater than 2.7 GPa, in line with the volatile content trend. These values correspond to a wide range of depths, from 60 km to 189 km.

However, it is important to note that this study cannot exactly report bubble nucleation pressures for kimberlite melts, but rather presents examples and trends of expected bubble nucleation pressures for the specific compositions studied. This is because carbon solubility depends on many factors, including the quantity of network-forming atoms in the melt (Si and Al) and the carbon and water contents (King and Holloway, 2002; Guillot and Sator, 2011). Real kimberlite melts encompass a range of silica contents (from 18 to 27 wt% SiO₂ (Giuliani et al., 2023)), and more silica-rich melts exhibit increased degassing behavior (Moussallam et al., 2015). As kimberlites ascend, they evolve from more carbonatitic-type magmas to more silica-rich compositions due to the incorporation and assimilation of new mineral phases. Variations in the chemical makeup of the magmas can lead to differences in how they behave with respect to volatile retention (Brooker et al., 2011), making it difficult to predict their degassing behavior. Thus, our reported bubble nucleation pressures illustrate an example pressure range over which bubble nucleation may begin in kimberlite melts, depending on the composition and the volatile concentration.

In general, our calculations predict bubble nucleation at pressures higher than previous experimental observations, which report degassing in kimberlite magmas at pressures corresponding to depths
of 4 to 8 km (Moussallam et al., 2015; Brooker et al., 2011). The main source of the discrepancy is most likely due to our procedure of counting the free species compared to the experimental determination of the bubbles. Our simulations treat the melt as a polymer and the gas phase as separate atomic clusters, equivalent to molecules. This prevents us from correctly identifying dissolved molecular gases in the melt. On the other hand, experiments reveal the formation of bubbles only above a certain size threshold, which is significantly larger than the size of our nanoscale-sized bubbles. The geometry of the simulations, i.e. slab vs bulk, and inherent finite-size effects would also add to the discrepancies between observations and modeling, but quantifying their effect is less straightforward.

Further causes are related to the exchange-correlation functional, as PBE may not reliably capture phase coexistence boundaries between the melt and vapor phases. Additionally, PBE tends to underestimate van der Waals interactions and may misrepresent the energetic landscape. Although more advanced functionals (e.g., SCAN or hybrid functionals) may offer improved accuracy, they are computationally prohibitive for the system sizes and timescales necessary to study vaporization processes in silicate melts. Thus, it is likely that in reality, bubble nucleation starts at pressures somewhere in between predictions from simulations and experimental observations.

After the initial nucleation of vapor, the percentage of carbon in nanobubbles increases with decreasing pressure. This increase is more gradual in the more volatile-rich compositions, which begin the nucleation process at greater depths within the Earth. The increase in carbon vaporization at higher water concentrations suggests a reduction in carbon solubility, consistent with previous experimental work that indicates that the presence of water encourages the formation of vapor at elevated pressures (Keppler, 2003).

The generation of carbon-bearing vapor at depth enhances magma buoyancy, as the exsolution of carbon-bearing gas increases the buoyancy force on the liquid, promoting ascent through the mantle. For magmas to erupt, the volume of exsolved volatiles (fluid or gas, depending on the pressure-temperature conditions) must constitute 70%-80% of the total volume of the magma-volatile system (Sparks, 1978). Calculating the volume of atoms in the vapor in the CRK there is 6.3% vapor at 0.61 GPa and 1.7% vapor at 0.56 GPa.

## 4. Conclusions

Kimberlite melts are transitional melts between silicic and carbonatitic melts, with their silicate networks influenced by volatile elements. Carbon and hydrogen, the most common volatile elements, have different but important effects on the structure and speciation of kimberlite melts.

We find that the coordination environments of magnesium, calcium, and iron readily transition from lower to higher coordination states with increasing pressure. Magnesium transitions from predominantly fourfold up to sevenfold coordination around 45 GPa, while calcium shifts from sixfold to tenfold coordination. Iron increases from tetrahedral to octahedral coordination above 30 GPa. The dominant coordination state for each of these cations never exceeds 50% of all cation-oxygen clusters, contributing to the structural diversity of the melt. All three elements form species that display short average lifetimes, ranging from 20 to 75 fs. Such brief lifetimes indicate that magnesium, calcium, and iron do not form long-lived, rigid coordination structures but instead contribute to the overall structural flexibility of the melt, rearranging in the melt as pressure increases. The flexible coordination states and short lifetimes of these elements reinforce their roles as network modifiers. The significant magnesium content in our simulated kimberlite melts, with a Mg:Si ratio of ~ 1.6, contributes to the depolymerized nature of kimberlite melts.

In contrast to magnesium, calcium, and iron, silicon exhibits a more rigid coordination environment, with more than 80% of silicon forming stable SiO₄ tetrahedra up to 9 GPa. With increasing pressure,

around 27 GPa, SiO₅ becomes prevalent. Aluminum, another network former, shows similar tetrahedral coordination at low pressures. SiO₄ dominance continues beyond 20 GPa, whereas AlO₄ is less prevalent than AlO₅ before 15 GPa across all compositions. SiO₄ polyhedra have some of the longest lifetimes in our simulations, reaching 7 picoseconds in the WPK composition at 3000 K, indicating their structural stability within the melt. The prevalence of SiO₄ tetrahedra and their long lifetimes contrasts sharply with the dynamic, short-lived coordination states of magnesium, calcium, and iron and indicates the role of silicon as a network former in kimberlite melts.

Carbon in kimberlite melts primarily bonds with oxygen, as 96% of bonds are C-O bonds. As pressure increases, carbon speciation shifts from isolated CO₂ species to carbonate ions. This transition occurs at 1 GPa in anhydrous melts and 3 GPa in hydrous compositions. The decrease in CO₂ abundance also varies with composition: in anhydrous melts, CO₂ falls to less than 10% of carbon species by 1 GPa, whereas in the WPK and WRK compositions, this threshold is reached at 1.8 and 3 GPa, respectively. Carbon also forms bonds with iron, other carbons and silicon but they represent only a small fraction of the total bonds that carbon is involved in.

The macroscopic properties of silicate melts, such as viscosity and diffusivity, are fundamentally controlled by the degree of polymerization of the silicate network. The degree of polymerization, quantified by the NBO/T ratio and the size of linked silicon and carbon polyhedra, is strongly influenced by the composition of the kimberlite melts and the volatile content. In the DK composition, NBO/T values range from 2.85 to 3.1 between 0 and 11 GPa, indicating a highly depolymerized structure even without volatiles. Dissolved volatiles like H₂O and CO₂ play opposing roles in influencing these properties. Water acts as a potent depolymerizing agent. It reacts with bridging oxygens (Si-O-Si) to form hydroxyl groups (Si-OH), effectively breaking down the silicate framework. This structural depolymerization is known to reduce melt viscosity and enhance the diffusivity of chemical species. Conversely, dissolved CO₂, existing primarily as carbonate complexes (CO₃²⁻), has a net polymerizing effect. Above 2 GPa, anhydrous carbonated kimberlites (CPK and CRK) show increases polymerization compared to the volatile-free melt (DK), as carbonate ions become part of the silicate network.

With decreasing pressure, during magma ascent, nanobubbles form at low densities, containing small atomic species that we assign to the gas phase. In dry kimberlite melts, the refractory nature of elements and high surface tension prevent significant vaporization. However, in carbonated and hydrated melts, gas species emerge, including CO, CO₂, H₂O, H₂, and OH. Water-bearing compositions exhibit increased vaporization. The pressure at which gas phases first appear depends on volatile content, with more volatile-rich compositions nucleating vapor-filled nanobubbles at higher pressures (1.7–3.9 GPa, corresponding to depths of 21–120 km). Carbon vaporization increases as the pressure decreases, with water promoting vapor formation by reducing the solubility of carbon. The exsolution of carbon-bearing vapor enhances magma buoyancy, facilitating ascent.

At the same time, large-scale exsolution of volatiles removes water's strong depolymerizing effect, leading to polymerization of the remaining melt. While this process increases viscosity, the magma's inherently ultrabasic, low-silica composition allows it to remain sufficiently fluid during ascent.

Although we have examined the detailed structural changes in carbon- and water-bearing kimberlite melts, investigating a broader range of kimberlite compositions is essential to fully understand how carbon and hydrogen influence melt structure. The principles established in this study can be extended to predict how the broad compositional spectrum of natural kimberlites will modulate their physical properties and degassing behavior. A higher silica and alumina content, characteristic of more evolved kimberlites, would create a more polymerized melt structure. This not only increases baseline viscosity but, more critically, decreases CO₂ solubility, shifting bubble nucleation to
shallower, lower-pressure regimes more consistent with experimental observations. Conversely, while total alkali content of kimberlites is lower than 5%, melts richer in alkalis (Na₂O, K₂O) would exhibit enhanced CO₂ solubility, retaining volatiles to greater depths. The magma's initial H₂O/CO₂ ratio also plays a kay role, as higher water content simultaneously lowers viscosity while promoting earlier CO₂ exsolution. Therefore, the ascent of a kimberlite is a dynamic process where progressive crustal assimilation continuously alters its chemistry, structure, and degassing threshold, implying that kimberlite eruption is driven by a protracted exsolution process rather than a single nucleation event at a specific depth.

### CRediT authorship contribution statement

**Ana Anzulović:** Writing – original draft, Methodology, Investigation, Formal analysis, Conceptualization. **Anne H. Davis:** Writing – review & editing, Methodology, Formal analysis. **Razvan Caracas:** Writing – review & editing, Software, Resources, Methodology, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

AA acknowledges support from the CompSci project, co-funded by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement N° 945371. We acknowledge support from the Research Council of Norway (NFR) via the Centers of Excellence funding scheme, project number 332523 (PHAB), and Project number 325567-HIDDEN to RC. The simulations were performed on resources provided by Sigma2 - the National Infrastructure for High-Performance Computing and Data Storage in Norway under grant NN9697K.

## Appendix A. Supplementary data

Appendix A contains the convergence of structural parameters for kimberlite compositions at 3000 K (Figures S1, S2, S3), mean-squared displacement data illustrating melt dynamics (Figure S4), lifetimes of carbon-oxygen species (Figure S5), relative abundance of coordination states for minor cations (Figure S6), and vapor species formed from water-rich kimberlite melt (Table T1).

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.gca.2025.10.016.

### Data availability

Data are available through Zenodo at https://doi.org/10.5281/zenodo.16400603.

### References

Bailey, D., 1980. Volatile flux, geotherms, and the generation of the kimberlite-carbonate-alkaline magma spectrum. Miner. Mag. 43 (330), 695-699.

Bajgain, S., Ghosh, D.B., Karki, B.B., 2015. Structure and density of basaltic melts at mantle conditions from first-principles simulations. Nat. Commun. 6 (1), 8578.

Becker, M., Roex, A.P.L., 2006. Geochemistry of South African on-and off-craton, group I and group II kimberlites: petrogenesis and source region evolution. J. Petrol. 47 (4), 673-703.

Blank, J.G., Brooker, R.A., 1994. Experimental studies of carbon dioxide in silicate melts: solubility, speciation, and stable carbon isotope behavior. Rev. Miner. 30, 157-157.

Blöchl, P.E., 1994. Projector augmented-wave method. Phys. Rev. B 50 (24), 17953.

Boulard, E., Gloter, A., Corgne, A., Antonangeli, D., Auzende, A.L., Perrillat, J.P., Guyot, F., Fiquet, G., 2011. New host for carbon in the deep earth. Proc. Natl. Acad. Sci. 108 (13), 5184-5187.

Brey, G., 1976. Carbon dioxide solubility and solubility mechanisms in silicate melts at high pressures. Contrib. Miner. Pet. 57, 215-221.

Brey, G.P., Bulatov, V.K., Girnis, A.V., Lahaye, Y., 2008. Experimental melting of carbonated peridotite at 6-10 GPa. J. Petrol. 49 (4), 797-821.

Brooker, R.A., Sparks, R.S.J., Kavanagh, J.L., Field, M., 2011. The volatile content of hypabyssal kimberlite magmas: some constraints from experiments on natural rock compositions. Bull. Volcanol. 73, 959-981.

Caracas, R., Kobsch, A., Solomatova, N.V., Li, Z., Soubiran, F., Hernandez, J.A., 2021. Analyzing melts and fluids from ab initio molecular dynamics simulations with the UMD package.

Caracas, R., Stewart, S.T., 2023. No magma ocean surface after giant impacts between rocky planets. Earth Planet. Sci. Lett. 608, 118014.

Davis, A.H., Caracas, R., 2024. Degassing of CO₂ triggers large-scale loss of helium from magma oceans. Commun. Earth Environ. 5, 344.

Davis, A.H., Solomatova, N.V., Campbell, A.J., Caracas, R., 2022. The speciation and coordination of a deep earth carbonate-silicate-metal melt. J. Geophys. Res.: Solid Earth 127 (3).

Davis, A.H., Solomatova, N.V., Caracas, R., Campbell, A.J., 2023. Carbon storage in earth's deep interior implied by carbonate-silicate-iron melt miscibility. Geochem. Geophys. Geosystems 24 (8).

de Koker, N.P., Stixrude, L., Karki, B.B., 2008. Thermodynamics, structure, dynamics, and freezing of Mg2SiO4 liquid at high pressure. Geochim. Cosmochim. Acta 72 (5), 1427-1441.

Dixon, J.E., Stolper, E.M., 1995. An experimental study of water and CO₂ solubilities in mid-ocean ridge basaltic liquids. Part II: applications to degassing. J. Petrol. 36 (6), 1633-1646.

Field, M., Scott Smith, B., 1999. Contrasting geology and near-surface emplacement of kimberlite pipes in southern Africa and Canada. In: Proceedings of the 7th International Kimberlite Conference, vol. 1, pp. 214-237.

Ghosh, D.B., Bajgain, S.K., Mookherjee, M., Karki, B.B., 2017. Carbon-bearing silicate melt at deep mantle conditions. Sci. Rep. 7 (1), 848.

Giuliani, A., Schmidt, M.W., Torsvik, T.H., Fedortchouk, Y., 2023. Genesis and evolution of kimberlites. Nat. Rev. Earth Environ. 4 (11), 738-753.

Guillot, B., Sator, N., 2011. Carbon dioxide in silicate melts: A molecular dynamics simulation study. Geochim. Cosmochim. Acta 75 (7), 1829-1857.

Hao, M., Zhou, W.Y., Hrubiak, R., Kenney-Benson, C., Kavanagh, J.L., Davis, W., Zhang, J.S., 2024. The ultralow viscosity of volatile-rich kimberlite magma: Implications for the water content of primitive kimberlite melts. Sci. Adv. 10 (36), 8550.

Harris, M., le Roex, A., Class, C., 2004. Geochemistry of the uintjiesberg kimberlite, South Africa: petrogenesis of an off-craton, group I, kimberlite. Lithos 74 (3-4), 149-165.

Henkelman, G., Arnaldsson, A., Jónsson, H., 2006. A fast and robust algorithm for Bader decomposition of charge density. Comput. Mater. Sci. 36 (3), 354-360.

Holloway, J.R., 1976. Fluids in the evolution of granitic magmas: consequences of finite carbon dioxide solubility. Geol. Soc. Am. Bull. 87 (10), 1513-1518.

Hoover, W.G., 1985. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 31 (3), 1695.

Karki, B.B., Ghosh, D.B., Karato, S.I., 2021. Behavior and properties of water in silicate melts under deep mantle conditions. Sci. Rep. 11 (1), 10588.

Keppler, H., 2003. Water solubility in carbonatite melts. Am. Mineral. 88 (11-12), 1822-1824.

King, P., Holloway, J., 2002. CO₂ solubility and speciation in intermediate (andesitic) melts: the role of water and composition. Geochim. Cosmochim. Acta 66 (9), 1627-1640.

Kohn, W., Sham, L.J., 1965. Self-consistent equations including exchange and correlation effects. Phys. Rev. 140 (4A), A1133.

Konschak, A., Keppler, H., 2014. The speciation of carbon dioxide in silicate melts. Contrib. Miner. Pet. 167, 1-13.

Kopylova, M., Matveev, S., Raudsepp, M., 2007. Searching for parental kimberlite melt. Geochim. Cosmochim. Acta 71 (14), 3616-3629.

Kresse, G., Furthmüller, J., 1996. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 6 (1), 15-50.

Le Roex, A.P., Bell, D.R., Davis, P., 2003. Petrogenesis of group I kimberlites from kimberley, South Africa: evidence from bulk-rock geochemistry. J. Petrol. 44 (12), 2261-2286.

Mitchell, R.H., 2013. Kimberlites: Mineralogy, Geochemistry, and Petrology. Springer Science & Business Media.

Mookherjee, M., Stixrude, L., Karki, B., 2008. Hydrous silicate melt at high pressure. Nature 452 (7190), 983-986.

Morizet, Y., Paris, M., Sifre, D., Di Carlo, I., Ory, S., Gaillard, F., 2017. Towards the reconciliation of viscosity change and carbon dioxide-induced polymerization in silicate melts. Chem. Geol. 458, 38-47.

Morizet, Y., Vuilleumier, R., Paris, M., 2015. A NMR and molecular dynamics study of carbon dioxide-bearing basaltic melts and glasses. Chem. Geol. 418, 89-103.
Moussallam, Y., Florian, P., Corradini, D., Morizet, Y., Sator, N., Vuilleumier, R., Guillot, B., Iacono-Marziano, G., Schmidt, B.C., Gaillard, F., 2016a. The molecular structure of melts along the carbonatite-kimberlite-basalt compositional joint: carbon dioxide and polymerisation. Earth Planet. Sci. Lett. 434, 129-140.

Moussallam, Y., Morizet, Y., Gaillard, F., 2016b. Water–CO₂ solubility in low SiO₂-melts and the unique mode of kimberlite degassing and emplacement. Earth Planet. Sci. Lett. 447, 151–160.

Moussallam, Y., Morizet, Y., Massuyeau, M., Laumonier, M., Gaillard, F., 2015. CO₂ dioxide solubility in kimberlite melts. Chem. Geol. 418, 198–205.

Mysen, B.O., 1990. Effect of pressure, temperature, and bulk composition on the structure and species distribution in depolymerized alkali aluminosilicate melts and quenched melts. J. Geophys. Res.: Solid Earth 95 (B10), 15733–15744.

Mysen, B., 2014. Water-melt interaction in hydrous magmatic systems at high temperature and pressure. Prog. Earth Planet. Sci. 1, 1–18.

Mysen, B.O., Virgo, D., 1980. The solubility behavior of carbon dioxide in melts on the join NaAlSi3O8- CaAl2Si2O8- carbon dioxide at high pressures and temperatures: a Raman spectroscopic study. Am. Mineral. 65 (11–12), 1166–1175.

Mysen, B.O., Virgo, D., Seifert, F.A., 1982. The structure of silicate melts: implications for chemical and physical properties of natural magma. Rev. Geophys. 20 (3), 353–383.

Ni, H., Keppler, H., 2013. Carbon in silicate melts. Rev. Miner. Geochem. 75 (1), 251–287.

Nosé, S., 1984. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 81 (1), 511–519.

Perdew, J.P., Burke, K., Ernzerhof, M., 1996. Generalized gradient approximation made simple. Phys. Rev. Lett. 77 (18), 3865.

Price, S., Russell, J., Kopylova, M., 2000. Primitive magma from the jericho pipe, NWT, Canada: constraints on primary kimberlite melt chemistry. J. Petrol. 41 (6), 789–808.

Priestley, K., McKenzie, D., Debayle, E., 2006. The state of the upper mantle beneath southern Africa. Tectonophysics 416 (1–4), 101–112.

Russell, J.K., Porritt, L.A., Lavallée, Y., Dingwell, D.B., 2012. Kimberlite ascent by assimilation-fuelled buoyancy. Nature 481 (7381), 352–356.

Sanloup, C., Drewitt, J.W., Konôpková, Z., Dalladay-Simpson, P., Morton, D.M., Rai, N., van Westrenen, W., Morgenroth, W., 2013. Structural change in molten basalt at deep mantle conditions. Nature 503 (7474), 104–107.

Sanloup, C., Hudspeth, J.M., Afonina, V., Cochain, B., Konôpková, Z., Lelong, G., Cormier, L., Cavallari, C., 2019. Polymerized 4-fold coordinated carbonate melts in the deep mantle. Front. Earth Sci. 7, 72.

Sanville, E., Kenny, S., Smith, R., Henkelman, G., 2007. Improved grid-based algorithm for bader charge allocation. J. Comput. Chem. 28 (5), 899–908.

Schoize, H., 1960. On the question of the difference between h₂O molecules and OH-groups in glasses and minerals. Naturwissenschaften 47, 226–227.

Solomatova, N., Caracas, R., 2019. Pressure-induced coordination changes in a pyrolicic silicate melt from ab initio molecular dynamics simulations. J. Geophys. Res.: Solid Earth 124 (11), 11232–11250.

Solomatova, N.V., Caracas, R., 2021. Buoyancy and structure of volatile-rich silicate melts. J. Geophys. Res.: Solid Earth 126 (2).

Solomatova, N., Caracas, R., Cohen, R., 2020. Carbon speciation and solubility in silicate melts. Carbon Earth's Inter. 179–194.

Sparks, R., 1978. Gas release rates from pyroclastic flows: a assessment of the role of fluidisation in their emplacement. Bull. Volcanol. 41, 1–9.

Sparks, R., 2013. Kimberlite volcanism. Annu. Rev. Earth Planet. Sci. 41, 497–528.

Sparks, R., Baker, L., Brown, R., Field, M., Schumacher, J., Stripp, G., Walters, A., 2006. Dynamical constraints on kimberlite volcanism. J. Volcanol. Geotherm. Res. 155 (1–2), 18–48.

Sun, N., Stixrude, L., de Koker, N., Karki, B.B., 2011. First principles molecular dynamics simulations of diopside (CaMgSi₂O₆) liquid to high pressure. Geochim. Cosmochim. Acta 75 (13), 3792–3802.

Tang, W., Sanville, E., Henkelman, G., 2009. A grid-based bader analysis algorithm without lattice bias. J. Phys.: Condens. Matter. 21 (8), 084204.

Van der Waals, J., 1873. On the Continuity of the Gaseous and Liquid States. Sijthoff.

Vuilleumier, R., Seitsonen, A.P., Sator, N., Guillot, B., 2015. Carbon dioxide in silicate melts at upper mantle conditions: Insights from atomistic simulations. Chem. Geol. 418, 77–88.

Wang, Y., Sakamaki, T., Skinner, L.B., Jing, Z., Yu, T., Kono, Y., Park, C., Shen, G., Rivers, M.L., Sutton, S.R., 2014. Atomistic insight into viscosity and density of silicate melts under pressure. Nat. Commun. 5 (1), 3241.

Weast, R., 1972. Handbook of Chemistry and Physics, fifty-three ed. Chemical Rubber Pub.

Wolf, G.H., McMillan, P.F., 1995. Pressure effects on silicate melt structure and properties. Rev. Miner. 32, 505–562.

Wyllie, P.J., 1980. The origin of kimberlite. J. Geophys. Res.: Solid Earth 85 (B12), 6902–6910.

Yao, X., Xie, C., Dong, X., Oganov, A.R., Zeng, Q., 2018. Novel high-pressure calcium carbonates. Phys. Rev. B 98 (1), 014108.