Atomistic simulations of molten carbonates: thermodynamic and transport properties of the Li2CO3-Na2CO3-K2CO3 system

TL;DR

This study develops a comprehensive atomistic simulation model for molten alkali carbonates, accurately predicting their thermodynamic and transport properties across conditions relevant to Earth's mantle and technological applications.

Contribution

It introduces a novel, validated classical force field for molten alkali carbonates applicable over wide temperature, pressure, and compositional ranges.

Findings

Good agreement with experimental data on thermodynamics and transport properties.

First reliable molecular model covering extensive conditions for molten alkali carbonates.

Provides insights into geochemical processes and potential technological uses.

Abstract

Although molten carbonates only represent, at most, a very minor phase in the Earth's mantle, they are thought to be implied in anomalous high-conductivity zones in its upper part (70-350 km). Besides the high electrical conductivity of these molten salts is also exploitable in fuel cells. Here we report quantitative calculations of their properties, over a large range of thermodynamic conditions and chemical compositions, that are a requisite to develop technological devices and to provide a better understanding of a number of geochemical processes. To model molten carbonates by atomistic simulations, we have developed an optimized classical force field based on experimental data of the literature and on the liquid structure issued from ab initio molecular dynamics simulations performed by ourselves. In implementing this force field into a molecular dynamics simulation code, we have evaluated the thermodynamics (equation of state and surface tension), the microscopic liquid structure and the transport properties (diffusion coefficients, electrical conductivity and viscosity) of molten alkali carbonates (Li2CO3, Na2CO3, K2CO3 and some of their binary and ternary mixtures) from the melting point up to the thermodynamic conditions prevailing in the Earth's upper mantle (~ 1100-2100 K, 0-15 GPa). Our results are in very good agreement with the data available in the literature. To our knowledge a reliable molecular model for molten alkali carbonates covering such a large domain of thermodynamic conditions, chemical compositions and physicochemical properties has never been published yet.

Tables

Table 1. Table 1: Intermolecular potential parameters and partial charges. Note that the 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)
Li	O	3.0 $\times {10}^5$	0.2228	1210	+0.82101
Na	O	11 $\times {10}^5$	0.2228	3000	+0.82101
K	O	9.0 $\times {10}^5$	0.2570	7000	+0.82101
O	O	5.0 $\times {10}^5$	0.252525	2300	$-$0.89429
C	O	0	1	0	+1.04085

Table 2. Table 2: Density and thermal expansion coefficient ( n = n 1100 ​ K − α ​ ( T − 1100 ) 𝑛 subscript 𝑛 1100 𝐾 𝛼 𝑇 1100 n=n_{1100\,K}-\alpha(T-1100) ) from MD simulations and comparison with the experimental literature: references Janz, 1988 ( ◆ ), Kojima, 2009 ( ▲ ) and Liu and Lange, 2003 ( ★ ).

	$n_{1100K}$ (g/cm3)		$\alpha$ (10-4g/cm3/K)
	MD	Exp	MD	Exp
\ceLi2CO3	1.796	1.792 ◆	4.251	3.729 ◆
		1.793 ▲		3.385 ▲
		1.795 ★		3.459 ★
\ceNa2CO3	1.991	1.986 ◆	5.515	4.487◆
		1.980 ▲		4.696 ▲
		1.988 ★		4.511 ★
\ceK2CO3	1.931	1.928 ◆	4.222	4.421 ◆
		1.925 ▲		3.563 ▲
		1.939 ★		4.615 ★

Table 3. Table 3: Parameters of the third-order Birch Murnaghan equation of state and T − P 𝑇 𝑃 T-P domain of validity for each pure end-member composition. See text and equations ( 4 ), ( 5 ) and ( 6 ) for details.

	$T_{ref}$	$n_{ref}^0$	${\alpha }_0$	${\alpha }_1$	$K_{ref}^0$	$b_1$	$b_2$	$K^{\prime 0}$	$T$	$P$
	(K)	(${\text{g/cm}}^3$)	(K-1)	(K-2)	(GPa)	(K-1)	(K-2)		(K)	(GPa)
Li2CO3	1100	1.80	-2.57$\times$10-4	2.41$\times$10-8	12.56	9.4$\times$10-4	8.5$\times$10-7	8.0	1100-1923	0-6
Na2CO3	1140	1.97	-2.64$\times$10-4	-4.61$\times$10-8	9.01	11.4$\times$10-4	12.7$\times$10-7	8.5	1100-2073	0-15
K2CO3	1190	1.89	-2.56$\times$10-4	-10.05$\times$10-8	5.96	12.0$\times$10-4	12.2$\times$10-7	7.6	1100-2073	0-15

Table 4. Table 4: Density and bulk modulus of the LKe and LNKe mixtures at T = 𝑇 absent T= 1100 K calculated either by MD simulation ( n M ​ D 0 subscript superscript 𝑛 0 𝑀 𝐷 n^{0}_{MD} , K M ​ D 0 subscript superscript 𝐾 0 𝑀 𝐷 K^{0}_{MD} ) or by using ideal mixing rules ( n M ​ D , m ​ i ​ x 0 subscript superscript 𝑛 0 𝑀 𝐷 𝑚 𝑖 𝑥 n^{0}_{MD,\,mix} , K M ​ D , m ​ i ​ x 0 subscript superscript 𝐾 0 𝑀 𝐷 𝑚 𝑖 𝑥 K^{0}_{MD,\,mix} ). Experimental data n E ​ x ​ p 0 subscript superscript 𝑛 0 𝐸 𝑥 𝑝 n^{0}_{Exp} , K E ​ x ​ p 0 subscript superscript 𝐾 0 𝐸 𝑥 𝑝 K^{0}_{Exp} are given for comparison. Liu and Lange ( 2003 ) Results issued from the MD model of Corradini et al. Corradini, Coudert, and Vuilleumier ( 2016 ) ( † ) are given in parenthesis.

	LKe		LNKe
$n_{Exp}^0$ (g/cm3)	1.88		1.91
$n_{MD}^0$ (g/cm3)	1.88	(1.85†)	1.91
$n_{MD,mix}^0$ (g/cm3)	1.87		1.91
$K_{Exp}^0$ (GPa)	8.9		9.1
$K_{MD}^0$ (GPa)	8.6	(6.81†)	8.9
$K_{MD,mix}^0$ (GPa)	8.6		8.8

Table 5. Table 5: First neighbors in molten carbonates at 1100 K: distance r 𝑟 r and coordination numbers N c subscript 𝑁 𝑐 N_{c} (with X=Li, Na or K and O and C are the oxygen and the carbon atom of \ce CO3, respectively). The first line reports values corresponding to the first maximum of the pair distribution function, the second line values evaluated at its first minimum.

	$T$	$P$	$n$	$r^{\text{X-X}}$	$N_c^{\text{X-X}}$	$r^{\text{C-C}}$	$N_c^{\text{C-C}}$	$r^{\text{X-C}}$	$N_c^{\text{X-C}}$	$N_c^{\text{C-X}}$	$r^{\text{X-O}}$	$N_c^{\text{X-O}}$	$N_c^{\text{O-X}}$
	(K)	(GPa)	(g/cm3)	(Å)		(Å)		(Å)			(Å)
\ceLi2CO3	1100	0	1.80	3.15	2.5	4.63	5.2	2.65	1.1	1.9	2.04	1.8	1.0
				4.75	12.4	6.22	14.3	4.03	4.5	9.0	3.03	5.9	3.9
\ceNa2CO3	1100	0	1.94	3.56	2.8	5.05	5.0	2.89	1.0	1.9	2.40	1.6	0.9
				5.38	13.7	6.90	15.0	4.48	4.7	9.4	3.48	6.8	4.6
\ceK2CO3	1100	0	1.89	3.94	2.9	5.40	4.8	3.25	1.1	2.3	2.73	1.5	1.2
				6.20	21.6	7.70	15.5	5.05	4.8	9.5	3.82	7.2	5.0

Table 6. Table 6: Parameters of the Arrhenius activation law ( 14 ) obtained by the interpolation of all MD simulation points (with X = Li, Na or K). P 𝑃 P is in GPa.

	$D_{0,X}$	$E_a^{D_X}$	$V_a^{D_X}$	$D_{0,C{O}_3}$	$E_a^{D_{C{O}_3}}$	$V_a^{D_{C{O}_3}}$
	(${10}^{-9}$m2/s)	(kJ/mol)	(cm3/mol)	(${10}^{-9}$m2/s)	(kJ/mol)	(cm3/mol)
Li2CO3	176	32	2.59-0.12P+0.0037P2	54	33	2.4-0.18P-0.0145P2
Na2CO3	171	37	4.2-0.25P+0.0081P2	82	39	4.0-0.24P+0.0079P2
K2CO3	135	34	5.2-0.43P+0.0179P2	78	38	5.3-0.48P+0.0210P2

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

	${\sigma }_0$	$E_a^{\sigma }$	$V_a^{\sigma }$	${\eta }_0$	$E_a^{\eta }$	$V_a^{\eta }$
	(S/m)	(kJ/mol)	(${\text{cm}}^3$/mol)	($\text{Pa}\cdot \text{s}$)	(kJ/mol)	(${\text{cm}}^3$/mol)
Li2CO3	3295.3	18	1.1+0.166P-0.0148P2	2.7$\times$10-4	24	2.2+0.01P
Na2CO3	2143.7	22	3.6-0.208P+0.0072P2	2.5$\times$10-4	27	3.5-0.10P
K2CO3	1101.6	17	3.9-0.241P+0.0076P2	2.2$\times$10-4	27	4.7-0.20P

Table 8. Table 8: Haven ratio and its ionic contributions. Note that the error bar on these quantities is of the order of the error bar on the electrical conductivity: ∼ similar-to \sim 5–10 %.

	$T$(K)	$P$ (GPa)	$\sigma$ (S/m)	$H$	$H_X$	$H_{{\mathrm{CO}}_3}$	$H_{{\mathrm{XCO}}_3}$
\ceLi2CO3	1100	0	430	1.04	0.04	-0.32	0.32
	1500	4.0	470	1.06	0.08	-0.34	0.32
\ceNa2CO3	1140	0	224	0.99	-0.24	-0.25	0.48
	1500	4.0	154	0.80	-0.28	-0.31	0.39
\ceK2CO3	1200	0	196	0.94	-0.38	-0.14	0.46
	1500	4.0	97	0.47	-0.38	-0.24	0.35

Equations

$$V_{\text{CO}}^{intra-{\mathrm{CO}}_3}(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 {ϵ}_0{r}_{\text{CO}}.$$

$$V_{\text{OO}}^{intra-{\mathrm{CO}}_3}(r_{\text{OO}})=A_{\text{OO}}^{intra-{\mathrm{CO}}_3}\exp (-r_{\text{OO}}/{\rho }_{\text{OO}}^{intra-{\mathrm{CO}}_3}).$$

$$V_{ij}^{inter}(r_{ij})=A_{ij}\exp (-r_{ij}/{\rho }_{ij})-C_{ij}/r_{ij}^6+q_i{q}_j/4\pi {ϵ}_0{r}_{ij},$$

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

$$n_T^0=n_{ref}^0{e}^{{\int }_{T_{ref}}^T-({\alpha }_0+{\alpha }_{1}T)\mathrm{d}T},$$

$$K_T^0=\frac{K_{ref}^0}{1+b_1(T-T_{ref})+b_2(T-T_{ref}{)}^2}.$$

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

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

$$\centering\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,\@add@centering$$

$$D_s$$

$$\sigma$$

$$\eta$$

$$\begin{split}\Pi_{\alpha\beta}(t)=\frac{1}{V}\sum_{i=1}^{N}&\Big{(}m_{i}v_{i\alpha}(t)v_{i\beta}(t)\\ &+\sum_{j>i}r_{ij\alpha}(t).F_{ij\beta}(t)\Big{)}\enspace,\end{split}$$

$$D_s(P,T)$$

$$\sigma (P,T)$$

$$\eta (P,T)$$

$$\begin{split}\sigma^{NE}=&\frac{e^{2}}{k_{\rm B}TV}\sum_{s}N_{s}z_{s}^{2}D_{s}\\ =&\frac{e^{2}}{k_{\rm B}TV}\lim_{t\to\infty}\frac{1}{6t}\sum_{s}z_{s}^{2}\sum_{i=1}^{N_{s}}\langle\big{(}\vec{\Delta}_{i}^{(s)}(t)\big{)}^{2}\rangle\enspace,\end{split}$$

$$\begin{split}\sigma=&\frac{e^{2}}{k_{\rm B}TV}\lim_{t\to\infty}\frac{1}{6t}\langle\Big{(}\sum_{s}\sum_{i=1}^{N_{s}}z_{s}\vec{\Delta}_{i}^{(s)}(t)\Big{)}^{2}\rangle\\ =&\frac{e^{2}}{k_{\rm B}TV}\lim_{t\to\infty}\frac{1}{6t}\Big{\{}\sum_{s}z_{s}^{2}\Big{(}\sum_{i=1}^{N_{s}}\langle\big{(}\vec{\Delta}_{i}^{(s)}(t)\big{)}^{2}\rangle\\ &+\sum_{i=1}^{N_{s}}\sum_{j\neq i}^{N_{s}}\langle\vec{\Delta}_{i}^{(s)}(t).\vec{\Delta}_{j}^{(s)}(t)\rangle\Big{)}\\ &+\sum_{s}\sum_{s^{\prime}\neq s}z_{s}z_{s^{\prime}}\Big{(}\sum_{i=1}^{N_{s}}\sum_{j=1}^{N_{s^{\prime}}}\langle\vec{\Delta}_{i}^{(s)}(t).\vec{\Delta}_{j}^{(s^{\prime})}(t)\rangle\Big{)}\Big{\}}\enspace,\end{split}$$

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

$$H_s=\underset{t\to \infty }{\lim }\frac{1}{6t}\frac{z_s^2}{{\sum }_s{N}_s{z}_s^2{D}_s}\sum _{i=1}^{N_s}\sum _{j\ne i}^{N_s}⟨{\vec{\Delta }}_i^{(s)}(t).{\vec{\Delta }}_j^{(s)}(t)⟩,$$

$$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}\sum _{i=1}^{N_s}\sum _{j=1}^{N_{s^{\prime }}}⟨{\vec{\Delta }}_i^{(s)}(t).{\vec{\Delta }}_j^{(s^{\prime })}(t)⟩.$$

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

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

Full text

To model molten carbonates by atomistic simulations, we have developed an optimized classical force field based on experimental data of the literature and on the liquid structure issued from ab initio molecular dynamics simulations performed by ourselves. In implementing this force field into a molecular dynamics simulation code, we have evaluated the thermodynamics (equation of state and surface tension), the microscopic liquid structure and the transport properties (diffusion coefficients, electrical conductivity and viscosity) of molten alkali carbonates (\ceLi2CO3, \ceNa2CO3, \ceK2CO3 and some of their binary and ternary mixtures) from the melting point up to the thermodynamic conditions prevailing in the Earth’s upper mantle ($\sim$ 1100–2100 K, 0–15 GPa). Our results are in very good agreement with the data available in the literature. To our knowledge a reliable molecular model for molten alkali carbonates covering such a large domain of thermodynamic conditions, chemical compositions and physicochemical properties has never been published yet.

molten carbonates, molecular dynamics simulations, equation of state, microscopic structure, transport properties, surface tension

I Introduction

Carbonate melts receive an increasing interest in both fundamental and applied fields either in the view of developing technological devices (e.g. fuel cell technology)Chery, Lair, and Cassir (2015a, b); Cassir, McPhail, and Moreno (2012); Cassir, Ringuedé, and Lair (2013); Lair et al. (2012) or for providing a better understanding of a number of geochemical processes (e.g. role of molten carbonates in the geodynamics of the Earth’s mantle).Gaillard et al. (2008); Dasgupta and Hirschmann (2010); Dasgupta (2013); Hammouda and Keshav (2015) The study of the physicochemical properties (thermodynamics, structure and transport) of, mainly alkali, molten carbonates, under atmospheric pressure has been initiated by the electrochemistry community a few decades ago.Janz (1988) These systems were depicted as highly conducting liquids of low viscosity. From a geosciences viewpoint the focus is more recent and mainly set on alkaline-earth compositions under high pressures, whose properties are so far poorly known.Gaillard et al. (2008); Sifré, Hashim, and Gaillard (2015); Kono et al. (2014) Indeed studies of these liquids at high $T-P$ require very specific devices adapted to the measurement a single observable. On the contrary, in a molecular dynamics simulation a number of properties can be obtained at once, even at extreme thermodynamic conditions. However molecular dynamics simulations require the knowledge of atomic interactions in the system under consideration. These interactions (the force field) can be deduced from electronic structure calculations (ab initio) or based on a more empirical approach (classical molecular dynamics or simply referred to as MD in the following). The former method, although a priori more accurate than the latter one is computationally much more demanding (by several orders of magnitude) and is thus inappropriate when the aim of the study is to span a large range of chemical composition and thermodynamic conditions. On the other hand, the relevance of a MD calculation relies on the quality of the implemented force field (FF) whose adjustment is system-specific. Since the first attempts by Janssen and TissenJanssen and Tissen (1990); Tissen and Janssen (1990, 1994) to model molten carbonates by two-body interactions, several modifications have been suggested. (Habasaki, 1990; Koishi et al., 2000; Costa, 2008; Ottochian et al., 2016; Wilding et al., 2016) However the FF of Janssen and Tissen (1990) has some inherent defects. For instance the calculated pressure is shifted by $\sim$15 kbar,Ottochian et al. (2016) at the 1-bar experimental density. A decisive improvement towards an accurate description of the properties of molten carbonates in quantitative agreement with the experiments came with the force fields (FF) of Vuilleumier et al. (2014) and Corradini *et al.*Corradini, Coudert, and Vuilleumier (2016). However these models are restricted in composition, to pure \ceCaCO3 and to the \ceLi2CO3–\ceK2CO3 (62:38 mol%) eutectic mixture respectively. For the sake of conciseness and clarity of the comparisons with the experimental literature, alkaline-earth compositions will be investigated in a companion paper.Desmaele (2017) Here we present MD simulation results that sketch a systematic and consistent model of molten carbonates in the \ceLi2CO3–\ceNa2CO3–\ceK2CO3 system over a large $T-P$ range ($\sim$ 1100–2100 K, 0–15 GPa). After presenting the details of the simulations, we discuss the thermodynamic properties (equation of state and surface tension), the liquid structure and the transport properties (diffusion coefficients, electrical conductivity and viscosity).

II Computational details

II.1 Ab initio molecular dynamics (AIMD)

The ab initio molecular dynamics (AIMD) simulations were based on the density functional theory (DFT) within the Born-Oppenheimer framework. They were performed thanks to the free QUICKSTEP/CP2K softwareVandeVondele et al. (2005) that applies a hybrid Gaussian plane-wave method.Lippert, Hutter, and Parrinello (1997) For valence electrons, we used a double-zeta valence plus polarization (DZVP) basis set optimized for molecules,VandeVondele and Hutter (2007) while core electrons were replaced 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 functionalBecke (1988); Lee, Yang, and Parr (1988) using a semi-empirical D3 dispersion correction scheme with a 30 Å  cutoff.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) Timestep was fixed to 0.5 fs, and simulation run was of the order of 13–48 ps.

The starting configurations were issued from classical MD runs using a guess force field. \ceNa2CO3 and \ceK2CO3 simulations were performed at state points corresponding to their atmospheric melting point,van Groos (1990); Klement and Cohen (1975) that is respectively (1140 K, 1.97 g/cm3) and (1190 K, 1.90 g/cm3), where densities were set to their experimental values.Liu and Lange (2003) The pressure calculated along both runs is 0 $\pm$ 0.5 GPa, as expected. To investigate the high pressure behavior, a third simulation was performed for \ceK2CO3 at 1773 K and 2.41 g/cm3 leading to a calculated pressure of 6 $\pm$ 1.5 GPa. Table S1 summarizes the simulation details of the three AIMD simulations. Note that the AIMD simulation of Corradini *et al.*Corradini, Coudert, and Vuilleumier (2016) for the \ceLi2CO3–\ceK2CO3 eutectic mixture uses similar methodological features, and hence it will be considered as a benchmark for Li-bearing compositions investigated here by MD.

II.2 Classical molecular dynamics (MD)

Classical MD simulations were carried out using the free DL_Poly_2 software. Smith and Forester (1996) The timestep was set to 1 fs and the studied systems were composed of $N\simeq$ 2000 atoms. To evaluate the equation of state (EoS) as discussed in section IV.1, calculations in the $NPT$ ensemble were performed. The total simulation time was 0.9 ns, including a 0.5 ns equilibration with a Nosé-Hoover thermostat, to reach a statistical accuracy on the density of about $\Delta n/n\leqslant\pm 1~{}\%$. For transport coefficients, production runs of 10 to 30 ns were performed in the $NVE$ ensemble. Structure data were extracted from the same simulations. Note that for a better comparison with the AIMD-generated pair distribution functions (see section III), some MD simulations were also specifically run on systems having the same density and number of atoms.

III Force field

The empirical force field is derived from an interaction potential summing the pair interactions between the atoms of a same carbonate molecule (intramolecular potentials $V_{ij}^{intra}$) and those between atom pairs belonging to different ions (intermolecular potentials $V_{ij}^{inter}$).

The carbonate molecule is featured as flexible and non-dissociative by introducing the following pair potential between C and O atoms:

[TABLE]

The force constant ${k_{\text{CO}}}=6118.17$ kJ/mol and the harmonic equilibrium distance ${r_{0,\text{CO}}}=$ 1.30 Å  are adjusted so as to obtain a mean C-O distance of 1.29 Å.Antao and Hassan (2010); Vuilleumier et al. (2014) Moreover the coulombic interactions are implemented in order to prevent the atomization of the molecule at high pressure. The oxygen atoms of a same molecule repel one another via the pair potential:

[TABLE]

The latter one ensures on average the planarity of the carbonate ion, where ${A_{\text{OO}}^{intra-\rm{CO}_{3}}}=2611707.2$ kJ/mol and ${\rho_{\text{OO}}^{intra-\rm{CO}_{3}}}=$ 0.22 Å. The intermolecular interactions are ruled by the pair potential:

[TABLE]

with $i,j$= O, C, Li, Na, K. The van der Waals part (first two terms) of this potential is set to zero between cations, thus enabling a straightforward transferability of the force field to any carbonate mixture in the Li–Na–K system. Parameters of the FF are summarized in Table 1. While the charges were constrained by the FF of Vuilleumier et al. (2014), the van der Waals interactions have been benchmarked on experimental volumetric studies and on AIMD structure data (note that dispersion interactions have little effect on the shape of the pair distribution functions, but affects strongly the density). Thus atmospheric density and compressibilityLiu and Lange (2003) are accurately reproduced by our model as discussed in section IV. The suitability of the force field was further constrained by fitting at best the pair distribution functions of AIMD simulations, obtained either by Corradini *et al.*Corradini, Coudert, and Vuilleumier (2016) for the \ceLi2CO3–\ceK2CO3 eutectic mixture of ratio 0.62:0.38 mol% (LKe), or specifically calculated by ourselves for \ceNa2CO3 and \ceK2CO3. The agreement between MD simulations and AIMD structure data is excellent (Figure S1). Notice that for LKe in particular, the fitting of the AIMD curves by MD is at least as good, as when using Corradini et al.’s FF (Figure 1 of reference Corradini, Coudert, and Vuilleumier, 2016). The microscopic structure is thoroughly described in section IV.2.

To be complete, and although the present FF is not aimed at studying the crystalline phase, we checked its ability to reproduce the experimentally determined structures of the monoclinic crystals of \ceLi2CO3, \ceNa2CO3 and \ceK2CO3 at room conditions. The densities calculated by anisotropic relaxation ($NST$ ensemble) of these structures are 2.12 g/cm3 for \ceLi2CO3, 2.50 g/cm3 for \ceNa2CO3 and 2.45 g/cm3 for \ceK2CO3 to be compared to 2.10, 2.54 and 2.44 g/cm3, respectivelyIdemoto et al. (1998); Arakcheeva et al. (2010); Gatehouse and Lloyd (1973) (see also Table S2 for lattice parameters).

IV Thermodynamics

IV.1 Equation of state

The density of molten alkali carbonates at atmospheric pressure has been reported by several experimental studies Janz and Lorenz (1961); Spedding (1970); Janz (1988); Hong-Min et al. (1991); Kojima (2009); Liu and Lange (2003); Kojima et al. (2003, 2008); Kojima (2009). For the pure end-members, \ceLi2CO3, \ceNa2CO3 and \ceK2CO3 a global consistency between these studies is found. Using the present FF, MD fits very well the above studies considering the small experimental uncertainties and the very limited temperature range ($\sim$ 100 K) investigated at 1 bar (Table 2). In this context notice that the model of Corradini *et al.*Corradini, Coudert, and Vuilleumier (2016) recovers reasonably the 1-bar density of \ceK2CO3, but it doesn’t perform as well for pure \ceLi2CO3, the deviations reaching $-$5 %.

Despite the geochemical interest for assessing the buoyancy of carbonated fluids in the mantle, no systematic measurements of liquid density at high pressure have been reported yet, except for the study of Dobson et al.Dobson et al. (1996), for \ceK2CO3 at 4 GPa. Besides the study of Dobson et al. has since been challenged by Liu *et al.*Liu, Tenner, and Lange (2007). These latter authors provided estimates of high pressure volumetric behavior for \ceK2CO3,Liu and Lange (2003); Liu, Tenner, and Lange (2007) by using the third-order Birch Murnaghan equation of state (BMEoS)Birch (1947) that is well-known for its accuracy in the modeling of the compressibility of molten and crystalline silicates Rigden, Ahrens, and Stolper (1989); Guillot and Sator (2007). It writes:

[TABLE]

where $n_{T}^{0}$ is the atmospheric density at temperature T, $K^{0}_{T}$ the bulk modulus and $K^{\prime 0}$ its pressure derivative at 1 bar. Such an equation expressing the density of the carbonated liquid at any thermodynamic condition is appealing, although $K^{\prime 0}$ has so far only been derived by indirect thermodynamic considerations (see below).

To determine the BMEoS of liquid \ceLi2CO3, \ceNa2CO3 and \ceK2CO3 from MD, their density was calculated along several isotherms at different pressures. The values of $n_{T}^{0}$ and $K^{0}(T)$ were fitted by using empirical temperature dependencies:

[TABLE]

and

[TABLE]

The reference temperature $T_{ref}$ is defined for numerical purposes and set to the atmospheric melting point for each composition (Table 3), $n_{ref}^{0}$ and $K_{ref}^{0}$ are the density and bulk modulus at 1-bar and $T_{ref}$. Values of the parameters $\alpha_{0}$, $\alpha_{1}$, $b_{1}$ and $b_{2}$ are given in Table 3.

The excellent agreement between calculated points and the fitted isotherms (Figures 1 and S2) illustrates the ability of the BMEoS to describe carbonate melts, at least on the investigated pressure-temperature range. Furthermore, the obtained BMEoS for \ceLi2CO3, \ceNa2CO3 and \ceK2CO3 match not only the experimental values at 1 bar, but also those of our AIMD calculations. However, AIMD simulations are performed at constant density with a relatively small system size, and display large pressure fluctuations which yield uncertainties, typically of the order of $\sim$ 1 GPa (inset of Figure 1).

The 1-bar bulk modulus $K^{0}_{T}$ (inverse of the isothermal compressibility) has been experimentally determined by sound-speed measurements for a number of mixtures.O’Leary, Lange, and Ai (2015) It is found in excellent accordance with the ones resulting from the parameterization of the BMEoS using MD (e.g. in Table 4). $K^{0}_{T}$ consistently decreases with the cationic radius (Table 3), unlike its P-dependence, $K^{\prime 0}$, that is almost invariant with chemical composition. Using a fusion curve analysis, Wang et al. (2016) have constrained $K^{\prime 0}$ for liquid K2CO3, at pressure below 5 GPa, to 14.4 instead of 7.6 by MD. This disagreement can be explained by the indirect method used by the authors to obtain $K^{\prime 0}$ from the Clapeyron curve. In this method (i) many thermodynamic data (which introduce uncertainties) are required, and (ii) the melting curve is not accurately known in the light of the differences between the studies of Klement and CohenKlement and Cohen (1975), Wang et al. (2016) and Li.Li (2015)

For carbonate mixtures, discrepancies on atmospheric density measurements can reach up to 5 %.Janz and Lorenz (1961); Spedding (1970); Janz (1988); Liu and Lange (2003); Kojima et al. (2003, 2008) The values from the study of Liu *et al.*Liu and Lange (2003) seem to be of first relevance since the authors have proven their reproducibility, have accounted for surface tension corrections and have tested their apparatus on standard liquids. The comparison of our MD results with their study shows the excellent transferability of our force field to the determination of the density of mixtures, as illustrated in Table 4 on some selected compositions, namely the eutectic ternary mixture (\ceLi2CO3, \ceNa2CO3 and \ceK2CO3 in proportions 43.5:31.5:25 mol%,(Maru, 1984) LNKe) and a binary eutectic mixture (\ceLi2CO3 and \ceK2CO3 in proportions 62:38 mol%, LKe). Note that the density of mixtures follows an ideal mixing rule as it was reported by Liu *et al.*Liu and Lange (2003):

[TABLE]

where $x_{i}$ is the molar fraction of species $i$ of molar mass ${M}_{i}$ and molar volume $\bar{V}_{i}$.

For the bulk modulus, an ideal mixing rule,

[TABLE]

is also verified (Table 4).

Therefore the use of the above mixing rules for the density and the bulk modulus, combined with the fact that $K^{\prime 0}$ varies very little for the three end-members, leads to an accurate estimation of the density of mixtures in the system \ceLi2CO3–\ceNa2CO3–\ceK2CO3 over a large range of $T-P$ conditions. For instance, considering the eutectic ternary mixture at 1500 K and 3 GPa, direct simulation gives a density of 2.15 g/cm3, while the use of equations (7) and (8), with $K^{\prime 0}$ in the range $7.6-8.5$ gives 2.16 g/cm3 after using equation (4).

IV.2 Evolution of the liquid structure with temperature and pressure

Despite their liquid nature, related to a long-range disorder, molten carbonates show a residual local ordering coming from the competition between ionic, dispersive and steric forces. Among these, the coulombic forces are essential as they impose an alternation of charges of opposite signs, as qualitatively confirmed by X-ray and neutron diffraction studies of molten carbonates.Zarzycki (1961); Kohara et al. (1998) Thus the pair distribution functions (PDF) of X$-$C and C$-$C (where C is the carbon atom of \ceCO3 and X=Li, Na or K) are nearly in opposition of phase (Figure S1). However the structure of molten carbonates is complicated further, in comparison to simple ionic liquids, by the peculiar shape of the carbonate ion and by the anion-cation asymmetry of charge. The anion is surrounded by cations that subsequently attract carbonate anions into a second coordination shell. In turn the anions encircle the cation by pointing one or two oxygen atoms in its direction (see Table 5). Actually the number of pointing O of each carbonate is barely sensitive to the size of the cation ($\sim$1.5 in each case, Table 5). The number density progresses as follows at 1 bar and 1100 K: $\rho_{\ce{Li2CO3}}=$ 24.4 mol/L $>$ $\rho_{\ce{Na2CO3}}=$ 18.7 mol/L $>$ $\rho_{\ce{K2CO3}}=$ 14.0 mol/L. Therefore some structural features can be described as a simple homothetic transformation upon volumetric shrinkage from \ceK2CO3 to \ceLi2CO3 (e.g. C$-$C pair on Figure 2). Interestingly, this simple behavior leads to a pseudo-ideal mixing rule of the cation-carbonate pair distribution functions. As a matter of fact the X$-$C PDFs associated with the pure end-members are quasi-identical to the ones observed in the corresponding mixtures (not shown).

According to the parametrization procedure we used, the PDFs from MD compare very well with the ones from AIMD (Figures S1a, S1b and S1c). Moreover the quality of the MD description is preserved at high $T-P$ as illustrated by the example of \ceK2CO3 (Figure S1d). Actually the agreement between MD and AIMD is even better under pressure, except for the cation-cation (X$-$X) pair. With regard to the structural changes occurring with increasing pressure and temperature, it is noteworthy that on the $T-P$ range under study the local structure in the first shell depends very little on the temperature. For example, for \ceK2CO3 at 4 GPa, the cation-oxygen $N_{c}^{\text{X-O}}$ coordination number is 8.3 at 1500 K and 8.0 at 2073 K. Consequently the evolution of the coordination numbers $N_{c}^{\text{C-C}}$, $N_{c}^{\text{X-C}}$, $N_{c}^{\text{X-O}}$ at increasing pressure and temperature can thus be assigned essentially to pressure (Figure 3). Surprisingly for all compositions, $N_{c}^{\text{X-O}}$ evolves very little with pressure, it increases somewhat before levelling off. The ratio $N_{c}^{\text{X-O}}$/$N_{c}^{\text{X-C}}$ plotted in the inset of Figure 3 is indicative of the orientation of the carbonate ions around the cations (there are 1.3 oxygen atoms pointing towards Li, $\sim$1.5-1.6 towards Na and K) and is virtually independent of pressure. Although the pressure evolution of the coordination number of the C$-$C pair is a bit steeper, there is nothing remarkable in this evolution. This is at variance with the hypothesized reminiscence in the liquid structure of crystalline polymorphs at the corresponding pressure.Hudspeth, Sanloup, and Kono (2018) The absence of a significant structural rearrangement of molten carbonates upon pressure should translate in smoothly varying transport coefficients over this pressure range.

Wilding et al. (2016) have studied the structure of \ceNa2CO3 at 1 bar and from 1100 to 1750 K, using MD and a FF adapted from Tissen and Janssen (1990). The authors assumed the formation of a low-dimensional network made of short carbonate-cation-carbonate chains argued by the presence of a shouldering in the first peak of $g_{\text{C-C}}(r)$ at 2.3 Å. Such a shouldering does not appear in our MD study, neither in the pair distribution function generated by AIMD. We believe that the aforementioned structural feature is due to the inaccuracy of the Janssen and Tissen (1990) force field, even when supplemented by an harmonic intramolecular potential for \ceCO3, as used by Wilding et al. (2016)

IV.3 Surface Tension

Surface tension can be evaluated by explicitly simulating the coexistence of the liquid with its vapor along a plane interface,Alejandre, Tildesley, and Chapela (1995); Aguado and Madden (2002); Salanne et al. (2007) and by calculating the long time limit of the average of some components of the stress tensor. More precisely, assuming that the $z$ axis is perpendicular to the interface, the surface tension is given by:Kirkwood and Buff (1949)

[TABLE]

where $L_{z}$ is the length of the simulation box and $\Pi_{xx}$, $\Pi_{yy}$ and $\Pi_{zz}$ are the diagonal components of the stress tensor (see equation (13) in next section). We studied several compositions at different temperatures, but for the sake of clarity, we chose to present the result of only a few of them to picture the main trends (Figure 4). Our model recovers the right order of magnitude of experimental measurements of the surface tension, but seems to overestimate it systematically by 10–20 %. Moreover the temperature dependencies are satisfying although one notices an overestimation of the activation energy in the case of the LNKe mixture. The hierarchy between the different compositions is respected: the bigger the cation, the smaller the surface tension. Thus the surface tension decreases when the asymmetry of size between the cation and the anion decreases.Gonzàlez-Melchor et al. (2005)

As mentioned above, the overestimation of the MD-calculated surface tension seems to be systematic (of the order of +10–15 % compared with Kojima (2009)). As our FF describes with accuracy the structure and the density as well as other properties (EoS, viscosity) which depend also strongly on cohesive forcesMarchand et al. (2011) (see next section), this observation is somehow a surprise and raises concerns on possible systematic errors due to the simulation parameters. Some authors have pointed out that away from the interface the bulk, instead of showing isotropic properties, could wrongly contribute a non-zero amount to equation (9).Sanmartín Pensado, Malfreyt, and Pádua (2009) We tested the possible effect of this anomalous contribution in the case of the \ceLi2CO3–\ceNa2CO3–\ceK2CO3 equimolar mixture at 900 K, by calculating the surface tension with a liquid depth increased by a factor 7 ($z$ dimension from $\sim$ 29 to 200 Å), all else being equal. This resulted in a value of $\gamma$ that is almost unchanged or maybe slightly decreased (from 239 $\pm$ 3 to 232 $\pm$ 6), which suggest a weak contribution of the bulk. We have tested several other parameters of the simulation which could affect the calculation of the surface tension namely, the size of the simulation box, the truncation cutoff of the van der Waals potential, and the parameters of the Ewald summation.Gonzàlez-Melchor et al. (2005) Although the calculated value changes a little with these parameters, the variations seem too weak to explain the discrepancy with the experimental data (calculated values differ by less than 5 %). Larger scaling tests with much greater simulations boxes would be needed to completely rule out the role of finite size effect. However, it could be that our FF is less suitable to describe interfacial (rather than bulk) properties as it possibly fails to render some polarization effects. In any case, the discrepancies with the experiments we obtained are similar in magnitude with previous MD studies of ionic liquids.Aguado and Madden (2002); Salanne et al. (2007) and this first attempt of simulating the surface tension of molten carbonates by MD leads to good qualitative predictions.

V Transport properties

Aiming at a systematic description of transport properties as a function of chemical composition, temperature and pressure, a large number of simulations were performed at different thermodynamic conditions for pure \ceLi2CO3, \ceNa2CO3 and \ceK2CO3, as well as for some binary and ternary compositions. Each run was carried out in the $NVE$ ensemble for a duration of 10–20 ns, which is sufficient to investigate the diffusive regime and to calculate accurately the self-diffusion coefficients $D_{s}$ of each chemical species ($s=$ Li, Na, K and CO${}_{3}^{2-}$), the electrical conductivity $\sigma$ and the viscosity $\eta$, defined respectively by:Allen and Tildesley (1989); Hess (2002)

[TABLE]

where the off-diagonal pressure tensor components $\Pi_{\alpha\beta}(t)$, with $\alpha,\beta=x,y,z$ is given by:

[TABLE]

where $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, $e$ the elementary charge, $m_{i}$ the mass of ion $i$ and $\mathbf{r}_{i}(t)$ its position, $v_{i\alpha}(t)$ the component $\alpha$ of its velocity and $F_{ij\beta}(t)$ is the component $\beta$ of the force exerted by ion $j$ on ion $i$, separated by distance $r_{ij}(t)$ at time $t$. The ensemble average \big{\langle}\dots\big{\rangle} is taken over many time origins. The parameter $z_{i}$ is the conduction charge of ion $i$ (i.e. $\sigma$ is effectively an ionic conductivity), chosen here as the formal charge, as it is usual for simple ionic liquids. Adams, McDonald, and Singer (1977) Note that the optimal choice for the conduction charges, either formal or partial charges, is still an open question.(Vuilleumier et al., 2014; Corradini, Coudert, and Vuilleumier, 2016) The estimation of the self-diffusion coefficient $D_{s}$ relies on an average over the $N_{s}$ ions of a specific species $s$ (sum $\sum_{i=1}^{N_{s}}$ in eq (10)), which ensures a small uncertainty of $\sim$1 %. On the contrary, conductivity and viscosity are collective quantities that necessarily require a large number of atoms (and long time runs) to be correctly estimated with an error bar of 5–10 % in this study. For this reason, AIMD estimations are today only crude and suffer from high uncertainties. As a matter of fact, until now only two classical MD studies have provided estimations of the viscosity and electrical conductivity of carbonate melts, namely \ceCaCO3Vuilleumier et al. (2014) and the \ceLi2CO3–\ceK2CO3 eutectic mixture (62:38 mol%) .Corradini, Coudert, and Vuilleumier (2016)

As usually observed for molten salts, the temperature dependence of the transport properties calculated by MD is well reproduced by an Arrhenius activation law. An activation volume $V_{a}^{X}$ is introduced to take into account the pressure dependence of the physical quantity $X=D_{s}$, $\sigma$ and $\eta$, as in references Vuilleumier et al., 2014; Corradini, Coudert, and Vuilleumier, 2016. The following expressions are used:

[TABLE]

where $E_{a}^{X}$ is the activation energy associated to the physical quantity $X$. The values of $X^{0}$, $E_{a}^{X}$ and $V_{a}^{X}$ were determined by fitting the MD data for the three end-members \ceLi2CO3, \ceNa2CO3 and \ceK2CO3 (see Figures S3, S4, S5, S6 and S7). These values are given in Tables 6 and 7, for $D_{s}$, $\sigma$ and $\eta$ and for the pressure and temperature range mentioned in Table 3.

A crucial advantage of MD simulations is that all these transport coefficients are calculated in the same theoretical framework, which allows one to check the validity of approximate relations between electrical conductivity and self-diffusion coefficients (Nernst-Einstein equation) and between self-diffusion coefficients and viscosity (Stokes-Einstein equation). Furthermore, the role played by each chemical species and the correlations between their contributions to the transport properties can be easily investigated, as we shall see below.

V.1 Self-diffusion coefficients

We have first estimated the ion self-diffusion coefficients from our AIMD simulations, which were long enough to reach the diffusive regime. Thus at 1 bar we found $D_{\rm Na}=$ 3.4 $\pm$ 0.3 $10^{-9}$ m2/s and $D_{\rm CO_{3}}=$ 1.5 $\pm$ 0.3 $10^{-9}$ m2/s in \ceNa2CO3 at 1140 K, and $D_{\rm K}=$ 4.0 $\pm$ 0.2 $10^{-9}$ m2/s and $D_{\rm CO_{3}}=$ 1.8$~{}\pm$ 0.3 $10^{-9}$ m2/s in \ceK2CO3 at 1190 K, values which compare well with those obtained from MD (3.5 and 1.4, and 4.6 and 1.8, respectively) a finding which emphasizes the reliability of our empirical force field.

The only diffusivity data in molten alkali carbonates have been provided by Mills and SpeddingSpedding and Mills (1965, 1966); Mills and Spedding (1966) who measured the self-diffusion coefficients of Na+, K+ and CO${}_{3}^{2-}$ ions in pure, binary and ternary alkali carbonate melts, at atmospheric pressure only. Note that the self-diffusion coefficient of Li+ has not been measured yet.

As illustrated in Figure 5, the agreement between MD and experimental results is satisfactory for $D_{\rm K}$ and $D_{\rm Na}$, at least in the pure \ceNa2CO3 melt and in the \ceNa2CO3–\ceK2CO3 eutectic mixture (NKe), but the diffusivity of the CO${}_{3}^{2-}$ anion measured in various melt compositions is systematically higher by a factor of 1.5–3 than the one calculated by MD or by AIMD.

The activation energy (32–39 kJ/mol) depends slightly on the ion species and on the melt composition (Table 6 and Figure 5), which is frequent in molten salts,(Spedding and Mills, 1965) whereas the experimental values are between $\sim$ 40 and 53 kJ/mol, with an error reportedly below 1 kJ/mol. Note that the activation energies ($\sim$ 36 kJ/mol) determined by the MD simulations of Corradini *et al.*Corradini, Coudert, and Vuilleumier (2016) for the molten LKe mixture are close to the ones we found in the same mixture (41 kJ/mol) as well as in the pure end-members (Table 6). However their values of $D_{\rm Li}$ and $D_{\rm K}$ are higher than ours, by a factor of $\sim$1.5.

General trends can be inferred from our MD simulations. By comparing the cation diffusivity in the three end-member compositions at a given temperature, the following hierarchy is found: $D_{\rm Li}>D_{\rm Na}>D_{\rm K}>D_{\rm CO_{3}}$, a behavior also found in experiments. As expected the largest ion CO${}_{3}^{2-}$ has the slowest diffusivity. This hierarchy of the self-diffusion coefficients with the ionic radius is also found in binary and ternary mixtures.

However upon mixing, the diffusion coefficient of an ion species may change with the melt composition. Hence, when diffusivity data are plotted as functions of the \ceLi2CO3 molar content in mixtures $x_{\rm Li_{2}CO_{3}}$, $D_{\rm Na}$ and $D_{\rm CO_{3}}$ display a maximum at the eutectic composition,Spedding and Mills (1965, 1966) as shown in Figure 6. According to Spedding and MillsSpedding and Mills (1965, 1966), this feature could be due to a high level of structural disorder at the eutectic point, which could also explain the low melting temperature for this composition. At variance, our MD simulations do not display such a composition effect. Furthermore the microscopic structure, as expressed by the PDFs, does not seem more disordered at the eutectic point (compare Figure S1a to Figures S1b and S1c). In fact the calculated self-diffusion coefficients of Na+, CO${}_{3}^{2-}$ (and also Li+) increase as the \ceLi2CO3 content increases in the binary mixture (Figure 6). A different situation is observed in the \ceLi2CO3–\ceK2CO3 mixture where $D_{\rm K}$ and $D_{\rm CO_{3}}$ are almost constant and exhibit a very slight minimum at 0.42, whereas $D_{\rm Li}$ is increasing with $x_{\rm Li_{2}CO_{3}}$. A consequence is that $D_{\rm Li}<D_{\rm K}$ when $x_{\rm Li_{2}CO_{3}}<0.42$ (see Figure  6a). Thus, the major effect of mixing is that Na+ and K+ hinder the motion of the light Li+ ion when $x_{\rm Li_{2}CO_{3}}$ is low. Notice that the mobility measurements of Yang et al. (1987) support the above findings. Those trends are also in accordance with previous numerical studies (Tissen and Janssen, 1990, 1994; Koishi et al., 2000; Costa, 2008; Ottochian et al., 2016) based on the force field by Janssen and TissenJanssen and Tissen (1990).

To our knowledge no diffusivity data are available under pressure. Our simulations predict an Arrhenius behavior according to equation (14). Thus the self-diffusion coefficients strongly decrease with pressure (see Figures S3, S4 and S5), a behavior also observed in the \ceCaCO3 melt.Vuilleumier et al., 2014 Concerning the comparison between MD and AIMD calculations we notice that at 6 GPa and 1773 K in the \ceK2CO3 melt, the self-diffusion of K+ and CO${}_{3}^{2-}$ calculated by MD are slightly lower than that obtained from AIMD: $D_{\rm K}=$ 3.6 and 5.0 $\pm$ 0.5 $10^{-9}$ m2/s and $D_{\rm CO_{3}}=$ 1.7 and 3.5 $\pm$ 0.5 $10^{-9}$ m2/s for MD and AIMD respectively, because densities are a little bit different in the two calculations (2.49 and 2.41 g/cm3 for MD and AIMD, respectively).

V.2 Electrical conductivity

The knowledge of the electrical conductivity of alkali carbonates is requisite for their industrial applications as electrolytes in fuel cell devices. Several experimental and numerical studies have be devoted to this issue. The measurements by Kojima et al. are broadly consistent with those reported by Janz (1988) for the end-members and various binaryKojima et al. (2007) and ternaryKojima et al. (2008) mixtures. Gaillard et al. (2008) also measured the electrical conductivity of the equimolar mixtures \ceNa2CO3–\ceK2CO3 (NK) and \ceLi2CO3–\ceNa2CO3–\ceK2CO3 (LNK) at atmospheric pressure. Moreover the electrical conductivity of the LKe eutectic mixture was measured by Lair et al. (2012). The electrical conductivity data at atmospheric pressure are usually fitted by an Arrhenius law (see equation (15)). However some deviations were noticed at low temperature ($<$ 973 K) by Kojima et al., and for this reason they introduced a quadratic function of temperature (in log$T$) to fit their data.

As shown in Figure 7, the electrical conductivity at 1 bar is negatively correlated to the cation size, as mentioned by Kojima et al. (2008) and Sifré *et al.*Sifré, Hashim, and Gaillard (2015) The agreement between MD calculations and the available 1-bar experimental dataKojima (2009) is good for \ceLi2CO3 and \ceK2CO3 (Figure 7), within the error bars of both the simulations and the measurements (expected to be of the order of 10–15 %). Nevertheless the simulations underestimate somewhat the conductivities with respect to the experimental values of KojimaKojima (2009). This discrepancy could be due to experimental inaccuracy. Indeed, the MD results are in very good agreement with the measurements of Gaillard et al. (2008) for the equimolar NK and LK systems and with those of Lair et al.Lair et al. (2012) for the LKe mixture, whose values are systematically below that of KojimaKojima (2009) (Figure 8).

The activation energy of calculated electrical conductivity is between 17 and 22 kJ/mol for the three end-members, which is slightly higher than the values found by Kojima (2009) (15–20 kJ/mol). In the LKe mixture at atmospheric pressure, $E_{a}^{\sigma}$ calculated by Corradini *et al.*Corradini, Coudert, and Vuilleumier (2016) ($E_{a}=20.22$ kJ/mol) is close to the ones measured by Kojima et al. (2007) and Lair et al. (2012) ($E_{a}=21.91$ kJ/mol for Kojima). As for the activation energies determined by Gaillard et al. (2008) in the equimolar \ceNa2CO3–\ceK2CO3 and \ceLi2CO3–\ceNa2CO3–\ceK2CO3 mixtures, they are higher, around 31 kJ/mol. Notice that the electrical conduction is a collective process including ion-ion correlations (see below). Therefore the activation energy associated to the electrical conductivity differs from that associated to the diffusion process (compare Table 7 and Table 6)).

Up to a few kbars, the electrical conductivity depends weakly on pressure (see activation volumes in Table 7), whereas at higher pressures it significantly decreases with pressure (see Figure S6).

Assuming that the ions move independently from each other (the ion-ion cross correlations are then neglected), the ”ideal” electrical conductivity, $\sigma^{NE}$ is given by the Nernst-Einstein equation:

[TABLE]

where $e$ is the elementary charge, $k_{\rm B}$ the Boltzmann constant and $N_{s}$, $z_{s}$ and $D_{s}$ are the number of ions, the conduction charge, equal to the formal charge here, and the self-diffusion coefficient of the chemical species $s=$Li, Na, K, and \ceCO3 respectively. In the right-hand side of equation (17), we made use of the definition (10) of the self-diffusion coefficient by writing $\vec{\Delta}_{i}^{(s)}(t)=\vec{r}_{i}^{(s)}(t)-\vec{r}_{i}^{(s)}(0)$.

The Haven ratio $H=\sigma/\sigma^{NE}$ is defined as the ratio of the exact electrical conductivity $\sigma$, here calculated using the Einstein relation (11) to the Nernst-Einstein conductivity (17). Note that the same electric charges are used (here, the formal ones) for evaluating $\sigma$ and $\sigma^{NE}$, therefore the value of $H$ is independent of the charges assigned to the ions.

The value of $H$ is usually less than 1 for molten salts, which is explained by the large contribution of the anticorrelations between ions of the same charge that tends to decrease the ideal conductivity $\sigma^{NE}$Kashyap et al. (2011). Simulations studies showed that $H\simeq 0.81$ for CaCO3 melt,Vuilleumier et al. (2014) $H\sim 0.7$ for \ceLi2CO3–\ceNa2CO3 mixturesOttochian et al. (2016) and $H=$0.5-1.0 for \ceLi2CO3–\ceK2CO3 mixturesCosta (2008). However, the Haven ratio was found greater than 1 ($H=1.8$) for LKe.Corradini, Coudert, and Vuilleumier (2016) Note that $H$ may depend on pressure and temperature. In this work, we find a Haven ratio which is slightly greater than 1 for \ceLi2CO3 (e.g. $H=1.04$ at the melting point) and lower than 1 for \ceNa2CO3 and \ceK2CO3 (e.g. $H=0.99$ and 0.94, respectively, at the melting point). This ratio seems to depend slightly on temperature, and more appreciably on pressure, especially for \ceNa2CO3 and \ceK2CO3. For instance, at 4 GPa and 1500 K (near the melting point) the $H$ ratio of \ceK2CO3 drops to 0.47 (Table 8). The amplitude of $H$ is often considered to be indicative of the degree of dynamic correlations between ions. So a value close to 1 should confirm the independence hypothesis of the Nernst-Einstein approximation. Nevertheless, we show below that such correlations are present in the melt even when $H\simeq 1$ and that the variations of $H$ with chemical composition and thermodynamic conditions are due to subtle changes in dynamic correlations between ionic motions.

To better understand the role of ion-ion correlations it is suitable to develop the exact Einstein expression of the conductivity given by equation (11). Note that in this expression the charge displacement is squared as opposed to equation (10) that is the sum of each squared ionic displacement. After a few algebra equation (11) leads to the following expression

[TABLE]

where $N_{s}$ is the number of ions of species $s$, $z_{s}$ is their charge and $\Delta_{i}^{(s)}(t)$ is the displacement of the ion $i$ of species $s$. The first term in the RHS of the last equation is the Nernst-Einstein conductivity (17). We then get the expression of the Haven ratio in terms of the correlations between ions:

[TABLE]

where $H_{s}$ expresses the average cross correlations (through a scalar product) between the displacements of ions of species $s$:

[TABLE]

and $H_{ss^{\prime}}$ expresses the average cross correlations between the displacement of an ion $i$ of species $s$ and that of an ion $j$ of another species $s^{\prime}$:

[TABLE]

If the trajectories of the ions are independent from each other, $H_{s}=H_{ss^{\prime}}=0$ and the Nernst-Einstein approximation is recovered (i.e. $H=1$). However, $H=1$ doesn’t imply that $H_{s}$ and $H_{ss^{\prime}}$ $\sim 0$. For instance in \ceLi2CO3 (see Table 8), $H$ is very close to 1 ($H=1.04$) although $H_{\rm LiLi}=0.04$, $H_{\rm CO_{3}}=-0.32$ and $H_{\rm LiCO_{3}}=0.32$. These results mean that the displacements of Li cations are mostly not correlated to each other, whereas those of \ceCO3 are negatively correlated to each other, i.e. two \ceCO3 have a high probability to move in opposite directions. At variance, the cation-anion correlation term $H_{\rm Li_{2}CO_{3}}=0.32$ is positive, meaning that cations and anions are also moving in the opposite direction (notice that in equation (21) $z_{s}z_{s^{\prime}}<0$). The net result is that $H_{\rm CO_{3}}$ and $H_{\rm LiCO_{3}}$ are canceling each other, leading to $H=1+H_{Li}=1.04$. A different situation is observed with \ceNa2CO3 and \ceK2CO3, for which the cancellation effect implies the three contributions $H_{\rm X}$, $H_{\rm CO_{3}}$ and $H_{\rm XCO_{3}}$. As a result the anion-cation correlations $H_{\rm XCO_{3}}$ contribute positively (0.48 for \ceNa2CO3 and 0.46 for \ceK2CO3), whereas the correlations between ions of the same species decrease $H$ by a value of $H_{{\rm X}}+H_{{\rm CO}_{3}}$ ($-$0.49 for \ceNa2CO3 and $-$0.52 for \ceK2CO3), but these two contributions almost balance each other providing a Haven ratio $H=1+H_{\rm X}+H_{\rm CO_{3}}+H_{\rm XCO_{3}}$ equal to $\sim~{}1$ (0.99 for \ceNa2CO3 and 0.94 for \ceK2CO3).

This approach can be carried out on other compositions, including binary and ternary mixtures. This issue will be more thoroughly discussed in a future paper. The conclusion we want to emphasize here is that the Nernst-Einstein relation leads to a reasonable estimation of the electrical conductivity, not because the underlying assumption (the ions move independently from one another) is valid, but because of a cancellation effect between ion-ion correlations. Therefore the Nernst-Einstein relation must be carried with great caution, in particular when it is used to extract the ion diffusion coefficients.

V.3 Viscosity

Viscosity of molten alkali carbonates at ambient pressure has been measured since the sixties and critically reviewed by JanzJanz (1988). Its value stands between 3 and 30 mPa.s, depending on composition and temperature, as compared to 1 mPa.s, the typical viscosity of alkali chloride melts (NaCl at 1100 K) and liquid water (0.8 mPa.s at 300 K) at atmospheric pressure.

Afterwards, Sato et al. performed accurate measurements (errors within 3 %) of the viscosity of \ceLi2CO3, \ceNa2CO3 and \ceK2CO3 (as well as \ceRb2CO3 and \ceCs2CO3) Sato (1999), and of \ceLi2CO3–\ceNa2CO3 and \ceLi2CO3–\ceK2CO3 binary meltsSato et al. (1999), by using the oscillating method. However this technique may present some drawbacks (dependence on the apparatus, long time analysis…). On the other hand, the measurement of low viscosities at high temperature with the rotation method is tricky, but do not suffer from these disadvantages.Kim et al. (2015) Using this approach with a high-temperature rheometer system, Kim et al. (2015) obtained results that confirm those of Sato et al. for \ceLi2CO3, the \ceLi2CO3–\ceK2CO3 eutectic mixture of ratio 0.43:0.57 mol% (LKe2) and for the LNKe mixture, which has the lowest melting temperature (668 K).Sato (1999); Sato et al. (1999) However Di Genova et al. (2016) have recently applied a similar but improved method, as they claimed, and found viscosities systematically higher, by 30 to 100 %, compared to previous studies.

Concerning our simulations, they satisfactorily reproduce the experimental data of Sato (1999) for \ceNa2CO3 and \ceK2CO3 and in particular the activation energy (see Figure 9). As mentioned above, there is a large disagreement between the experimental study of Di Genova et al. (2016) and the one of Sato.Sato (1999) In particular the viscosity measured by Di Genova et al. for \ceLi2CO3 is twice the one measured by Sato, and leads to a surprisingly large difference between the viscosity of \ceLi2CO3 and those of \ceNa2CO3 and \ceK2CO3. Indeed, such a contrasting behavior between the small Li+ cation and the bigger Na+ and K+ cations, is not observed for other transport coefficients. As for our simulations of \ceLi2CO3, the calculated activation energy is compatible with the data but the magnitude of the viscosity is $\sim$ 15 % lower than that measured by Sato et al. (1999) and Kim et al.Kim et al. (2015) This result could be explained by the neglecting of some electronic effects (e.g. polarizability) in our model. Still, with regard to the viscosity of \ceNa2CO3 and \ceK2CO3 the calculated values are very close to the experimental ones (within the experimental uncertainties).

On the inset of Figure 9 the measured and calculated viscosity of the \ceLi2CO3–\ceNa2CO3–\ceK2CO3 eutectic mixture is also shown. A good agreement is found with both the reference data of Janz (1988) and Ejima (1987) and the recent measurements of An et al.An et al. (2016), emphasizing the ability of our force field to model the contribution of all the ions in mixtures. With regard to the LKe mixture, the viscosity estimated by our MD simulations is 10–30 % higher than the one measured by Sato et al.,Sato et al. (1999) with an activation energy equal to 42.6 kJ/mol (as compared to 34 kJ/mol), whereas the one calculated by Corradini *et al.*Corradini, Coudert, and Vuilleumier (2016) is by 20–30 % lower (with a lower activation energy of 29.6 kJ/mol ).

Like the other transport coefficients, the evolution of viscosity with temperature at high pressure is conveniently reproduced by an Arrhenius activation law given by equation (16), the parameters of which are given in Table 7. As shown in Figure S7 for low and moderate pressures ($P<$ 1 GPa), the viscosity depends almost exclusively on temperature. However beyond a few GPa, the weight of the activation volume $V_{a}^{\eta}$ is no longer negligible. Very recently, Stagno et al. (ress) reported viscosity measurements for \ceNa2CO3 in the 1.7–4.6 GPa pressure range, that are in agreement with our calculations (deviation within 0–50 %, Figure S7b).

Interestingly, the activation energy differs from that of electrical conductivity (Table 7). The viscosity is indeed dominated by the motion of slow ions, as opposed to the electrical conductivity which is controlled by the fastest species. Still the viscosity can be described fairly well by the phenomenological Stokes-Einstein equation:

[TABLE]

using 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, and $d$ is the hydrodynamic diameter of the diffusing ions. By choosing $d=1.6$, 2.2 and 2.4 Å  for \ceLi2CO3, \ceNa2CO3 and \ceK2CO3 respectively, we found that the behavior of the MD-calculated viscosity with pressure and temperature is correctly reproduced by equation (22) as demonstrated by Figure 10.

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

[TABLE]

where $A=$ 5.375, 2.615 and 1.825 for \ceLi2CO3, \ceNa2CO3 and \ceK2CO3 respectively (Figure S8). A comparable relation between these two transport properties has been also highlighted from experimental data of various melt compositions.Sifré, Hashim, and Gaillard (2015); Vuilleumier et al. (2014) The simple relations between $D$, $\sigma$ and $\eta$ can be very useful to infer one of these quantities from the others. They also show that these transport properties are essentially ruled by the diffusion of all ionic species in carbonate melts. This is in contrast with silicate melts whose viscosity is governed by the relaxation time of the Si-O bonds network, whereas the conductivity is mainly controlled by the mobility of network modifier elements, at least at low temperature near the melting point.

VI Conclusion

We have studied the thermodynamics (equation of state and surface tension), the microscopic structure and the transport properties (diffusion coefficients, electrical conductivity and viscosity) of molten alkali carbonates (\ceLi2CO3, \ceNa2CO3, \ceK2CO3 and some of their binary and ternary mixtures) in the $T-P$ range: $\sim$ 1100–2100 K, 0–15 GPa. For this we have developed an empirical force field using benchmark data from experiments (density and compressibility of the end-members at 1 bar) and from AIMD simulations (microscopic structure of four liquids). The density and compressibility data for binary and ternary mixtures that were not used in the FF parameterization procedure are also very well reproduced. We also showed for the end-members that the density of the crystalline phase at room conditions could be obtained with a good accuracy. Furthermore, the global agreement between the MD results and the full set of experimental data (thermodynamics ans transport coefficients) is very satisfying. This gives consistent evidence for the ability of the FF to be transferable in terms of composition and thermodynamic conditions. In this respect we believe it is a significant improvement compared to the previously existing force fields for molten alkali carbonates.

Using this force field, the equation of state of the end-members was first evaluated and modeled by a third-order Birch-Murnaghan formula. Next we have shown that mixtures behave ideally regarding the density and the compressibility. The analysis of the PDFs showed that the microscopic structure is that of an ionic liquid and the only modification we could identify upon increasing pressure was a compaction effect. As for the transport properties (diffusion coefficients, electrical conductivity and viscosity), they evolve smoothly (Arrhenius-like) over the studied $T-P$ domain.

Moreover, for the first time the surface tension of molten carbonates was calculated using MD. Although the calculated values are systemically higher than the experimental values published in the literature, the discrepancy is moderate ($\sim+15$ %), as observed in previous numerical studies dealing with ionic liquids, and could be due to computational uncertainties (size effects).

Finally we discussed the significance of the Nernst-Einstein and the Stokes-Einstein equations, that relate the diffusion coefficients to the electrical conductivity and to the viscosity, respectively. We showed that, by using ad hoc parameters, both formula lead to reasonable values. However they are not always representative of the transport mechanism itself. In particular we have demonstrated that the fairly good predictions provided by the Nernst-Einstein equation result from a partial canceling of interionic dynamic correlations.

Supplementary Material

See supplementary material for comparisons of the PDFs issued by MD and by AIMD simulations, a summary of all calculated properties with their uncertainties, plots of the evolution of these properties with $T$ and $P$, convergence plots of the correlation functions and details on the AIMD simulations including an example input file.

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).