Noble gases in carbonate melts: constraints on the solubility and the surface tension by molecular dynamics simulation

TL;DR

This study uses molecular dynamics simulations to quantify noble gas solubility in carbonate melts under various pressures, revealing how composition and pressure influence noble gas behavior relevant to Earth's mantle geochemistry.

Contribution

It provides the first computational estimates of noble gas solubility in carbonate melts across a range of pressures, enhancing understanding of noble gas partitioning in Earth's deep interior.

Findings

Noble gas solubility decreases with larger atomic size and alkaline-earth content.

Solubility in carbonate melts remains comparable to silicates at mantle pressures.

Surface tension between melts and noble gases increases with pressure.

Abstract

Although they are rare elements in the Earth's mantle, noble gases (NG) owe to their strongly varying masses contrasting physical behaviors making them important geochemical tracers. When partial melting occurs at depth, the partitioning of NGs between phases is controlled by a distribution coefficient that can be determined from the solubility of the NGs in each phase. Here we report quantitative calculations of the solubility of He, Ne, Ar and Xe in carbonate melts based on molecular dynamics simulations. The NG solubilities are first calculated in K2CO3-CaCO3 mixtures at 1 bar and favorably compared to the only experimental data available to date. Then we investigate the effect of pressure (up to 6 GPa), focusing on two melt compositions: a dolomitic one and a natrocarbonatitic one (modeling the lava of Ol Doinyo Lengai). The solubility decreases with the amount of alkaline-earth cation in the melt and with the size of the noble gas. In the natrocarbonatitic melt, Henry's law is fulfilled at low pressures (up to ~ 0.1 GPa). At higher pressures the solubility levels off or even starts to diminish smoothly. In contrast, in molten dolomite the effect of pressure is negligible on the studied P range (3-6 GPa). At the pressures of the Earth's mantle, the solubilities of noble gases in carbonate melts are still of the same order of magnitude as the ones in molten silicates (~ 1-10 mol%). This suggests that carbonatitic melts at depth are not preferential carriers of noble gases, even if the dependence with the melt composition is not negligible and has to be evaluated on a case-by-case basis. Finally we evaluate the surface tension at the interface between carbonate melts and noble gases and its evolution with pressure. Whatever the composition of the melt and of the NG phase, the surface tension increases when P increases from 0 to 6 GPa.

Tables

Table 1. Table A.1: Buckingham parameters for the NG-NG interactions. Guillot and Sator ( 2012 )

NG	$A_{ij}$ (kJ/mol)	${\rho }_{ij}$ (Å)	$C_{ij}$ (Å6/mol)
He	132917.0	0.2051	109.84
Ne	684325.4	0.2083	523.19
Ar	2947863.0	0.2485	5607.64
Xe	4474435.5	0.2940	27142.12

Table 2. Table A.2: Intermolecular potential parameters and partial charges for the carbonate melts. Intramolecular repulsion parameters between oxygen atoms are A OO i ​ n ​ t ​ r ​ a − CO 3 = 2611707.2 superscript subscript 𝐴 OO 𝑖 𝑛 𝑡 𝑟 𝑎 subscript CO 3 2611707.2 A_{{\rm{OO}}}^{intra-{\rm{CO}_{3}}}=2611707.2  kJ/mol and ρ OO i ​ n ​ t ​ r ​ a − CO 3 = 0.22 superscript subscript 𝜌 OO 𝑖 𝑛 𝑡 𝑟 𝑎 subscript CO 3 0.22 \rho_{{\rm{OO}}}^{intra-{\rm{CO}_{3}}}=0.22  Å. Desmaele  et al. ( 2019a , b )

$i$	$j$	$A_{ij}$ (kJ/mol)	${\rho }_{ij}$ (Å)	$C_{ij}$ (Å6/mol)	$q_i$ (e)
Mg	O	243 000	0.24335	1 439	+1.64202
Ca	O	200 000	0.2935	5 000	+1.64202
Na	O	1 100 000	0.2228	3 000	+0.82101
K	O	900 000	0.2570	7 000	+0.82101
C	O	0	1	0	+1.04085
O	O	500 000	0.252525	2 300	$-$0.89429

Table 3. Table A.3: Lennard-Jones parameters for noble gas - carbonate melt interactions with ϵ italic-ϵ \epsilon in kJ/mol and σ 𝜎 \sigma in Å. Guillot and Sator ( 2012 ); Aubry, Sator, and Guillot ( 2013 )

	${ϵ}_{\mathrm{He}}$	${\sigma }_{\mathrm{He}}$	${ϵ}_{\mathrm{Ne}}$	${\sigma }_{\mathrm{Ne}}$	${ϵ}_{\mathrm{Ar}}$	${\sigma }_{\mathrm{Ar}}$	${ϵ}_{\mathrm{Xe}}$	${\sigma }_{\mathrm{Xe}}$
Mg	1.554	2.076	2.896	2.125	3.897	2.455	4.068	2.746
Ca	1.142	2.465	2.142	2.513	3.389	2.843	3.943	3.134
Na	0.523	2.552	0.984	2.600	1.596	2.930	1.891	3.221
K	0.351	3.008	0.663	3.056	1.232	3.386	1.603	3.677
O	0.525	2.81	0.671	2.89	1.131	3.19	1.317	3.44

Table 4. Table B.1: Solubility parameters of NGs in melts: \ce K2CO3– \ce CaCO3 mixtures with different molar ratios of \ce K2CO3 at 1173 K and 1 bar and natrocarbonatite (Na 1.1 K 0.18 Ca 0.36 CO 3 , NKC) at 1600 K, calculated by MD using the TPM.

	${\rho }_m$	${\Gamma }_m^{\mathrm{He}}$	${\Gamma }_m^{\mathrm{Ne}}$	${\Gamma }_m^{\mathrm{Ar}}$	${\Gamma }_m^{\mathrm{Xe}}$
	(mol/L)	($\times {10}^3$)	($\times {10}^3$)	($\times {10}^3$)	($\times {10}^3$)
KC
$x_{K_2{\mathrm{CO}}_3}=0$	24.08	$2.69\pm 0.03$	$1.74\pm 0.04$	$0.27\pm 0.01$	$0.059\pm 0.05$
$x_{K_2{\mathrm{CO}}_3}=0.25$	20.12	$6.70\pm 0.02$	$4.43\pm 0.02$	$0.96\pm 0.02$	$0.26\pm 0.02$
$x_{K_2{\mathrm{CO}}_3}=0.33$	19.39	$6.65\pm 0.03$	$4.24\pm 0.03$	$0.80\pm 0.02$	$0.17\pm 0.03$
$x_{K_2{\mathrm{CO}}_3}=0.4$	19.67	$7.29\pm 0.05$	$4.60\pm 0.04$	$0.83\pm 0.02$	$0.17\pm 0.02$
$x_{K_2{\mathrm{CO}}_3}=0.5$	17.66	$8.63\pm 0.03$	$5.47\pm 0.03$	$1.08\pm 0.03$	$0.26\pm 0.03$
$x_{K_2{\mathrm{CO}}_3}=1$	13.77	$16.00\pm 0.07$	$10.62\pm 0.06$	$2.66\pm 0.07$	$0.79\pm 0.03$
NKC
$P=1$ bar	17.64	$27.97\pm 0.05$	$20.90\pm 0.05$	$6.95\pm 0.06$	$3.22\pm 0.06$
$P=0.1$ GPa	17.98	$22.91\pm 0.04$	$16.66\pm 0.04$	$4.94\pm 0.03$	$2.09\pm 0.03$

Table 5. Table C.1: Solubility of NGs in melts: natrocarbonatite (Na 1.1 K 0.18 Ca 0.36 CO 3 , NKC) at 1600 K and dolomite (Ca 0.5 Mg 0.5 CO 3 , CM) at 1623 K, calculated using the EIM.

	$P$	$X\text{<sub>He</sub>}$	$X\text{<sub>Ne</sub>}$	$X\text{<sub>Ar</sub>}$	$X\text{<sub>Xe</sub>}$
	(GPa)	(mol%)	(mol%)	(mol%)	(mol%)
NKC	0.1	1.5 $\pm$ 0.05	0.84 $\pm$ 0.1	0.28 $\pm$ 0.1	0.16 $\pm$ 0.1
	1.0	7.5 $\pm$ 0.5	5.2 $\pm$ 0.20	1.7 $\pm$ 0.1	0.5 $\pm$ 0.15
	2.0	9.0 $\pm$ 0.2	6.6 $\pm$ 0.4	1.8 $\pm$ 0.10	1.0 $\pm$ 0.15
	4.5	8.3 $\pm$ 0.2	7.3 $\pm$ 0.3	1.95 $\pm$ 0.15	1.4 $\pm$ 0.15
	6.0	7.5 $\pm$ 0.5	7.4 $\pm$ 0.3	1.35 $\pm$ 0.1	1.6 $\pm$ 0.15
CM	3.0	3.1 $\pm$ 0.4	2.7 $\pm$ 0.4	0.45 $\pm$ 0.05	0.12 $\pm$ 0.03
	4.5	3.1 $\pm$ 0.4	3.0 $\pm$ 0.3	0.55 $\pm$ 0.05	0.17 $\pm$ 0.06
	6.0	2.7 $\pm$ 0.2	2.9 $\pm$ 0.4	0.55 $\pm$ 0.1	0.15 $\pm$ 0.05

Table 6. Table C.2: Surface tension between NG fluids and melts: natrocarbonatite (Na 1.1 K 0.18 Ca 0.36 CO 3 , NKC) at 1600 K and dolomite (Ca 0.5 Mg 0.5 CO 3 , CM) at 1623 K.

	$P$	$\gamma \text{<sub>He</sub>}$	$\gamma \text{<sub>Ne</sub>}$	$\gamma \text{<sub>Ar</sub>}$	$\gamma \text{<sub>Xe</sub>}$
	(GPa)	(mN/m)	(mN/m)	(mN/m)	(mN/m)
NKC	0.1	206 $\pm$ 2	182 $\pm$ 3	171 $\pm$ 3	195 $\pm$ 5
	1.0	229 $\pm$ 7	202 $\pm$ 7	225 $\pm$ 6	222 $\pm$ 6
	2.0	242 $\pm$ 4	250 $\pm$ 8	259 $\pm$ 4	252 $\pm$ 5
	4.5	304 $\pm$ 4	277 $\pm$ 4	314 $\pm$ 10	266 $\pm$ 5
	6.0	363 $\pm$ 10	309 $\pm$ 15	326 $\pm$ 6	293 $\pm$ 8
CM	3.0	398 $\pm$ 7	381 $\pm$ 9	357 $\pm$ 6	364 $\pm$ 8
	4.5	444 $\pm$ 9	388 $\pm$ 7	406 $\pm$ 12	382 $\pm$ 11
	6.0	468 $\pm$ 10	435 $\pm$ 9	413 $\pm$ 11	407 $\pm$ 10

Equations

$$X=\frac{\frac{{\rho }_v}{{\rho }_m}\cdot \frac{{\Gamma }_m}{{\Gamma }_v}}{1+\frac{{\rho }_v}{{\rho }_m}\cdot \frac{{\Gamma }_m}{{\Gamma }_v}},$$

$$X=P/k_H,$$

$${\mu }^{ex,i}=-k_{B}T\ln ⟨{\mathrm{e}}^{-{\Psi }_i/k_{B}T}{⟩}_0,$$

$$\begin{array}{rl}\Gamma & ={\mathrm{e}}^{-{\mu }^{ex}/k_{B}T}=⟨{\text{e}}^{-\Psi /k_{B}T}⟩\\ & =\frac{1}{N_{\text{test}}}\sum _i^{N_{\text{test}}}{\text{e}}^{-{\Psi }_i/k_{B}T},\end{array}$$

$$\begin{split}&\gamma=\underset{t\to\infty}{\lim}\gamma(t)\,,\,\\ &\text{with}\,\gamma(t)=\frac{L_{z}}{2}\Big{\langle}\Pi_{zz}-\frac{1}{2}\big{(}\Pi_{xx}+\Pi_{yy}\big{)}\Big{\rangle}\enspace,\end{split}$$

$$V_{ij}(r_{ij})=A_{ij}\exp (-r_{ij}/{\rho }_{ij})-C_{ij}/r_{ij}^6.$$

$$V_{ij}(r_{ij})=4{\epsilon_{ij}}\Big{(}\big{(}{\sigma_{ij}}/{{r}_{ij}}\big{)}^{12}-\big{(}{\sigma_{ij}}/{{r}_{ij}}\big{)}^{6}\Big{)}\enspace,$$

$${\mu }_v={\mu }_m.$$

$${\mu }_{v,m}^{IG}=C+k_{B}T\ln {\rho }_{v,m}^g,$$

$$\frac{{\rho }_m^g}{{\rho }_v^g}={\mathrm{e}}^{-\left({\mu }_m^{ex}-{\mu }_v^{ex}\right)/k_{B}T}=\frac{{\Gamma }_m}{{\Gamma }_v},$$

$$X=\frac{\frac{{\rho }_v}{{\rho }_m}\cdot \frac{{\Gamma }_m}{{\Gamma }_v}}{1+\frac{{\rho }_v}{{\rho }_m}\cdot \frac{{\Gamma }_m}{{\Gamma }_v}},$$

$$X=P{\Gamma }_m/{\rho }_m{k}_{B}T,$$

Full text

Here we report quantitative calculations of the solubility of He, Ne, Ar and Xe in carbonate melts based on molecular dynamics simulations. The NG solubilities are first calculated in \ceK2CO3-\ceCaCO3 mixtures at 1 bar and favorably compared to the only experimental data available to date. Then we investigate the effect of pressure (up to 6 GPa), focusing on two melt compositions: a dolomitic one and a natrocarbonatitic one (modeling the lava emitted at Ol Doinyo Lengai). The solubility decreases with the amount of alkaline-earth cation in the melt and with the size of the noble gas. In the natrocarbonatitic melt, Henry’s law is fulfilled at low pressures (up to $\sim 0.1$ GPa). At higher pressures (a few GPa) the solubility levels off or even starts to diminish smoothly (for He at $P>2$ GPa and Ar at $P>4$ GPa). In contrast, in molten dolomite the effect of pressure is negligible on the studied $P$ range ($3-6$ GPa). At the pressures of the Earth’s mantle, the solubilities of noble gases in carbonate melts are still of the same order of magnitude as the ones in molten silicates ($10^{0}-10^{1}$ mol%). This suggests that carbonatitic melts at depth are not preferential carriers of noble gases, even if the dependence with the melt composition is not negligible and has to be evaluated on a case-by-case basis.

Finally we evaluate the surface tension at the interface between carbonate melts and noble gases and its evolution with pressure. Whatever the composition of the melt and of the NG phase, the surface tension increases (by a factor $\sim 2$) when $P$ increases from 0 to 6 GPa. This behavior contrasts with the situation occurring when \ceH2O is in contact with silicate melts (then surface tension drops when pressure increases to a few GPa).

molten carbonates, noble gases, solubility, surface tension

I Introduction

By the time of accretion, noble gases were already present in proto-Earth material.Marty (2012) Their concentration in the mantle has since evolved through the competing effect of volcanic degassing and radioactive decay. Although they are inert species their strongly varying masses in the series from He to Rn confer them contrasting physical behaviors in the Earth’s mantle. Hence these elements (and their isotopes even more so) constitute important geochemical tracers.Moreira (2013) More specifically, noble gases are incompatible in mantle rocks and tend to partition into the liquid phase when partial melting occurs at depth. In the presence of two immiscible liquids, the partioning is determined by a distribution coefficient, that is related to the ratio of the solubility of the noble gas in each phase. For example a carbonatite-silicate immiscibility is currently the most probable scenario for the genesis of the natrocarbonatites emitted by Ol Doinyo Lengai in Tanzania Fischer et al. (2009). In a more general perspective, the noble gas systematics provides an insight into the past and present dynamics of the carbon-bearing phases (silicates and/or carbonates) in the Earth’s mantle (e.g. the carbon cycle). Hence solubility data under high pressure and in different magmatic liquids, are requisite in order to study these phenomena. Whereas many studies have been devoted to the solubility of noble gases in molten silicates (mainly at low pressures),White, Brearley, and Montana (1989); Montana et al. (1993); Carroll and Stolper (1993); Chamorro-Perez, Gillet, and Jambon (1996); Schmidt and Keppler (2002); Bouhifd, Jephcoat, and Kelley (2008) very little is known concerning their behavior in carbonate melts which are yet of particular interest.Dasgupta and Hirschmann (2006, 2010); Burnard, Toplis, and Medynski (2010); Hammouda and Keshav (2015)

To date, the only experimental constraints are given by Burnard *et al.*Burnard, Toplis, and Medynski (2010), who measured the solubility of noble gases in quenched carbonate melts. In this approach the concentration in the liquid is assumed to be the same as in the glass resulting from the quenching. Moreover the quenching is considered to be fast enough to avoid crystal nucleation. The study reported the solubility of He and Ar in some \ceK2CO3-\ceCaCO3 mixtures at 1 bar between 1123 and 1223 K. The authors also attempted to study Mg-bearing carbonate melts, but could not succeed in quenching them into a glass.

By contrast, in a molecular dynamics (MD) simulation the liquid phase can always be studied, regardless of its composition (and of the pressure/temperature condition). The relevance of the liquid properties calculated from a MD simulation relies on the quality of the implemented interaction potential or force field (FF). The adjustment of the FF generally is system-specific. In a previous study we have presented a FF to model carbonate melts in the \ceCaCO3–\ceMgCO3–\ceK2CO3–\ceNa2CO3–\ceLi2CO3 system, and demonstrated that it leads to an accurate reproduction of their thermodynamics, structure and transport properties Desmaele et al. (2019a, b). Here we perform molecular dynamics simulations based on this FF to study the solubility of noble gases (from He to Xe) into \ceK2CO3-\ceCaCO3 mixtures (for comparison with the data of Burnard *et al.*Burnard, Toplis, and Medynski (2010)) and into dolomitic and natrocarbonatitic melts to investigate the evolution with pressure (up to 6 GPa). First we present in section II the method to evaluate the noble gas solubility. The results are discussed in section III. The surface tension of carbonate melts in contact with a noble gas fluid is evaluated and discussed in section IV.

II Method

II.1 Force Field

To describe the interactions within the carbonate melt we used the interatomic FF proposed by Desmaele *et al.*Desmaele et al. (2019a, b) The interactions between noble gas (NG) atoms were modeled by the FF used by Guillot and SatorGuillot and Sator (2012). These Buckingham potentials (see Eq. (A.1) in Appendix A) were adjusted on the accurate potential energy surface derived by Tang and ToenniesTang and Toennies (2003). As for the interactions between the NG and the carbonate melt elements, they are described by two contributions, a NG-cation and a NG-carbonate one, both of them being modeled by Lennard-Jones (LJ) potentials (Eq. (A.2)). For the NG-cation interactions we used the LJ parameters determined by Guillot and Sator (2012) for the {silicate + NG} system. For the interaction between the NG and the oxygen atoms of the carbonate ions, we retained the parameters established by Aubry *et al.*Aubry, Sator, and Guillot (2013) for the system {silicate + CO2 + NG} and set the interactions between the NG and the carbon atom of the carbonate anion to zero, as the carbon is considered to be screened by the surrounding oxygen atoms Desmaele et al. (2019a, b) (in fact, the carbon atom is deeply embedded into the oxygen electronic clouds). The choice of such a force field is based on the assumption that the interaction potentials between NGs and melts cited above are transferable from a melt to another (e.g. \ceCO2-bearing silicate to carbonate). It is justified a posteriori by the agreement between the calculated solubilities and the solubility data of Burnard *et al.*Burnard, Toplis, and Medynski (2010) (see Section III).

II.2 Solubility calculations

II.2.1 Test Particle Method (TPM)

The solubility, expressed as a molar fraction, of a noble gas in a melt is given by (see Appendix B.1):

[TABLE]

where $\rho_{m}$ and $\rho_{v}$ are the number densities of the pure melt and of the vapor phase (noble gas), respectively, and $\Gamma_{m}$ and $\Gamma_{v}$ are the solubility constants of the noble gas in the melt and in its pure phase, respectively.

At low pressure $P$, the gas phase can be considered as ideal, allowing to retrieve Henry’s law from Eq. (1):

[TABLE]

where $k_{H}=\rho_{m}k_{B}T/\Gamma_{m}$ is the Henry constant. The solubility of noble gases in the melt was calculated by the Test Particle Method (TPM) introduced by Widom.Widom (1963) This method allows to determine the excess chemical potential of a solute in a solvent, as a function of the potential energy distribution seen by the solute inserted in the solvent (at infinite dilution), namely

[TABLE]

where $\mu^{ex,i}$ is the excess chemical potential of a solute $i$ (e.g. a noble gas), $k_{B}$ is the Boltzmann constant, $T$ is the temperature, $\Psi_{i}$ is the interaction energy between the solute $i$ and the solvent (e.g. a carbonate melt) in which it is solvated, and $\langle\cdots\rangle_{0}$ means that a statistical averaging is done on the solvent molecules only, at a given ($T,P$) condition. Notice that the solute particle acts as a ghost particle, and that many insertions of the latter particle in the solvent configurations are needed to accurately evaluate $\mu^{ex}$ (for further details see Guillot and Sator (2012) and Appendix B.2).

In practice a set of microscopic configurations of the melt is first generated by MD. Then the insertion of the solute particle (the noble gas atom) is attempted many times into the MD-generated configuration. For each insertion of the solute particle the interaction energy $\Psi$ is calculated.

The solubility parameter $\Gamma$ is then evaluated by the averaging over all the attempted insertions

[TABLE]

where $N_{\text{test}}$ is the number of attempted insertions. According to this method, the solubility parameter $\Gamma$ can be evaluated in the two coexisting phases: the carbonate melt and the noble gas phase (providing $\Gamma_{m}$ and $\Gamma_{v}$), and the solubility of the NG can be easily obtained from Eq. (1).

II.3 Explicit Interface Method (EIM)

As an alternative to the TPM, the solubility of a gas in a melt can be calculated from a numerical experiment which consists in simulating explicitly the equilibrium between the two phases (for further details see Guillot and Sator (2011, 2012)). The solubility is then simply obtained by counting the average number of noble gas atoms in the melt, once equilibrium is reached between the two phases, i.e. when the liquid is saturated in gas (Figure 1).

II.4 Complementarity between the two methods

Depending on the thermodynamic conditions of interest, one of the two methods presented above is the most appropriate one to calculate the solubility. Thus in the framework of the TPM, the solubility parameter $\Gamma_{m}$ is theoretically defined at infinite dilution, which means that it is only meaningful for low concentrations of NG in the melt. On the contrary the explicit interface method (EIM) becomes useful when the pressure of the gas phase is sufficiently high for the simulated melt to accommodate a significant number of NG atoms. In brief, the TPM will be used at low pressures and the EIM at high pressures.

In this context we have checked that the two methods give results that are consistent with one another. For example in natrocarbonatite (Na1.1K0.18Ca0.36CO3) at 1600 K and 0.1 GPa: $X^{TPM}_{\text{Ne}}=0.78\pm 0.01$ mol% (using Eq. (1) with $\Gamma_{m}=16.66\pm 0.04\times 10^{-3}$, $\Gamma_{v}=820.5\pm 0.1\times 10^{-3}$, $\rho_{m}=17.98$ mol/L and $\rho_{v}=6.947$ mol/L) and $X^{EIM}_{\text{Ne}}=0.84\pm 0.10$ mol%. Note that the rather large uncertainty on $X^{EIM}_{\text{Ne}}$ (see Figure 1) is due to the relatively low pressure (0.1 GPa) considered here.

III Solubility

III.1 Low pressure

First we focus on the mixtures studied by Burnard *et al.*Burnard, Toplis, and Medynski (2010), namely \ceK2CO3–\ceNa2CO3 mixtures at 1173 K and 1 bar. The solubilities calculated by MD by the means of the TPM, and the ones measured by Burnard *et al.*Burnard, Toplis, and Medynski (2010) are reported on Figure 2. For some mixtures Burnard *et al.*Burnard, Toplis, and Medynski (2010) have made several measurements of the solubility of helium. The error bars on each measurement is low and comprised within the symbol of Figure 2. However the data dispersion for the He solubility at a given melt composition suggests that the uncertainty on the experimental results is actually greater than 50 %. According to the authors, this could be due to the fast diffusion of the He gas out of the glass upon quenching. This He loss is likely systematic although its magnitude is varying. So the solubility data given by Burnard *et al.*Burnard, Toplis, and Medynski (2010) for He should be considered as a lower limit. Hence the results of our simulations for He are compatible with the measurements of Burnard *et al.*Burnard, Toplis, and Medynski (2010). In contrast, for Ar the agreement between the TPM results and the experimental data is better. This is consistent with the assumption made for He (gas loss), because Ar being a heavier element it less easily escapes out of the quenched sample.

Beside He and Ar, we report the solubility for Ne and Xe. For a given composition of the liquid, the solubility increases as the atomic radius of the NG decreases. This behavior is similar to the one observed in silicates, where the solubility is related to the ionic porosity Guillot and Sator (2012). Such a solvation mechanism can be termed as entropic, as it is mainly related to the structural disorder and cavity formation in the liquid. In contrast, in polar solvents (e.g. water) the solvation of NG is enthalpically driven with a solubility that increases with the polarizability of the solute Smith and Kennedy (1983); Guillot and Guissani (1993) (i.e. with the size of the NG).

In their study, Burnard *et al.*Burnard, Toplis, and Medynski (2010) reported that the solubility is hardly sensitive to the composition of the melt. We believe that the high uncertainties on their data as well as the narrowness of the studied composition range ($0.2<x_{\ce{K2CO3}}\leq 0.5$) are accountable for this observation. On the contrary MD simulations reveal that the composition of the melt has a strong effect (especially for $x_{\ce{K2CO3}}<0.25$ and $x_{\ce{K2CO3}}>0.4$), regardless of the noble gas species that is considered. In fact the solubility decreases continuously when the content of calcite increases. For instance the solubility of He contrasts by an order of magnitude between the two end-members \ceK2CO3 end \ceCaCO3 (see Figure 2).

III.2 Evolution with pressure

Simulating a biphasic system (EI method) enables to calculate the NG solubility under high pressures. We focus on two carbonate melts: molten dolomite (\ceCa_0.5Mg_0.5CO3) at 1673 K and a \ceNa2CO3–\ceK2CO3–\ceCaCO3 mixture (in proportions 55, 9 and 36 mol%, respectively) at 1600 K, modeling the natrocarbonatite emitted at Ol Doinyo Lengai.Keller and Zaitsev (2012); Desmaele et al. (2019b) Figure 3 reports the solubility of He, Ne, Ar and Xe calculated at pressures ranging from 3 to 6 GPa for dolomite, and from 0.1 to 6 GPa for the natrocarbonatite. The solubilities in the natrocarbonatite were also calculated at 1 bar using the TPM.

Over the whole studied pressure range, the solubility of a given NG is greater in the natrocarbonatite than in molten dolomite. This is consistent with the decrease of the solubility when increasing the molar fraction of alkaline-earth cation as observed at low pressure in \ceK2CO3-\ceCaCO3 mixtures (Figure 2). The observation that the solubility is negatively correlated with the size of the NG still stands at high pressure (at least up to 3 GPa). However beyond 3 GPa, the solubilities of Ne and He come close to each other in natrocarbonatite and even tend to cross at about 6 GPa (He shows a solubility maximum at about 2 GPa), features which are also observed in dolomite. In the same way, the solubility of Xe in natrocarbonatite increases steadily and becomes as high as the one of Ar at 6 GPa, the latter leveling off at about 1 GPa.

In the natrocarbonatite melt, the solubility of NGs first increases quasi linearly with pressure, so Henry’s law is fulfilled up to a few kbar (Figure 3). Above $\sim 1$ GPa the solubility generally levels off and it even goes through a maximum for He at $P\sim 2$ GPa and for Ar at $P\sim 4$ GPa. As for Xe, its solubility is quasi linear with $P$ (up to $\sim 2$ GPa). In molten dolomite, the solubility of noble gases barely depends on pressure between 3 and 6 GPa. In comparison the solubility of noble gases in silicate melts increases up to $\sim 5$ GPa. At higher pressures, some experimental studies report a drop of the solubilityBouhifd, Jephcoat, and Kelley (2008) and other do not.Schmidt and Keppler (2002); Niwa et al. (2013) In their simulation study Guillot and Sator (2012) show a good agreement with the experimental results up to $\sim 5$ GPa and predict a plateau value for the solubility with a slowly decreasing behavior at higher pressures. This trend is similar to what we observe in molten carbonates (Figure 3), although the solubility plateau occurs at a lower pressure for carbonates. In any case at the pressures of the Earth’s upper mantle, our results point out that the solubility of the nobles gases has the same order of magnitude in molten carbonates and in molten silicates Guillot and Sator (2012). This suggests that noble gases would not partition massively in a carbonated phase. However, seeing the significant dependence of the NG solubilities with the composition of the melt (a factor of $\sim 3$ for Ar in dolomitic versus natrocarbonatitic melts, see the inset of Figure 3) the investigation of the partition coefficients between carbonate and silicate deserves a further studying. Still, considering that carbonatite melts only represent a minor fraction of magmatic liquids, they can unlikely be the main carrier of noble gases.

IV Surface Tension

To our knowledge, there are no data of the pressure evolution of the surface tension at the interface between a noble gas phase and molten carbonates, despite the interest in a geochemical perspective. For example if we consider a gas-saturated magmatic melt, it is the liquid/gas surface tension that controls the formation, the growth and the coalescence of gas bubbles below a certain supersaturation pressure Mangan and Sisson (2005).

By simulating a biphasic system in MD as shown in Figure 1 for the calculation of the solubility by the EIM, the surface tension $\gamma$ can be calculated simultaneously for the same computational cost. It is given by the long time limit of the average of the diagonal components of the stress tensor:Kirkwood and Buff (1949)

[TABLE]

where $L_{z}$ is the length of the simulation box perpendicular to the interface and $\Pi_{xx}$, $\Pi_{yy}$ and $\Pi_{zz}$ are the diagonal components of the stress tensor evaluated in the entire simulation box (see Desmaele et al. (2019a) for more details).

Figure 4 plots the surface tension between the natrocarbonatitic (or dolomitic) melt and the noble gas phase (He, Ne, Ar and Xe), calculated at different pressures. For natrocarbonatite at 1600 K the surface tension was also calculated at the interface with a vacuum ($190~{}\pm~{}3$ mN/m). This value compares well with the ones obtained at the interface with noble gases at 0.1 GPa ($171-206$ mN/m).

Irrespective of the nature of the melt and of the noble gas considered, the surface tension increases with increasing pressure. In contrast, it is known that the surface tension of silicate melts in contact with \ceH2O decreases with pressure.Colucci, Battaglia, and Trigila (2016) The behavior of carbonate melts may somehow appear as counter-intuitive because with increasing pressure the density of the gas phase approaches that of the melt, thus the energetic cost to create the interface should decrease. However the surface tension involves subtle mechanisms which are likely pressure-dependent. Among other factors, the amount of gas solubilized in the melt which increases with pressure (see previous section) could be responsible for a non-negligible contribution (positive or negative) to the surface stress. On the one hand, it could tend to destabilize the melt (and thus decrease $\gamma$), but on the other hand the gas phase in contact with the melt can act like a piston and rigidify the interface (and thus increase $\gamma$). In any case, the solubilized gas is likely not the only factor in action as the rather large differences of solubility observed between the four gases (almost a factor of ten between He and Xe, see Figure 3) are not retrieved for the surface tension (see Figure 4).

There is also a distinct effect of the melt composition on the surface tension: For a given noble gas species and at a given pressure, the surface tension of dolomite is greater than the one of the natrocarbonatite by $\sim$ 50%. This increase of the surface tension with the amount of alkaline-earth cations in the melt is also observed at the interface with a vacuumDesmaele (2017) (compare the values of $\gamma$ for the natrocarbonatite and for the \ceCaCO3 melt on Figure 4).

As for the effect of the gas composition, it is trickier to decipher. At a given pressure and for a given melt composition, there is no systematic trend of the surface tension as a function of the size of the noble gas atom (except for He in dolomite, see Figure 4). Moreover, when pressure is increased the hierarchy between noble gases seems to modify somewhat. It is possible that the uncertainties on the calculated values are greater than the ones we estimate from the fluctuations ($\sim 5\%$, see Figure C.1 in the appendix) of the running average of Eq. (5).

V Conclusion

To complete and go beyond the precursory study of Burnard *et al.*Burnard, Toplis, and Medynski (2010), the solubility of noble gases (He, Ne, Ar and Xe) in carbonate melts was calculated by molecular dynamics simulations. These simulations used empirical interaction potentials whose accuracy was previously demonstrated in studying the thermodynamic and transport properties of carbonate melts and the solubility of noble gases in silicate melt.Desmaele (2017); Desmaele et al. (2019a, b); Guillot and Sator (2012)

The NG solubilities were first calculated in \ceK2CO3-\ceCaCO3 mixtures at 1 bar and the results are in a fair agreement with the data of Burnard *et al.*Burnard, Toplis, and Medynski (2010) once considered the uncertainties on the experimental values. Then we investigated the effect of pressure (up to 6 GPa), focusing on two melt compositions: a dolomitic one and a mixture modeling the carbonatitic lava emitted at Ol Doinyo Lengai. We observed that the solubility decreases with the amount of alkaline-earth cation in the melt and with the size of the noble gas (entropy-driven solubility). Concerning the solubility in the natrocarbonatitic melt, Henry’s law is fulfilled at low pressures ($\sim 0.1$ GPa). At higher pressures (a few GPa) the solubility levels off or even starts to diminish smoothly (for He at $P>2$ GPa and Ar at $P>4$ GPa). In contrast, in molten dolomite the effect of pressure is negligible on the studied $P$ range ($3-6$ GPa). At the pressures of the Earth’s mantle, the solubilities of noble gases in carbonate melts are still of the same order of magnitude as the ones in molten silicates, a finding in agreement with the ratio 4He/40Ar measured in the gas emitted from Ol Doinyo Lengai crater in 2005, which is close to the mantle value (see Fischer et al. (2009)). Furthermore this finding also suggests that carbonatitic melts at depth cannot be preferential carriers of noble gases.

Finally we provided some insight into the surface tension at the interface between carbonate melts and noble gases. With increasing pressure ($P$ from 0 to 6 GPa), the surface tension increases (by a factor $\sim 2$), whatever the composition of the melt and of the NG phase. This is in strong contrast with the effect of \ceH2O on the surface tension of silicate melts which drops when pressure increases to a few GPa.

Acknowledgements.

The research leading to these results has received funding from the Région Ile-de-France and the European Community’s Seventh Framework Program (FP7/2007-2013) under Grant agreement (ERC, N∘ 279790).

Appendix A Force Field

The interactions between two noble gas (NG) atoms $i$ and $j$ were modeled by the FF used by Guillot and Sator (2012) and consisting of Buckingham potentials:

[TABLE]

Values of the parameters for these interactions are collected in Table A.1.

For the interactions within the carbonate melts we used the FF developed in Desmaele et al. (2019a, b). The intramolecular potential energy associated with a carbonate molecule anion consists of an oxygen-oxygen repulsive potential as in Eq. (A.1) and of a carbon-oxygen potential (harmonic stretching + coulombic interaction): $V_{\text{CO}}^{intra}(r_{\text{CO}})=\frac{1}{2}{k_{\text{CO}}}(r_{\text{CO}}-{r_{0,\text{CO}}})^{2}+{q_{\text{O}}q_{\text{C}}}/{4\pi\epsilon_{0}r_{\text{CO}}}$. Moreover, two elements $i$ and $j$, with $i,j$= Na, K, Ca, Mg, O and C (with O and C not belonging to a same carbonate group) interact through a pair potential: $V_{ij}(r_{ij})={A_{ij}}\exp(-{r}_{ij}/{\rho_{ij}})-C_{ij}/r_{ij}^{6}+{q_{i}q_{j}}/{4\pi\epsilon_{0}{r}_{ij}}$, that is a sum of a van der Waals (Buckingham-like) and of a coulombic term. All the parameters for the melt-melt interactions are summarized in Table A.2.

For the NG-melt interactions we used Lennard-Jones potentials:

[TABLE]

as adjusted by Guillot and Sator (2012) and by Aubry *et al.*Aubry, Sator, and Guillot (2013) (Table A.3).

Appendix B The Test Particle Method

B.1 Thermodynamics

Calculating the solubility of a noble gas within a melt consists in considering a thermodynamic equilibrium between the gas and the melt (in contact with each other), expressed by the equality of the chemical potentials of the gas in both phases:

[TABLE]

The chemical potentials can be decomposed as the sum $\mu_{v,m}=\mu_{v,m}^{IG}+\mu_{v,m}^{ex}$, of an ideal gas part

[TABLE]

where $\rho^{g}$ is the number density of the noble gas phase in its pure phase ($\rho_{v}^{g}$) or in the melt ($\rho_{m}^{g}$) and $C$ is a constant that is equal in both phases (gas and melt), and an excess part $\mu^{ex}$ ($\mu_{v}^{ex}$ in the gas phase, $\mu_{m}^{ex}$ in the melt). Introducing Eq. (B.4) into Eq. (B.3) leads to:

[TABLE]

where $\rho_{m}^{g}$ and $\rho_{v}^{g}$ are the number densities of the gas in the melt and in the vapor phase, respectively, and $\Gamma_{m}$ and $\Gamma_{v}$ are the solubility parameters of the gas in the two phases.

The solubility of the gas in the melt, expressed as the number of moles of gas in the melt divided by the number of moles of melt and of noble gas, is given by

[TABLE]

where $\rho_{v}$ is the number density of the gas phase and $\rho_{m}$ is the number density of the melt.

At low pressure, the gas phase can be considered as ideal: $\mu_{v}^{ex}=0$ (i.e. $\Gamma_{v}=1$ and $P~{}=~{}\rho_{v}k_{B}T$) thus $\rho_{m}^{g}=\Gamma_{m}\rho_{v}^{g}=\Gamma_{m}P/k_{B}T$. Then from equation (B.6) it follows:

[TABLE]

which is nothing but Henry’s law. In this case, only the solubility parameter of the noble gas in the melt, $\Gamma_{m}$, has to be calculated.

B.2 Simulation details

Classical MD simulations of the carbonate melt ($N\simeq$ 1000 atoms) were carried out using the DL_POLY 2 software Smith and Forester (1996), with a timestep of 1 fs. The simulations were performed in the $NVE$ ensemble with an equilibration run of 0.5 ns, followed by a production run of at least 10 ns. Insertion of the test particle in the simulation box was attempted every 1000 timesteps (1 ps) on a grid of width $\rm{d}L$, with $L/\rm{d}L=10$ to 50 (i.e. the grid has $10^{3}=1000$ to $50^{3}=125000$ nodes, Figure B.2 evidences that convergence is reached for $L/\rm{d}L\geq 30$). In practice, if the distance between the node and a solvent atom is too small (within a defined cutoff radius $r_{min}$, see Figure B.1), the insertion is rejected and a null contribution is added to the average of Eq. (4) (because the value of $\mathrm{e}^{-\Psi/k_{B}T}$ in Eq. (4) is vanishingly small for this event). Note that when the density of the solvent is high or when the inserted particle has a large van der Waals radius, the rate of rejection for the insertion is high and uncertainties may become important(e.g. Ne, Ar and Xe in \ceK2CO3, see Figure 2). In this study, the solubility calculations using this method have uncertainties of a few percents. For more detailed applications with MD, see Guillot and Sator (2012) and Aubry *et al.*Aubry, Sator, and Guillot (2013).

B.3 Data

Appendix C The Explicit Interface Method

C.1 Simulation details

The EIM simulations consisted in modeling the solvent ($N\simeq$ 2000 atoms of melt) in contact with a NG phase in a parallelepipedic simulation box. The number of NG atoms in the box ($N\simeq$ 300 - 900 atoms) was chosen so that the vapor phase was large enough for the consecutive periodic images of the melt to not interact with one another.

An equilibration run was first performed in the $NPT$ ensemble (using a Nosé-Hoover thermostat) for 0.9 ns (including 0.5 ns to equilibrate the temperature) allowing to reach an accuracy on the density value of $\Delta n/n\sim\pm 1\%$ for the two coexisting phases. To evaluate the solubility $X$, simulations were performed in the $NVE$ ensemble with an equilibration run of 0.5 ns, followed by a production run of 10 ns. Configurations were extracted every 1 ps to determine the density profiles and calculate the surface tension between the NG fluids and the carbonate melts (Figure C.1).

C.2 Data