The MgCO$_3$-CaCO$_3$-Li$_2$CO$_3$-Na$_2$CO$_3$-K$_2$CO$_3$ Carbonate Melts: Thermodynamics and Transport Properties by Atomistic Simulations

TL;DR

This study uses atomistic simulations to investigate the thermodynamics and transport properties of complex carbonate melts relevant to Earth's mantle, extending force fields to include Mg and Ca components and validating against experimental data.

Contribution

The paper introduces a new force field for Mg- and Ca-bearing carbonate melts and demonstrates its accuracy in modeling their thermodynamic and transport properties across wide temperature and pressure ranges.

Findings

Force field accurately models MgCO3 and CaCO3 melts.

Simulation results agree well with experimental data.

First study to explore transport properties of magnesite and dolomite melts.

Abstract

Atomistic simulations provide a meaningful way to determine the physico-chemical properties of liquids in a consistent theoretical framework. This approach takes on particular usefulness for the study of molten carbonates, in a context where thermodynamic and transport data are crucially needed over a large domain of temperatures and pressures (to ascertain the role of these melts in geochemical processes) but are very scarce in the literature, especially for the calco-magnesian compositions prevailing in the Earth's mantle. Following our work on Li2CO3-Na2CO3-K2CO3 melts, we extend our force field to incorporate Ca and Mg components. The empirical interaction potentials are benchmarked on the density data available in the experimental literature (for the crystals and the K2Ca(CO3)2 melt) and on the liquid structure issued from ab initio molecular dynamics simulations. Molecular dynamics simulations are then performed to study the thermodynamics, the microscopic structure, the diffusion coefficients, the electrical conductivity and the viscosity of molten Ca, Mg-bearing carbonates up to 2073 K and 15 GPa. Additionally, the equation of state of a Na-Ca-K mixture representative of the lavas emitted at Ol Doinyo Lengai (Tanzania) is evaluated. The overall agreement between the MD results and the existing experimental data is very satisfying and provides evidence for the ability of the force field to accurately model any MgCO3-CaCO3-Li2CO3-Na2CO3-K2CO3 melt over a large T - P range. Moreover it is the first report of a force field allowing to study the transport properties of molten magnesite (MgCO3) and molten dolomite (CaMg(CO3)2).

Tables

Table 1. Table 1: Intermolecular Buckingham parameters and partial charges. Note that for electroneutrality considerations q C = + 1.04085 subscript 𝑞 C 1.04085 q_{{\rm{C}}}=+1.04085  e. Intramolecular Born repulsion parameters between oxygen atoms of a same carbonate ion 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  Å.

$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
Li	O	300 000	0.2228	1 210	+0.82101
Na	O	1 100 000	0.2228	3 000	+0.82101
K	O	900 000	0.2570	7 000	+0.82101
O	O	500 000	0.252525	2 300	$-$0.89429

Table 2. Table 2: Density and bulk modulus of \ce MgCO3, \ce CaCO3 and \ce CaCO3-bearing mixtures at T = 1100 𝑇 1100 T=1100  K calculated either by MD simulations ( n M ​ D 0 superscript subscript 𝑛 𝑀 𝐷 0 n_{MD}^{0} , K M ​ D 0 superscript subscript 𝐾 𝑀 𝐷 0 K_{MD}^{0} ) or by using ideal mixing rules ( n M ​ D , m ​ i ​ x 0 superscript subscript 𝑛 𝑀 𝐷 𝑚 𝑖 𝑥 0 n_{MD,\,mix}^{0} , K M ​ D , m ​ i ​ x 0 superscript subscript 𝐾 𝑀 𝐷 𝑚 𝑖 𝑥 0 K_{MD,\,mix}^{0} ). Extrapolation from experiments n E ​ x ​ p , m ​ i ​ x 0 superscript subscript 𝑛 𝐸 𝑥 𝑝 𝑚 𝑖 𝑥 0 n_{Exp,\,mix}^{0} , K E ​ x ​ p , m ​ i ​ x 0 superscript subscript 𝐾 𝐸 𝑥 𝑝 𝑚 𝑖 𝑥 0 K_{Exp,\,mix}^{0} are given for comparison. Liu and Lange ( 2003 ); O’Leary, Lange, and Ai ( 2015 ) Note that for \ce MgCO3 an additional assumption relative to the cationic mass is needed to estimate n E ​ x ​ p , m ​ i ​ x 0 superscript subscript 𝑛 𝐸 𝑥 𝑝 𝑚 𝑖 𝑥 0 n_{Exp,\,mix}^{0} (see Hurt and Lange ( 2019 ) ). In parenthesis are given the MD results from the models of Vuilleumier  et al. ( 2014 ) ( † ) and Hurt and Wolf ( 2018 ) ( ‡ ).

	\ceMgCO3	\ceCaCO3	\ceCaMg(CO3)2	\ceK2Ca(CO3)2	Natro
$n_{MD}^0$ (g/cm3)	2.49	2.44 (2.47†, 2.44‡)	2.49	2.14	2.11
$n_{MD,mix}^0$(g/cm3)	–	–	2.46	2.12	2.11
$n_{Exp,mix}^0$ (g/cm3)	2.45	2.49	2.47	2.13	2.13
$K_{MD}^0$ (GPa)	15.0	17.6 (15.7†, 15.1‡)	12.8	9.1	11.4
$K_{MD,mix}^0$ (GPa)	–	–	–	8.6	10.4
$K_{Exp,mix}^0$ (GPa)	–	18.7	–	8.9	10.4

Table 3. Table 3: Parameters of the third-order Birch-Murnaghan equation of state and T − P 𝑇 𝑃 T-P domain of validity. Natro refers to the natrocarbonatite melt of composition Na 1.1 K 0.18 Ca 0.36 CO 3 (see text). The reference temperature T r ​ e ​ f subscript 𝑇 𝑟 𝑒 𝑓 T_{ref} is defined for numerical purposes and set to the melting temperature for each composition, n r ​ e ​ f 0 superscript subscript 𝑛 𝑟 𝑒 𝑓 0 n_{ref}^{0} and K r ​ e ​ f 0 superscript subscript 𝐾 𝑟 𝑒 𝑓 0 K_{ref}^{0} are the density and bulk modulus at 1-bar and T r ​ e ​ f subscript 𝑇 𝑟 𝑒 𝑓 T_{ref} . Other empirical parameters ( α 0 subscript 𝛼 0 \alpha_{0} , α 1 subscript 𝛼 1 \alpha_{1} , b 1 subscript 𝑏 1 b_{1} and b 2 subscript 𝑏 2 b_{2} ) are defined by n T 0 = n r ​ e ​ f 0 ​ e ∫ T r ​ e ​ f T − ( α 0 + α 1 ​ T ) ​ d ​ T superscript subscript 𝑛 𝑇 0 superscript subscript 𝑛 𝑟 𝑒 𝑓 0 superscript 𝑒 superscript subscript subscript 𝑇 𝑟 𝑒 𝑓 𝑇 subscript 𝛼 0 subscript 𝛼 1 𝑇 d 𝑇 n_{T}^{0}=n_{ref}^{0}e^{\int_{T_{ref}}^{T}-(\alpha_{0}+\alpha_{1}T)\,\mathrm{d}T} and K T 0 = K r ​ e ​ f 0 / ( 1 + b 1 ​ ( T − T r ​ e ​ f ) + b 2 ​ ( T − T r ​ e ​ f ) 2 ) superscript subscript 𝐾 𝑇 0 superscript subscript 𝐾 𝑟 𝑒 𝑓 0 1 subscript 𝑏 1 𝑇 subscript 𝑇 𝑟 𝑒 𝑓 subscript 𝑏 2 superscript 𝑇 subscript 𝑇 𝑟 𝑒 𝑓 2 K_{T}^{0}={K_{ref}^{0}}/\big{(}1+b_{1}(T-T_{ref})+b_{2}(T-T_{ref})^{2}\big{)} . See Desmaele  et al. ( 2019 ) for details on the fitting method.

	\ceMgCO3	\ceCaCO3	\ceCaMg(CO3)2	Natro
$T$ (K)	$1873-2073$	$1100-2073$	$1653-2073$	$823-1600$
$P$  (GPa)	$0-15$	$0-15$	$0-15$	$0-6$
$T_{ref}$ (K)	1873	1623	1653	1073
$n_{ref}$ (${\text{g/cm}}^3$)	2.23	2.25	2.26	2.14
${\alpha }_0$ (K-1)	$-1.50\times {10}^{-4}$	$-1.54\times {10}^{-4}$	$-1.23\times {10}^{-4}$	$-2.36\times {10}^{-4}$
${\alpha }_1$ (K-2)	$8.22\times {10}^{-8}$	$0.49\times {10}^{-8}$	$-7.65\times {10}^{-8}$	$-0.20\times {10}^{-8}$
$K_{ref}^0$(GPa)	10.98	12.74	12.13	11.72
$b_1$ (K-1)	$5.0\times {10}^{-4}$	$6.1\times {10}^{-4}$	$3.7\times {10}^{-4}$	$10.2\times {10}^{-4}$
$b_2$ (K-2)	$2.0\times {10}^{-7}$	$1.5\times {10}^{-7}$	$5.0\times {10}^{-7}$	$7.7\times {10}^{-7}$
$K^{\prime 0}$	9.5	7.7	8.5	8.0

Table 4. Table 4: Parameters of the Arrhenius activation law ( 7 ) obtained by the interpolation of all MD simulation points. P 𝑃 P is in GPa.

	\ceMgCO3	\ceCaCO3	\ceCaMg(CO3)2 (left: X=Mg, right: X=Ca)
$D_{0,\text{X}}$ (${10}^{-9}$m2/s)	117	122	164	150
$E_a^{D_{\text{X}}}$ (kJ/mol)	54	54	58	57
$V_a^{D_{\text{X}}}$ (cm3/mol)	$2.5-0.07P+0.0016P^2$	$4.6-0.28P+0.01006P^2$	$3.5-0.13P+0.0029P^2$	$4.6-0.28P+0.0100P^2$
$D_{0,{\text{CO}}_3}$ (${10}^{-9}$m2/s)	83	86	100
$E_a^{D_{{\text{CO}}_3}}$(kJ/mol)	54	51	54
$V_a^{D_{{\text{CO}}_3}}$(cm3/mol)	$2.6-0.11P+0.0036P^2$	$3.3-0.08P+0.0007P^2$	$2.1-0.08P+0.0085P^2$

Table 5. Table 5: Parameters of the Arrhenius activation laws ( 8 ) and ( 9 ) for the electrical conductivity and the viscosity, obtained by the interpolation of all MD simulation points. P 𝑃 P is in GPa.

	\ceMgCO3	\ceCaCO3	\ceCaMg(CO3)2
${\sigma }_0$ (S/m)	3842.5	2593.7	4838.7
$E_a^{\sigma }$(kJ/mol)	38	34	42
$V_a^{\sigma }$ (cm3/mol)	1.0+0.076P-0.0036P2	2.3+0.048P-0.0041P2	2.7+0.090P-0.0043P2
${\eta }_0$ (Pa$\cdot$s)	3.9$\times$10-4	2.8$\times$10-4	1.3$\times$10-4
$E_a^{\eta }$ (kJ/mol)	37	39	50
$V_a^{\eta }$ (cm3/mol)	2.1-0.02P	3.6-0.09P	3.4-0.05P

Table 6. Table 6: Haven ratio and its ionic contributions, where X denotes all cations. Note that the error bar on these quantities is of the order of the error bar on the electrical conductivity: ∼ 5 − 10 similar-to absent 5 10 \sim 5-10 %. For mixtures, the notation X 1 (= Mg, Na or K) and X 2 (= Ca or Mg) is the same as the one appearing in the composition name (first column, with M = = Mg, C = = Ca, N = = Na and K = = K).

	$T$(K)	$P$ (GPa)	$\sigma$ (S/m)	$H$	$H_{X_1}$	$H_{X_2}$	$H_{{\mathrm{CO}}_3}$	$H_{X_1-X_2}$	$H_{X-{\mathrm{CO}}_3}$
\ceMgCO3	1823	3.0	245	0.92	$-0.14$	–	$-0.32$	–	0.38
\ceCaCO3	1773	1.0	220	0.81	$-0.22$	–	$-0.35$	–	0.38
MC	1573	3.0	112	0.85	$-0.08$	0	$-0.36$	$-0.03$	0.34
NM	1573	3.0	186	0.87	$-0.04$	$-0.04$	$-0.29$	$-0.16$	0.40
KM	1573	3.0	75	0.72	$-0.10$	$-0.03$	$-0.23$	$-0.22$	0.29
KC	1573	3.0	95	0.75	$-0.10$	$-0.05$	$-0.25$	$-0.20$	0.35

Equations

$$n_{mix}=\frac{{\sum }_i{x}_i{M}_i}{{\sum }_i{x}_i{\bar{V}}_i},$$

$$\begin{split}P=&\cfrac{3}{2}K^{0}(T)\Bigg{(}\Big{(}\cfrac{n}{n_{T}^{0}}\Big{)}^{7/3}-\Big{(}\cfrac{n}{n_{T}^{0}}\Big{)}^{5/3}\Bigg{)}\\ &\times\left(1-\cfrac{3}{4}(4-K^{\prime 0})\Bigg{(}\Big{(}\cfrac{n}{n_{T}^{0}}\Big{)}^{2/3}-1\Bigg{)}\right)\enspace,\end{split}$$

$$K_{mix}=\frac{{\sum }_i{x}_i{\bar{V}}_i}{{\sum }_i{x}_i{\bar{V}}_i/K_i}.$$

$$D_s$$

$$\sigma$$

$$\eta$$

$$D_s(P,T)$$

$$\sigma (P,T)$$

$$\eta (P,T)$$

$${\sigma }^{NE}=\frac{e^2}{k_{\mathrm{B}}TV}\sum _s{N}_s{z}_s^2{D}_s,$$

$$H=\frac{\sigma }{{\sigma }^{NE}}=1+\sum _s{H}_s+\sum _s\sum _{s^{\prime }\ne s}{H}_{s{s}^{\prime }},$$

$$\begin{array}{rl}{H}_s=\underset{t\to \infty }{\lim }& \frac{1}{6t}\frac{z_s^2}{{\sum }_s{N}_s{z}_s^2{D}_s}\\ & \times \sum _{i=1}^{N_s}\sum _{j\ne i}^{N_s}⟨{\vec{\Delta }}_i^{(s)}(t)\cdot {\vec{\Delta }}_j^{(s)}(t)⟩,\\ \mathrm{where}& {\vec{\Delta }}_i^{(s)}(t)={\vec{r}}_i^{(s)}(t)-{\vec{r}}_i^{(s)}(0)\end{array}$$

$$\begin{array}{rl}{H}_{s{s}^{\prime }}=\underset{t\to \infty }{\lim }& \frac{1}{6t}\frac{z_s{z}_{s^{\prime }}}{{\sum }_s{N}_s{z}_s^2{D}_s}\\ & \times \sum _{i=1}^{N_s}\sum _{j=1}^{N_{s^{\prime }}}⟨{\vec{\Delta }}_i^{(s)}(t)\cdot {\vec{\Delta }}_j^{(s^{\prime })}(t)⟩.\end{array}$$

$${\eta }_{SE}=\frac{k_{\mathrm{B}}T}{2\pi Dd},$$

$$\sigma =\frac{A}{{\eta }^{0.8}},$$

Full text

Following our work on \ceLi2CO3–\ceNa2CO3–\ceK2CO3 melts,Desmaele et al. (2019); Desmaele (2017) we extend our force field to incorporate Ca and Mg components. The empirical interaction potentials are benchmarked on the density data available in the experimental literature (for the crystals and the \ceK2Ca(CO3)2 melt) and on the liquid structure issued from ab initio molecular dynamics simulations. Molecular dynamics simulations are then performed to study the thermodynamics, the microscopic structure, the diffusion coefficients, the electrical conductivity and the viscosity of molten Ca, Mg-bearing carbonates up to 2073 K and 15 GPa. Additionally, the equation of state of a Na–Ca–K mixture representative of the lavas emitted at Ol Doinyo Lengai (Tanzania) is evaluated.

The overall agreement between the MD results and the existing experimental data is very satisfying and provides evidence for the ability of the force field to accurately model any \ceMgCO3–\ceCaCO3–\ceLi2CO3–\ceNa2CO3–\ceK2CO3 melt over a large $T-P$ range. Moreover it is the first report of a force field allowing to study the transport properties of molten magnesite (\ceMgCO3) and molten dolomite (\ceCaMg(CO3)2).

molten carbonates, molecular dynamics, equation of state, microscopic structure, transport properties, diffusion coefficients, electrical conductivity, viscosity, carbonatite

I Introduction

Whether it is for the applied or the fundamental fields in which carbonate melts are implied (electrochemistry, geochemistry), having a reliable model for their alkali, alkaline-earth mixtures (and their end-members) is essential. For instance the longevity of molten carbonate fuel cells, usually based on alkali eutectic mixtures, is known to be improved by the addition of alkaline-earth cations,Lair et al. (2012) in a small enough amount that doesn’t compromise the efficiency of the cell.Chery, Lair, and Cassir (2015a, b); Cassir, McPhail, and Moreno (2012); Cassir, Ringuedé, and Lair (2013); Lair et al. (2012) From the viewpoint of geosciences, carbonate melts, although they constitute a very minor phase of the Earth’s mantle (the most part is silicate), are thought to have important implications in the lowering of the melting temperature of silicate rocks, the deep carbon cycle or the high conductivity anomalies observed in the 70–200 km depth region.Dasgupta and Hirschmann (2006); Gaillard et al. (2008); Dasgupta and Hirschmann (2010); Dasgupta (2013); Hammouda and Keshav (2015) On Earth, evidence of a past volcanic activity induced by carbonatitic lava (\ceSiO2 $<20$ wt%) is given by petrology.Mitchell (2005); Woolley and Kjarsgaard (2008) Nowadays, Ol Doinyo Lengai in Tanzania is the only active volcano to produce these remarkable (low temperature, low viscosity) carbonatitic lavas. The carbonatitic melt is mainly composed of a \ceNa2CO3–\ceK2CO3–\ceCaCO3 mixture (in proportions 55:9:36 mol% according to Keller and Zaitsev (2012)), called natrocarbonatite, whereas the majority of other inventoried carbonatites are of calco-magnesian composition.Woolley and Church (2005) In addition, in the Earth’s mantle the most abundant carbonate compositions are calcite, \ceCaCO3, predominating at shallow depth and magnesite, \ceMgCO3, progressively taking over as depth increases. As a consequence their mixture and particularly dolomite, \ceCaMg(CO3)2, are of major interest. For pressure of a few kbars and beyond, these Ca,Mg-carbonates form stable liquid phases over a wide temperature domain.Spivak et al. (2012); Solopova et al. (2013, 2015) But at atmospheric pressure they break down into \ceCO2 + oxide (CaO, MgO) at temperatures below their melting point (decarbonation occurs for $P$ below $\sim$1 GPa for calcite and $\sim$3 GPa for magnesite and dolomite).Irving and Wyllie (1975); Suito et al. (2001) Consequently only a few data are available in the literature for calco-magnesian compositions (most of which are alkali-bearing mixtures stable at low pressure),Gaillard et al. (2008); Kojima (2009); Lair et al. (2012) and very few at high pressures.Dobson et al. (1996); Kono et al. (2014); Sifré, Hashim, and Gaillard (2015) Thus for molten \ceMgCO3 there is simply no data to our knowledge.Hurt and Lange (2019)

In this context where experimental data are sparse in terms of thermodynamic conditions, chemical composition and physical properties, a noticeable advantage of molecular dynamics simulations (classical and ab initio) is that a variety of properties (structure, equation of state and transport coefficients) can be computed in the same theoretical framework. However, one has to be aware that collective quantities such as the ionic conductivity and the viscosity are calculated from slowly converging time correlation functions, that necessarily require (to be accurate) a rather large number of atoms and long time runs. For this reason, ab initio molecular dynamics simulations (AIMD), deriving from electronic density calculations, are either restricted to the study of thermodynamics and structure propertiesZhang and Liu (2015); Li et al. (2017) or only provide crude estimates of transport coefficients, especially when it is question to establish their evolution with pressure and temperature.Vuilleumier et al. (2014); Du et al. (2018) As for classical molecular dynamics simulations (MD), based on a more empirical approach, very few studies have provided estimates of the viscosity and electrical conductivity of carbonate melts, namely \ceCaCO3,Vuilleumier et al. (2014) the \ceLi2CO3-\ceK2CO3 eutectic mixture (62:38 mol%)Corradini, F.-X., and Vuilleumier (2016) and more recently the \ceLi2CO3–\ceNa2CO3–\ceK2CO3 mixtures and end-members.Desmaele et al. (2019) Obviously, the scarcity of thermodynamic dataLiu and Lange (2003); Hudspeth, Sanloup, and Kono (2018); Hurt and Lange (2019) does not help for developing empirical potentials. As a consequence, very few force fields (FF) are available for Ca-Mg and most of them were developed for crystals Yuen, Lister, and Nyburg (1978); Dove et al. (1992); Pavese et al. (1992, 1996); Fisler, Gale, and Cygan (2000); Archer et al. (2003); Raiteri et al. (2010) based on the Born model of solidsBorn and Huang (1954). A first step towards an accurate description of the \ceCaCO3 melt was made by Genge *et al.*Genge, Price, and Jones (1995) by adapting the FF of Dove et al. (1992) (for the crystal). The resulting FF reproduces quite well the liquid structure issued from a AIMD calculation published long afterVuilleumier et al. (2014), but the calculated pressure is overestimated by $\sim+15$ kbar, compared to the recently published equation of state of Zhang and Liu (2015) based on AIMD calculations. Moreover the melting temperature is underestimated. Both these features point toward a too weakly cohesive melt, that may result from the absence of dispersion interaction in the FF.Desmaele et al. (2019) In another attempt Hurt and Wolf (2018) adapted the model of Archer et al. (2003) and proposed a FF fitted on crystalline properties for carbonates in the \ceCaCO3–\ceSrCO3–\ceBaCO3 system.

For \ceMgCO3, in absence of experimental data, Hurt and Wolf (2018) only list in their study some FF parameters, but no explanation on the way they were derived and no result is given for this composition. In the present study, we have developed an empirical force field to describe the thermodynamics and transport properties of \ceMgCO3–\ceCaCO3 melts.

In a previous study devoted to the \ceLi2CO3–\ceNa2CO3–\ceK2CO3 system, we have demonstrated the ability of molecular dynamics (MD) simulations to reproduce experimental data for alkali carbonate melts.Desmaele et al. (2019) Here we extend the latter study to incorporate alkaline-earth components (Mg, Ca). In our study of alkali carbonates we evidenced the need to take into account in the FF the dispersion interactions, which were not included in the FF previously developed by our team for \ceCaCO3.Vuilleumier et al. (2014) Thus for the sake of consistency the interactions between the carbonate anions and the alkali or alkaline-earth cations are treated on the same footing in the present FF. In section II the details of the simulations are presented. In section III the thermodynamic properties (equation of state) and the liquid structure are reported, whereas in section IV we evaluate the transport coefficients (diffusion coefficients, electrical conductivity and viscosity). Furthermore, the heuristicity of the Nernst-Einstein equation (relating the electrical conductivity and diffusion coefficients) and of the Stokes-Einstein equation (relating the viscosity and diffusion coefficients) is evaluated and commented.

II Method

II.1 Computational details

II.1.1 Ab initio molecular dynamics (AIMD)

Two AIMD simulations were run for this study, one for molten magnesite (\ceMgCO3) and one for molten dolomite (\ceCaMg(CO3)2). They were based on the density functional theory (DFT) within the Born-Oppenheimer approximation using the freely available QUICKSTEP/CP2K software VandeVondele et al. (2005a) that applies a hybrid Gaussian/plane-wave method Lippert, Hutter, and Parrinello (1997). For the study of \ceCaMg(CO3)2 we used for the valence electrons of carbon and oxygen a triple- zeta valence plus polarization basis set optimized for molecules (TZV2P)VandeVondele et al. (2005b). Otherwise we used a double-zeta plus polarization basis set (DZVP) VandeVondele and Hutter (2007). Core electrons were treated by the Goedecker-Teter-Hutter (GTH) norm conserving pseudo-potentials Goedecker, Teter, and Hutter (1996); Hartwigsen, Goedecker, and Hutter (1998); Krack (2005). The cutoff for the electronic density was set to 700 Ry. Exchange and correlation interactions were accounted for by the gradient corrected BLYP functional Becke A. (1988); Lee, Yang, and Parr R. (1988) using a semi-empirical D3 dispersion correction scheme with a cutoff $\geq\sqrt{3}L$, where $L$ is the length of the simulation box.Grimme et al. (2010) All DFT calculations were run in the $NVT$ ensemble with the temperature set constant by a Nosé-Hoover thermostat. Nosé (1984a, b)

For \ceMgCO3 and for \ceCaMg(CO3)2 the simulations were run with 200 atoms for 20 ps and with 640 atoms for 12 ps at state points (1873 K, 2.49 g/cm3) and (1773 K, 2.25 g/cm3), leading to an average pressure of 4.5 $\pm$ 1.5 and 1.8 $\pm$ 1.0 GPa, respectively. For information, congruent melting is known to occur at 2.7 GPa and 1850 K for magnesiteIrving and Wyllie (1975) and at 2.7 GPa and 1653 K for dolomite.Irving and Wyllie (1975) As for molten calcite (\ceCaCO3), the AIMD simulations of Vuilleumier et al. (2014), based on the same DFT approach, were used as benchmark.

II.2 Classical molecular dynamics (MD)

Classical MD simulations were carried out using the DL_POLY 2 software Smith and Forester (1996), with a timestep of 1 fs. Density calculations were performed in the $NPT$ ensemble with a Nosé-Hoover thermostat for a simulation time of $t\simeq$ 0.9 ns (including a 0.5 ns equilibration run) allowing to reach an accuracy on the density value of $\Delta n/n\sim\pm 1\%$. To evaluate the transport coefficients, simulations were performed in the $NVE$ ensemble with an equilibration run of 0.5 ns, followed by a production run of 10 to 30 ns. Structural data were extracted from the same simulations. All simulations had a system size $N\simeq$ 2000 atoms, except for specific calculations that used the number of atoms as the AIMD simulations in order to compare at best the pair distribution functions obtained by both approaches.

II.3 Force Field

The force field is composed of interionic pair potentials as presented in a previous paperDesmaele et al. (2019). The intramolecular part (interactions within a carbonate molecule) contains an oxygen-oxygen term: $V_{\text{OO}}^{intra-{\rm{CO}_{3}}}({r}_{\text{OO}})={A_{\text{OO}}^{intra-{\rm{CO}_{3}}}}\exp(-{r}_{\text{OO}}/{\rho_{\text{OO}}^{intra-{\rm{CO}_{3}}}})$, and a carbon-oxygen term: $V_{\text{CO}}^{intra-{\rm{CO}_{3}}}(r_{\text{CO}})={k_{\text{CO}}}(r_{\text{CO}}-{r_{0,\text{CO}}})^{2}/2+{q_{\text{O}}q_{\text{C}}}/{4\pi\epsilon_{0}r_{\text{CO}}}$, with a force constant ${k_{\text{CO}}}=6118.17$ kJ/mol and an harmonic equilibrium distance $r_{0,\text{CO}}=$1.30 Å  adjusted so as to obtain a mean C–O distance of 1.29 Å, as found by X-ray diffraction measurements and AIMD simulations Antao and Hassan (2010); Vuilleumier et al. (2014); Desmaele et al. (2019). Two ions $i$ and $j$, with $i,j$= Li, Na, K, Ca, Mg, O and C (with O and C not belonging to a same carbonate group) interact through the following intermolecular potential: $V_{ij}^{inter}({r}_{ij})={A_{ij}}\exp(-{r}_{ij}/{\rho_{ij}})-{C_{ij}}/{{r}_{ij}^{6}}+{q_{i}q_{j}}/{4\pi\epsilon_{0}{r}_{ij}}$. Table 1 recaps the previously published parameters for the Li–Na–K meltsDesmaele et al. (2019), and provides the parameters for Mg and Ca. The latter parameters were adjusted so as to reproduce at best (i) the density measurements of crystalline \ceMgCO3, \ceCaCO3 and \ceCaMg(CO3)2 at 300 K and 1 bar, the density and the compressibility of the \ceK2CO3–\ceCaCO3 melt at 1 bar and 1100–1200 K, and (ii) the microscopic structure, in the form of atomic pair distribution functions (PDFs), issued from AIMD simulations of the \ceCaCO3, \ceMgCO3 and \ceCaMg(CO3)2 melts.

At variance with alkali carbonates, alkaline-earth carbonates do not form stable melts under atmospheric pressure Irving and Wyllie (1975); Suito et al. (2001); Buob et al. (2006); Spivak et al. (2012); Shatskiy, Litasov, and Palyanov (2015) and no reliable measurement of their density under high pressure has been published yet. However molten calcium carbonate is stable at 1 bar in mixtures with alkalis. In particular, the density of the equimolar \ceK2CO3–\ceCaCO3 mixture has been measured by Liu and Lange (2003) using an advanced double-bob Archimedean method overcoming artifacts due to surface tension and to the mass of the bob, and which reproduced with a great accuracy the density of many liquid standards.Janz (1988) As shown on Fig. 1 the agreement between the density measurements of Liu and Lange (2003) on the K–Ca melt and our MD calculation is excellent. Note that Dobson *et al.*Dobson et al. (1996) also measured the density of this mixture but the values they report are by $\sim$ 5% lower than the ones from Liu and Lange.Liu and Lange (2003) Additionally, Liu and Lange (2003) have suggested that the density of molten carbonates virtually has an ideal behavior regarding composition Liu and Lange (2003); Desmaele et al. (2019):

[TABLE]

where $x_{i}$ is the molar fraction of species $i$ of molar mass ${M}_{i}$ and molar volume $\bar{V}_{i}$. Although this approximation is fairly good for purely alkali mixtures, it seems less accurate for mixtures containing both alkali and alkaline-earth cations (see Fig. 1 and next section). Still this mixing rule allows to extrapolate the density of an hypothetical \ceCaCO3 liquid at 1 bar.Liu and Lange (2003) In a MD simulation such a liquid is metastable, meaning that it is possible to evaluate its density straightforwardly, at variance with the experiments. Satisfactorily the density calculated by MD (e.g. 2.44 g/cm3 at 1100 K) is fairly compatible with the extrapolation proposed by Liu and Lange (2003) (2.49 g/cm3 at 1100 K).

For \ceMgCO3 no measurement of density, not even in a mixture, has been published yet. But based on the density of various carbonates (in the \ceLi2CO3–\ceNa2CO3–\ceK2CO3–\ceRb2CO3–\ceCs2CO3-\ceCaCO3–\ceSrCO3–\ceBaCO3 system) and considering how this property evolve as a function of the cationic radius, Hurt and Lange (2019) suggested that an hypothetical \ceMgCO3 melt at 1 bar and 1100 K would have a density of 2.45 g/cm3. This guess is within 2% of the density we obtained by simulating the metastable liquid (2.49 g/cm3). So by introducing in Eq. (1) the density of molten calcite (from Liu and Lange (2003)) and that of magnesite (from Hurt and Lange (2019)) at the chosen reference state of 1 bar and 1100 K gives a density of 2.47 g/cm3 for the \ceCaMg(CO3)2 melt, in close match with the value of 2.49 g/cm3 as obtained by MD (Table 2).

As the densities of molten magnesite and calcite proposed by Hurt and Lange (2019) are merely an approximate estimation, we chose to test the accuracy of our FF by calculating the density of the crystalline phases that are well constrained even under high $T-P$. Markgraf and Reeder (1985); Fiquet, Guyot, and Itie (1994); Zhang et al. (1997); Ross (1997); Redfern, Wood, and Henderson (1993) We chose to focus on the calcite structure (rhombohedral) as it is a common polymorph to \ceMgCO3, \ceCaCO3 and \ceCaMg(CO3)2.Hazen et al. (2013). For magnesite at 300 K up to 4 GPa (Fig. 2), the calculated density is in excellent agreement with the experimental values, well within the scattered data of the experimental literature. Moreover MD simulations reproduce very well the compressibility: $K_{T}$ = 125 GPa as compared to 117 GPa according to Ross.Ross (1997) We believe this should lead to a reliable evaluation of densities for the liquid phase. We also calculated the density of calcite (\ceCaCO3) at room conditions and obtained 2.63 g/cm3, against 2.71 and 2.72 g/cm3 according to Redfern *et al.*Redfern, Wood, and Henderson (1993) and Fiquet et al.Fiquet, Guyot, and Itie (1994), respectively, that is within $\sim$3%. For comparison, at the same conditions AIMD yielded 2.67 g/cm3.Vuilleumier et al. (2014) For dolomite, the agreement between MD ($n=2.81$ g/cm3) and the X-ray diffraction data ($n=2.84$ g/cm3)Fiquet, Guyot, and Itie (1994) is even better with a deviation of $-1$%. The densities were also calculated by anisotropic relaxation of the crystal structure (MD simulations in the $NST$ ensemble). The accuracy of the MD-calculated values deteriorates slightly for magnesite ($n=3.04$ g/cm3 instead of 3.01, see Fig. 2) but improves for dolomite ($n=2.84$ g/cm3) and remains unchanged for calcite. Tables S1 and S2 in the supplementary material detail the lattice parameters and the calculated density values up to 4 GPa.

III Thermodynamics

III.1 Equation of State

Although there exist no measured data of the density of calco-magnesian carbonate melts, an estimation can be made by assuming that carbonate liquids mix linearly with respect to carbonate components (ideal mixing assumption). As discussed above, this allowed Liu and Lange (2003) to estimate the density of molten \ceCaCO3 at 1-bar. Now, to infer the high pressure properties of these liquids, constraining the 1-bar compressibility is also needed. Then the density of molten carbonates can be accurately modeled by the third-order Birch Murnaghan equation of state (BMEoS) Birch (1947):

[TABLE]

where $n_{T}^{0}$ is the atmospheric density at temperature $T$, $K_{T}^{0}$ the bulk modulus (inverse of the compressibility) and $K^{\prime 0}$ its pressure derivative at 1 bar (note that it is a constant under the thermodynamic conditions of this study). After fitting our MD data, we propose an equation of state for \ceMgCO3, \ceCaMg(CO3)2 and for a natrocarbonatite (Na1.10K0.18Ca0.36CO3) modeling the carbonatitic lavas from Ol Doinyo Lengai (see Table 3). O’Leary *et al.*O’Leary, Lange, and Ai (2015) measured the compressibility of the alkali end-members and several mixtures of the system \ceCaCO3–\ceNa2CO3–\ceK2CO3–\ceLi2CO3, including mixtures containing various ratios of \ceCaCO3. Then they extrapolated the 1 bar compressibility of molten \ceCaCO3, that is a metastable liquid at this pressure, by using an ideal mixing rule:

[TABLE]

The comparison of our MD results with their study is good as shown in Table 2. In particular our results are in a better agreement than when using the FFs developed by Vuilleumier et al. (2014) and by Hurt and WolfHurt and Wolf (2018). Then O’Leary *et al.*O’Leary, Lange, and Ai (2009) used fusion curve analysis to constrain $K^{\prime 0}=7\pm 1$, still for \ceCaCO3. Note that this value is closer to the one produced using our FF ($K^{\prime 0}=8.2$, as in Vuilleumier et al. (2014)) than using the one of Hurt and Wolf (2018) ($K^{\prime 0}=10.3$). Inserting this value, as well as the 1-bar density from Liu and Lange (2003) and the 1-bar compressibility from O’Leary *et al.*O’Leary, Lange, and Ai (2009) in Eq. (2) we built an experimental compression curve for \ceCaCO3 at 1500 K (see the inset in Fig. 3b). This estimation is very close to the one issued from our MD simulations.

Based on AIMD simulations, Zhang and Liu (2015) proposed an equation of state for \ceCaCO3. However, because of the LDA approximation they used to compute exchange-correlation energies, an approximation which tends to overestimate the density (at a given $P)$, the authors have applied a rescaling method. The compression curve given by their equation of state at 1500 K is plotted in the inset of Fig. 3b for \ceCaCO3. It fits our MD compression curve quite well although there remains a slight deviation at low pressure. Moreover, the compressibility provided by our model is slightly lower than the one from Zhang and Liu (2015). However that may be, at 1500 K, we get $K_{1500K}^{0}=13.7$ GPa, just like O’Leary *et al.*O’Leary, Lange, and Ai (2015), whereas Zhang and Liu (2015) give a value of 11.9 GPa. As for the AIMD calculations we performed (this study and Vuilleumier et al. (2014)), they are based on the GGA approximation and include a (semi-empirical) correction for dispersion forces. The introduction of dispersion interactions enhances interionic cohesion in the liquid, but maybe not sufficiently, a feature which could account for the slight remaining difference between MD and AIMD ($P_{MD}<P_{AIMD}$ in general, see Fig. 3).

If the number density, $\rho$, is considered (rather than mass density $n$), it goes as follow at 3 GPa and 1873 K: $\rho_{\ce{MgCO3}}=$ 40.1 mol/L $>$ $\rho_{\ce{CaCO3}}=$ 25.2 mol/L, which is consistent with the observation made for alkali meltsDesmaele et al. (2019), that the number density is negatively correlated to the cation size. However, the bulk modulus is the same in magnesite and calcite ($K_{\ce{MgCO}_{3}}^{0}=K_{\ce{CaCO}_{3}}^{0}=11.0$ GPa), at variance with alkali melts where $K^{0}$ increases with increasing cation radius. This suggests that the compressibility of calco-magnesian carbonate melts is not similar to that of a hard sphere system and instead is dominated by the coulombic repulsion between divalent cations. If we consider the compressibility at 1 bar and 1100 K of the metastable \ceMgCO3 and \ceCaCO3 melts ($K_{\ce{MgCO}_{3}}^{0}=14.0$ GPa and $K_{\ce{CaCO}}^{0}=17.6$ GPa), it is greater by a factor $\sim$ 2 than that of alkali carbonates at the same $T-P$ conditions. However, the compressibilities calculated in corresponding states (i.e. near the melting point, 1823 K and 3 GPa for \ceMgCO3, and 1623 K and 1 GPa for \ceCaCO3) are very similar ($K_{1823\,\ce{K}}^{3\,\ce{GPa}}=7.5$ GPa for \ceMgCO3 and $K_{1623\,\text{K}}^{1\,\text{GPa}}=6.4$ GPa for \ceCaCO3) for both alkali and alkaline-earth carbonates.

Knowing the equations of state for \ceMgCO3, \ceCaCO3, \ceNa2CO3, \ceK2CO3 and \ceLi2CO3 and some mixtures (dolomite and natrocarbonatite) it is worth checking how these latter melts behave regarding ideality (Table 2). For that we consider, for instance, the density of molten dolomite at 1 bar and at 3 GPa, let’s say at 1873 K. The comparison between the raw MD data and the ideal mixing rules (Eqs. (1) and (3)) using the data of Table 3 gives: $n_{T}^{0}=n_{T}^{0,\,mix}=2.19$ g/cm3 and $K_{T}^{0}=K_{T}^{0,\,mix}=11.0$ GPa, and $n_{T}^{3\,\text{GPa}}=n_{T}^{3\,\text{GPa},\,mix}=2.55$ g/cm3. Proceeding the same way for the natrocarbonatite at 1 bar and 1073 K we get $n_{T}^{0}=2.11$ g/cm3 and $n_{T}^{0,mix}=2.13$ g/cm3, and $K_{T}^{0}=11.4$ GPa and $K_{T}^{0,\,mix}=10.4$ GPa. Hence it comes out that dolomite is an ideal mixture on a certain $P$-domain, whereas natrocarbonatite slightly deviates from ideality. Most likely the latter finding can be explained by the fact that the natrocarbonatite includes both uni- and divalent cations. In practice, as the deviation from ideality seems to be fairly small, the EoS of the end-members of the Mg–Ca–Li–Na–K system can be used to estimate with confidence the density of any mixture.

III.2 Structure

In the $T-P$ range of this study ($T\leq 2013$ K and $P\leq 15$ GPa), the stability of the internal structure of the carbonate ion has been pointed out by neutron diffraction measurementsKohara et al. (1998), in situ X-ray diffraction measurementsHudspeth, Sanloup, and Kono (2018) and by AIMD simulations (Ref. Vuilleumier et al., 2014; Desmaele et al., 2019 and this study). Moreover the carbonate anion, unlike the \ceSiO4 units in silicates, cannot share a covalent bond with other atoms. Therefore carbonate melts display a structure contrasting with that of most geological melts (silicate) that are always polymerized to a certain degree. This is why carbonate melts share with molten salts (e.g. NaCl) a high ionic diffusivity and a low viscosity.

Until now, among the alkaline-earth carbonates, only the structure of calcite melt has been investigated by classical molecular dynamics simulationsGenge, Price, and Jones (1995); Vuilleumier et al. (2014); Hurt and Wolf (2018) and by AIMD simulations.Vuilleumier et al. (2014); Du et al. (2018) Fig. 4 shows the PDFs for \ceMgCO3 and \ceCaCO3 at corresponding states (i.e. near melting point: $T=1873$ K, $P\sim 2$ GPa and $n=2.49$ g/cm3 for \ceMgCO3; $T=1623$ K, $P\sim 0.5$ GPa and $n=2.30$ g/cm3 for \ceCaCO3). Note that, as expected, they are in very good agreement with the PDFs issued from AIMD (compare the plain lines with the dotted lines). The first peak of $g_{\text{CO}}(r)$ at $r_{\text{CO}}=$ 1.29 Å  corresponds to the three oxygen atoms bonded to a same carbon atom in a carbonate unit. Between 1.7 Å  and 2.1 Å  $g_{\text{CO}}(r)$ = 0 because no C–O bond dissociation occurs during the simulation. On the other hand, the first peak of $g_{\text{OO}}$ represents the O–O intramolecular distance ($r_{\ce{OO}}=2.22$ Å) in a carbonate ion. The first minimum of $g_{\ce{OO}}$ is non-zero, meaning that the distance between two oxygen atoms of two different \ceCO32- units may be as short as the O–O intramolecular distance. Concerning the anion-anion correlations, each \ceCO32- unit is surrounded by $12-15$ other carbonate ions at $r_{\ce{CC}}\simeq 4-6$ Å  on average. The first C–C peak on Fig. 4a is broad with a shouldering on its low-$r$ flank ($3.5-4$ Å), that is especially noticeable for \ceMgCO3, and a small bump at $\sim~{}5-7$ Å  in the region of the first minimum, this time more pronounced in \ceCaCO3. Each cation (Mg2+ or Ca2+) ion is surrounded by 6 carbonate groups (against $4.5-4.8$ for alkali), at a mean cation-carbon distance of $3-4$ Å  in \ceMgCO3 and $3.15-4.5$ Å  in \ceCaCO3 (note that these distances are the same in \ceCaMg(CO3)2 at 1773 K and 2.25 g/cm3). As for the ratio of the coordination numbers $N_{c}^{\text{X-O}}$/$N_{c}^{\text{X-C}}$ (where X = Mg or Ca), indicative of the orientation of the carbonate ions around the cations, it is 1.1 for Mg and 1.3 for Ca (as compared to $1.3-1.6$ for alkali cations). In comparison to pair distribution functions for alkali carbonate melts, the cation-anion PDF $g_{\ce{XC}}(r)$ for alkaline-earth carbonates are more simple as they show no shoulder on the first peak (see Fig.4b and compare with Fig. 2 in Desmaele et al. (2019)). Moreover, the second peak of $g_{\ce{XC}}(r)$ is broader, especially for Mg. It can also be noted that the shape of the PDFs is less sensitive to the size of the cation. Looking at the cation-cation PDFs, the only difference between Mg and Ca is the amplitude, greater for the Ca–Ca pair with a coordination number around $\sim 10$ instead of $\sim 7$ for Mg–Mg. It is noteworthy that the PDFs of \ceCaMg(CO3)2 are intermediate between those of \ceCaCO3 and \ceMgCO3. In fact, most structural features observed in molten calcite and magnesite are similar and can be interpreted as a simple homothetic transformation upon volumetric change from \ceCaCO3 to \ceMgCO3. Thus the pair distribution functions X–C, X–O and X–X in \ceCaMg(CO3)2 are almost identical to the corresponding ones observed in \ceCaCO3 and \ceMgCO3 (Fig. 4b) and this remains true even at high pressures (Figs. S1, S2 and S3 in the supplementary material). Under pressure (and up to 12 GPa), the PDFs shift progressively towards lower distances, reflecting the melt compaction (see Figs. S1, S2 and S3). Moreover the average number of \ceCO32- anions around cations increases slightly under pressure, from 6 to 7. In contrast, the coordination number of \ceCO32- around \ceCO32- ($\sim 12-15$) does not evolve with pressure, neither does the $N_{c}^{\text{X-O}}$/$N_{c}^{\text{X-C}}$ ratio.

IV Transport Properties

For \ceMgCO3, \ceCaCO3 and \ceCaMg(CO3)2, a series of simulations ($\sim$ 15) was performed at different thermodynamic conditions, with a duration long enough to reach the diffusive regime ($10-20$ ns). From each run we calculated accurately (see below) the self-diffusion coefficients $D_{s}$ of each chemical species $s=$ Ca, Mg and CO${}_{3}^{2-}$, the electrical conductivity $\sigma$ and the viscosity $\eta$, given by Allen and Tildesley (1989); Hess (2002):

[TABLE]

where $\Pi_{\alpha\beta}(t)$ refers to the off-diagonal pressure tensor components ($\alpha,\beta=x,y,z$, see references Allen and Tildesley (1989); Hess (2002) for more details), $N$ is the total number of ions in the simulation box of volume $V$, $N_{s}$ the number of ions of species $s$, $k_{{\rm B}}$ the Boltzmann constant and $e$ the elementary charge, $m_{i}$ is the mass of ion $i$, $\mathbf{r}_{i}(t)$ its position, $v_{i\alpha}(t)$ the component $\alpha$ of its velocity, $F_{ij\beta}(t)$ is the component $\beta$ of the force exerted by ion $j$ on ion $i$, separated by a distance $r_{ij}(t)$ at time $t$. The conduction charge $z_{i}$ is taken as the formal charge, which is usual for simple ionic liquids. Adams, McDonald, and Singer (1977) Note from the equations above that the self-diffusion coefficient $D_{s}$, resulting from an average over the $N_{s}$ ions of a specific species $s$ (sum $\sum_{i=1}^{N_{s}}$ in Eq. (4)), are calculated with a great accuracy (calculation uncertainty of 1% in this study). On the other hand, the electrical conductivity and the viscosity, as collective observables, are trickier to estimate due to the slow convergence of the corresponding time correlation functions. In the following the values we present for these quantities have an error bar within $5-10$%.

The temperature and pressure dependence of the transport coefficients can be modeled by an Arrhenius activation law

[TABLE]

where $E_{a}^{X}$ is the activation energy associated to the physical quantity $X$ and $V_{a}^{X}$ is an activation volume accounting for its pressure dependence.Vuilleumier et al. (2014); Corradini, F.-X., and Vuilleumier (2016); Desmaele et al. (2019). The values of $X^{0}$, $E_{a}^{X}$ and $V_{a}^{X}$ were determined by fitting the molecular dynamics data for the three melts \ceMgCO3, \ceCaCO3 and \ceCaMg(CO3)2 (see Figs. 5, 6, 7, 10, 11). They are given in Table 4 for $D_{s}$ and in Table 5 for $\sigma$ and $\eta$, for the pressure and temperature range mentioned in Table 3.

IV.1 Diffusion coefficient

The diffusivity of calco-magnesian carbonate melts has never been measured and it is only recently that estimations have been provided by the AIMD simulations of Vuilleumier et al. (2014). According to their calculations the diffusion coefficients of Ca and CO${}_{3}^{2-}$ in \ceCaCO3 along its melting curve (up to 12 GPa) are of comparable magnitude with the ones in purely alkali melts.Desmaele et al. (2019) The agreement between the data of the present MD study and the AIMD simulations of Vuilleumier et al. (2014) is very good at 0.5 and 4.5 GPa, a bit less at 12 GPa where the values of MD are below those of AIMD (Fig. 5). We have also evaluated the diffusion coefficients from our AIMD simulations of \ceMgCO3 at 4.5 GPa and \ceCaMg(CO3)2 at 2 GPa. In the two cases they are greater than the values issued from MD (Figs. 6 and 7). This is consistent with the slight discrepancy between the two models for the equation of state. Indeed at a given $T,P$ point, the densities yielded by MD simulations are systematically greater than the ones yielded by AIMD (Fig. 3). This means that the free volume of diffusion is smaller in the MD model. By performing a MD simulation of \ceMgCO3 at the same density as AIMD (1873 K and 2.49 g/cm3), and recalculating the coefficients from this run, we get larger values (namely $D\textsubscript{Mg}=3.11$ and D\textsubscript{CO{}_{3}}=2.13~{}10^{-9} m2/s, instead of 1.85 and 1.30), but still below the values from AIMD: $4.4\pm 0.5$ and $3.5\pm 0.5~{}10^{-9}$ m2/s, respectively. This could indicate that our empirical force field fails to some extent to describe cohesive forces in every detail. However there is no evidence that AIMD sketches them much more accurately.

The activation energy ($51-58$ kJ/mol) depends little on the ion species (Ca, Mg or \ceCO3). But it is perceptibly higher than that for alkali melts at 1 bar in which the coulombic forces are weaker.Desmaele et al. (2019) The magnitude of the diffusion coefficients are also slightly smaller in Ca-Mg carbonate melts than in their alkali counterparts. Interestingly if we consider the \ceK2CO3-\ceCaCO3 mixture (at 1 bar and $1100-1200$ K) $D$Ca is lower than $D$K (by a factor $\sim$2) and very close to $D$CO3. Moreover, compared to pure \ceK2CO3, the presence of the divalent cation decreases $D$CO3 and $D$K by a factor 2 and increases their activation energies by $\sim+50\%$.

IV.2 Electrical Conductivity

As molten salts, carbonate melts are characterized by a high electrical conductivity in the range of $10^{1}-10^{3}$ S/m, depending on composition, temperature and pressure, which is up to two orders of magnitude more conductive than silicate melts at the same thermodynamic conditions.Gaillard et al. (2008); Sifré et al. (2014); Sifré, Hashim, and Gaillard (2015) The knowledge of the electrical conductivity of alkali carbonates is crucial for their industrial applications as electrolytes in fuel cell devices. Many experimental and numerical studies have been devoted to this issue. Hence, the electrical conductivity of the end-members and of binary and ternary mixtures of the system \ceLi2CO3–\ceNa2CO3–\ceK2CO3 have been abundantly documented.Janz (1988); Kojima et al. (2007, 2008); Desmaele et al. (2019) Because the addition of small amounts of alkaline-earth carbonates improves the performance of fuel cell devices, in particular in terms of durability,Kojima et al. (2003) the electrical conductivity of the Li–K (62–38 mol%) molten carbonate was measured by impedance spectroscopy and found to decrease linearly with small amounts of \ceCaCO3 Lair et al. (2012). The electrical conductivity of carbonate melts is not only important for industrial applications, it is also of fundamental interest to understand the conductivity anomalies in the asthenosphere of the Earth’s mantle. To address this question, Gaillard et al. (2008) measured the electrical conductivity of binary and ternary mixtures in the Na–K–Ca system at atmospheric pressure. With regard to Ca-bearing mixtures, MD is in quantitive agreement with their values (deviation of 30% at most), but with a slightly higher activation energy (Fig. 8). Following the work of Gaillard et al. (2008), Sifré *et al.*Sifré, Hashim, and Gaillard (2015) studied (up to 3 GPa) Ca and Mg-bearing carbonate compositions: \ceCaCO3, a natural dolomite ($\sim$ \ceCaMg(CO3)2), \ceK2Mg(CO3)2, \ceK2Ca(CO3)2 and \ceNa2Mg(CO3)2. These studiesGaillard et al. (2008); Sifré, Hashim, and Gaillard (2015) show that the electrical conductivity depends slightly on the chemical composition (see Figs. 8 and 9). The smaller the cation and the lower its charge, the higher is the electrical conductivity.Kojima et al. (2008); Sifré, Hashim, and Gaillard (2015). As a consequence, the addition of \ceCaCO3 or \ceMgCO3 in a pure alkali carbonate reduces somewhat the conductivity.

Fig. 9 reports the calculated and measured conductivities for \ceCaCO3 and for the above mentioned Ca or Mg-bearing mixtures at 3 GPa. The agreement is very good for dolomite and for the \ceNa2Mg(CO3)2 mixture. Note that for dolomite at 3 GPa and 1800 K, Yoshino et al. (2012) reported an electrical conductivity of 105 S/m, which is almost twice lower than the values of this study and from Sifré et al..Sifré et al. (2014) Concerning \ceCaCO3, the calculated and the measured electrical conductivities overlap within uncertainty (e.g. at 3 GPa and 1923 K MD yields $210\pm 10$ S/m and Sifré *et al.*Sifré, Hashim, and Gaillard (2015) $230\pm 25$ S/m). The agreement between MD and experiments is also good for \ceK2Ca(CO3)2, although it slightly degrades towards low temperatures. Most striking is our disagreement on the \ceK2Mg(CO3)2 mixtures (by a factor $\sim 2-3$, Fig. 9), although our force field satisfactorily reproduced the behavior of other Mg-containing mixtures. An explanation could be that this composition is the most dissymetric one of our study (a small divalent cation coexists with a large monovalent cation in equal proportions). On the other hand, an experimental bias cannot be excluded because this mixture is known for being glass-forming and for easily decomposing Dobson et al. (1996).

As already emphasized in previous studies, the increase of conductivity with temperature is well fitted by an Arrhenius law with activation energies ranging from 34 to 42 kJ/mol at $P=1$ atm (Table 5). Note that Gaillard et al. (2008) found $E_{a}^{\sigma}$ in the range 30 to 35 kJ/mol in alkali-bearing melts and Desmaele et al. (2019) $\sim$ 20 kJ/mol for purely alkali melts by MD. The conductivity decreases weakly with pressure, which can be accounted for by an increase of the activation energy with $P$. Hence, we calculated that this energy is between 42 and 51 kJ/mol at 3 GPa depending on the melt composition, whereas Sifré *et al.*Sifré, Hashim, and Gaillard (2015) reported values from 37 to 48 kJ/mol. For dolomite, these authors found an activation energy of 48 kJ/mol, while our calculated value is 41 kJ/mol. That is a fairly good agreement considering that the experimental error on each measurement is about 10%. As for Yoshino et al. (2012) they reported a value of 38 kJ/mol at 3 GPa. With regard to the calcite melt (\ceCaCO3), the activation energy calculated by MD (41 kJ/mol at 3 GPa) is much lower than the value reported by Sifré *et al.*Sifré, Hashim, and Gaillard (2015) (76 kJ/mol at the same pressure), but consistent with the ones they reported for other carbonate compositions (in the range $37-48$ kJ/mol).

The electrical conductivity estimated from AIMD simulations (this study for \ceCaMg(CO3)2 and Vuilleumier et al. (2014) for \ceCaCO3) are reported on Figs. 10. For \ceCaCO3 the agreement between MD and AIMD is good at 0.5 GPa but diminishes at higher pressures. As evoked in the case of diffusion coefficients, most of these discrepancies may be related to the small difference between the equation of state provided by the two models. At a given $T-P$, the density is smaller in the AIMD simulations, hence the ionic mobility and the electrical conductivity are higher. For dolomite, given the large uncertainty on the AIMD value ($\sim 50$%), it can only be stated that the results of the two simulations models (MD and AIMD) are compatible.

Another advantage of the MD approach is that phenomenological relations for transport properties can be tested. It can be shown that the electrical conductivity $\sigma$ and the diffusion coefficients are linked by a generalized Nernst-Einstein relation, which takes into account the cross correlations between ionic motions, $\sigma=H\sigma^{NE}$ where

[TABLE]

with $N_{s}$ , $z_{s}$ and $D_{s}$ are the number, the formal charge and the self diffusion coefficient of ions of atomic species $s$. The Nernst-Einstein relation assumes that ions move independently from each others. The Haven ratio:

[TABLE]

accounts for the average cross correlations (through a scalar product) between the displacements of ions of species $s$

[TABLE]

and the average cross correlations between the displacement of an ion $i$ of species $s$ ($\vec{\Delta}_{i}^{(s)}$) and that of an ion $j$ of another species $s^{\prime}$ ($\vec{\Delta}_{j}^{(s^{\prime})}$)

[TABLE]

The Nernst-Einstein equation is recovered for $H=1$, although it doesn’t imply that $H_{s}$ and $H_{ss^{\prime}}\simeq 0$, as these two terms can cancel each other (see Desmaele et al. (2019)). For instance in \ceMgCO3 (for further details see Table 6), $H$ is close to 1 ($H=0.92$) although its decomposition gives $H_{\rm Mg}=-0.14$, $H_{\rm CO_{3}}=-0.32$ and $H_{\rm Mg-CO_{3}}=0.38$. Furthermore the fact that both $H_{\rm Mg}$ and $H_{\rm CO_{3}}$ contributions are negative implies that ions of a same species have a high probability of moving towards opposite directions (see Eq. (12)), which decreases the conductivity. Regarding the cation-anion correlation term, $H_{\rm Mg-CO_{3}}=0.38$, it is positive and almost exactly cancels out the $H_{\rm Mg}+H_{\rm CO_{3}}$ sum. Because Mg and CO3 ions have charges of opposite signs (see Eq. (13)), $H_{\rm Mg-CO_{3}}>0$ is indicative of an anti-correlation of the displacements of the cations and of the anions, that is to say, they are, on average, moving towards opposite directions, which increases the conductivity. For \ceCaCO3, an anti-correlation of all the ionic displacements is also observed: a positive contribution of the anion-cation correlation $H_{\rm Ca-CO_{3}}$ is not fully overbalanced by the correlations between ions of a same species leading to a Haven ratio less than 1 ($H=1+H_{\rm Ca}+H_{\rm CO_{3}}+H_{\rm Ca-CO_{3}}=0.81$).

We carried on this approach with some binary mixtures (see Table 6), including dolomite. It is worth noticing that the displacements of cations of the same species are mostly not correlated to one another ($H_{\rm X_{1}}$ and $H_{\rm X_{2}}$ both $\sim 0$, where X1, X2 = Mg, Ca, Na or K). However the cross correlations between cationic displacements (as expressed by the term $H_{\rm X_{1}-X_{2}}$) seem to depend on the charge of the cations. Indeed, in dolomite $H_{\rm X_{1}-X_{2}}$ is close to 0 ($H_{\rm Mg-Ca}=-0.03$), while in the investigated alkali-alkaline earth mixtures it ranges between -0.16 and -0.22. As for the other terms ($H_{\rm CO_{3}}$ and $H_{\rm X-CO_{3}}$), they give opposite contributions to $H$ (see Table 6). So, although the Nernst-Einstein equation yields, for the melts studied here, a reasonable estimation of the electrical conductivity ($+10-30$% from the exactly calculated one, $\sigma$), it provides no information on the relevance of its underlying assumption (the ions move independently from each other). In fact, as we show it here, its usefulness relies on a cancellation effect between the different ion-ion correlations. As the dependence of this cancellation effect with composition, $T$ or $P$ is non-trivial, we recommend a circumspect use of the Nernst-Einstein approximation.

IV.3 Viscosity

Very little is known on the pressure dependence of the viscosity of molten carbonates. Based on molten salt data reviewed by Janz,Janz (1988) Wolff (1994) roughly estimated an order of magnitude ($\sim$ 100 mPa$\cdot$s) for the viscosity of calcium-rich carbonatites at 973 K and ambient pressure. The ex-situ measurements at mantle pressures ($P=3$ GPa) made by Sykes et al.Sykes, Baker, and Wyllie (1992) seem to greatly overestimate the viscosity of dolomitic melts (e.g. for 70:30 mol% \ceCaCO3–\ceMgCO3, $\eta=600$ mPa.s at 1473 K and 1 GPa) in light of the very high activation energy they obtained Kono et al. (2014), probably because of an incomplete melting of the sampleJones, Dobson, and Genge (1995). On the other hand, Dobson et al. (1996) performed the first in-situ measurements of the viscosity of \ceK2Ca(CO3)2 and \ceK2Mg(CO3)2 melts, up to 5.5 GPa, by using the falling sphere method with X-ray radiography. In this pioneering experiment, the temperature was difficult to control and the relative error on viscosity was typically 50%, due to incomplete melting of the sample and to convection effects at the highest temperatures. Besides the limited frame rate used to capture the images of the falling sphere greatly reduces the accuracy of the final velocity measurement, which is crucial for determining the viscosity Kono et al. (2014). Nevertheless it was found that Ca, Mg-bearing carbonates under high $T-P$ have a viscosity similar to that of alkali carbonates at 1 bar ($6-36$ mPa$\cdot$s in the range $2.5-5.5$ GPa). This qualitative observation is consistent with the results of our MD calculations (see Fig. 11). However we do not agree with the assertion of the authorsDobson et al. (1996) that the effect of pressure on viscosity is negligible in the pressure range investigated (up to 5.5 GPa) and we think that this observation results from the large error bars of the study ($\sim$ 50% on $\eta$ and $\pm$ 0.5 GPa on $P$, see also the discrepancy with the first results of Jones et al.Jones, Dobson, and Genge (1995)). We calculated the viscosity of \ceK2Ca(CO3)2 at 0, 2.5 and 3 GPa and the viscosity of \ceK2Mg(CO3)2 at 3 GPa. The viscosity depends strongly on the temperature by following an Arrhenius law with an activation energy a little higher for \ceK2Mg(CO3)2 than for \ceK2Ca(CO3)2. Even at these moderate pressures it is obvious that the activation energy depends on pressure (Fig. 12), in contrast with the observation of Dobson et al. (1996). At a given pressure, we obtained similar activation energies for the two compositions \ceK2Mg(CO3)2 and \ceK2Ca(CO3)2 ($\sim$ 80 kJ/mol at 3 GPa). However these activation energies differ from the ones calculated for the end-members \ceCaCO3 and \ceMgCO3 ($43-49$ kJ/mol), for dolomite (60 kJ/mol) and for alkali melts ($31-39$ kJ/mol)Desmaele et al. (2019) at the same pressure. This illustrates that the viscosity not only depends on the temperature, but also on the chemical composition. Recently, Kono et al. (2014) have reported viscosities of calcite and dolomitic (Mg0.40Fe0.09Ca0.51CO3) melts at temperatures just above the melting point up to 6.2 GPa, by using the falling sphere method with a powerful technique of ultra-fast synchrotron X-ray imaging. The measured viscosities ($6-10$ mPa.s, with an error of 9%) are of the same order of magnitude as those of Dobson et al. (1996), although the composition and the temperature differ in the two studies. The values of Kono et al. (2014) (at $0.9-6.2$ GPa and $1653-2063$ K for calcite, $3.0-5.3$ GPa and $1683-1783$ K for dolomite) are reported on Fig. 11. For calcite a satisfying agreement (within $10-15$%, i.e. within the overlap of error bars) is found with our calculations, better than with the values of Vuilleumier et al. (2014) For dolomite the agreement, although comprised between 3 and 30%, is reasonable once considered some anomalous trend in the experimental data. Indeed the viscosities measured by Kono et al. (2014) at 3.9 GPa and $1000/T\sim 0.6$ is lower than the one measured at 3.0 GPa at the same temperature, whereas the viscosity is expected to increase upon increasing pressure. Incidentally, we do not think that the small content (9 mol%) of \ceFeCO3 in the experimental composition has a significant contribution to the viscosity and could account for the MD-experiment discrepancy.

For natrocarbonatite (Natro, Na1.1K0.18Ca0.36CO3) at 1 bar and 823 K (near eruption conditions at the Ol Doinyo Lengai volcano), we calculated a viscosity of 68 $\pm$ 10 mPa$\cdot$s. We believe this is a more reliable value than the estimation given by Treiman and Schedl (1983) using molten salt data (5 mPa$\cdot$s), which seems rather low for a Ca-bearing mixture at relatively low temperature. For example compared with the alkali ternary eutectic mixture (Li0.435Na0.315K0.25CO3) the viscosity of natrocarbonatite is greater by a factor $\sim 3$ according to our MD calculations.Desmaele et al. (2019) This is consistent with the observations made by both experiments and simulations, of a decrease of the ionic conductivity upon addition of the alkaline-earth cation Ca.Gaillard et al. (2008); Lair et al. (2012) Interestingly the calculated viscosity is much lower than the ones measured by Norton and Pinkerton (1997) for several carbonatitic melts issued from the eruption of Ol Doinyo Lengai in 1988. At 823 K, the measured viscosities vary a lot ($10^{2}-10^{4}$ mPa/s) as a result of varying compositions (with or without a silicate component), crystallinity (presence of crystals of various sizes) and vesicularity (presence of \ceCO2 bubbles). Consequently, for a composition close to the one we studied (Na1.1K0.18Ca0.36CO3), they reported a viscosity as high as $\sim$ 600 mPa$\cdot$s. At variance with experiments where a complete melting is uneasy to assess, the MD value (68 $\pm$ 10 mPa$\cdot$s) is representative of carbonate melts. Besides, this oil-like viscosity is consistent with field observations where effusion of a very fluid lava is seen (see the videos of Fischer et al. (2009) with the link of reference NSF, 2009).

To obtain a phenomenological description of the viscosity of simple liquids, the Stokes-Einstein equation is often used

[TABLE]

where $d$ is the diameter of the diffusing particle and $D$ its diffusion coefficient. In the present case, the latter is assumed to be given by the arithmetic mean of the diffusion coefficients, $D=\sum_{s}x_{s}D_{S}/\sum_{s}x_{s}$ where $x_{s}$ is the molar fraction of ion of species $s$ in the melt. We chose $d=3.2$, 3.5 and 3.4 Å  for \ceMgCO3, \ceCaCO3 and \ceCaMg(CO3)2 respectively. The behavior of the MD-calculated viscosity with pressure and temperature was very well reproduced by Eq. (14) using these parameters (Fig. 13). To calculate the hydrodynamic diameter in a more grounded framework we used the following expression: $d=x_{\text{X}}^{2}d_{\text{XX}}+x_{\text{CO}_{3}}^{2}d_{\text{CO}_{3}\text{CO}_{3}}+2x_{\text{X}}x_{\text{CO}_{3}}d_{\text{X}\text{CO}_{3}}$, where $x_{\text{X}}$ and $x_{\text{CO}_{3}}$ are the molar fractions of cation X and of anion CO3, respectively, and $d_{\text{XX}}$, $d_{\text{CO}_{3}\text{CO}_{3}}$ and $d_{\text{X}\text{CO}_{3}}$ the cation-cation, carbonate-carbonate and cation-carbonate distances issued from the closest approach distances indicated by the corresponding PDFs (in fact, the $d$ parameters correspond to the distances at which the integral of the PDFs is equal to 1). This ansatz leads to $d=3.1$, 3.4 and 3.3 Å  for \ceMgCO3, \ceCaCO3 and \ceCaMg(CO3)2 respectively. These values are almost identical ($-0.1$ Å) to the ones obtained by fitting the values of $d$ to Eq. (14).

Finally, the electrical conductivity $\sigma$ (in S/m) and the viscosity $\eta$ (in Pa.s) can be related by the following empirical formula:

[TABLE]

where $A=$ 4.80, 3.25 and 3.58 for \ceMgCO3, \ceCaCO3 and \ceCaMg(CO3)2 respectively (Fig. 14). A comparable relation between these two transport coefficients has been also highlighted experimentally and computationally for various melt compositions. Sifré, Hashim, and Gaillard (2015); Vuilleumier et al. (2014); Desmaele et al. (2019).

V Conclusion

Following our previous work on the \ceLi2CO3–\ceNa2CO3–\ceK2CO3 melts, we have studied the thermodynamics, the microscopic structure and the transport properties (diffusion coefficients, electrical conductivity and viscosity) of Ca and Mg-bearing carbonate melts up to 2073 K and 15 GPa. For that we have developed an empirical force field benchmarked on data from experiments (density of the rhombohedral crystal phases at 300 K and up to 4 GPa, and the density and the compressibility of the \ceK2Ca(CO3)2 melt) and from AIMD simulations (microscopic structure of five liquids). The density and compressibility, evaluated for the metastable melts of \ceMgCO3 and \ceCaCO3 at 1100 K and 1 bar, are in very good agreement with the estimates of the experimental literature. Moreover we have shown that alkaline-earth carbonate mixtures behave ideally regarding the density and the compressibility. Based on the example of a Na–Ca–K melt (natrocarbonatite), the assumption of an ideal behavior for alkali-alkaline-earth mixtures seems a little less accurate but still reasonable.

The equations of state of carbonate melts with a composition of prior interest in the study of the Earth’s mantle (\ceMgCO3, \ceCaCO3 and \ceCaMg(CO3)2), as well as that modeling the natrocarbonatite emitted at the Ol Doinyo Lengai volcano (Na1.10K0.18Ca0.36CO3), were evaluated and modeled by a third-order Birch-Murnaghan formula. Covering a large $T$, $P$ range, these data may help in the debate on the geodynamics of carbonatitic melts relative to silicate melts.Liu, Tenner, and Lange (2007)

The analysis of the PDFs associated to the dolomitic melt showed an ideal behavior of the microscopic structure, as the PDFs in \ceCaMg(CO3)2 are similar to the corresponding ones in \ceCaCO3 and in \ceMgCO3. Moreover, no major modification of the structure was observed with increasing pressure.

As for the transport properties (diffusion coefficients, electrical conductivity and viscosity), they evolve smoothly (Arrhenius-like) over the studied $T-P$ domain. Thus at the high $T-P$ conditions of the Earth’s mantle, Mg and Ca-bearing molten carbonates keep the main features of molten salts, namely the ions are highly mobile. For example the viscosity is comprised in the range $1-100$ mPa$\cdot$s, that is in between the one of water and the one of olive oil at room conditions. However, the composition has a non negligible effect on transport coefficients: calco-magnesian melts are systematically more viscous than their alkali counterparts. These results may provide new insights into magmatic processes implying carbonatitic melts.Treiman and Schedl (1983); Treiman (1989); Wolff (1994).

Finally we discussed the reliability of the phenomenological Nernst-Einstein and the Stokes-Einstein equations, that relate the diffusion coefficients to the electrical conductivity and to the viscosity, respectively. These relations are often used to infer one of the transport coefficient from another one that has been measured. According to the present study, both formulas lead to reasonable values. However the underlying assumptions of the relations are not always representative of the transport mechanism itself. We have shown that the fairly good predictions provided by the Nernst-Einstein equation (which assumes that the ions move independently from one another) result from a partial cancellation of interionic dynamic correlations whose dependence with composition, $T$ or $P$ is non-trivial. As a consequence we recommend a circumspect use of the Nernst-Einstein approximation.

In summary, the overall agreement between the results of the MD simulations using this force field and the full set of available experimental data (thermodynamics and transport coefficients) provides evidence of the ability of our FF to describe with accuracy the properties of any melt in the \ceMgCO3–\ceCaCO3–\ceLi2CO3–\ceNa2CO3–\ceK2CO3 system.

Supplementary Material

See supplementary material for the density of crystal phases calculated by anisotropic relaxation of the rhombohedral structures up to 4 GPa, the evolution of the PDFs issued by MD with increasing pressure (up to 15 GPa), a summary of all calculated properties with their uncertainties.

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). The authors acknowledge GENCI for HPC resources (Grant No. 2015-082309).