A molecular simulation study of Cs-Cl and Cs-F ion pairs in hydrothermal fluids

Xi Zhang·Xiandong Liu·Tianhua Wang·Xiancai Lu·Rucheng Wang

Abstract Magmatic-hydrothermal processes play an important role in the transport,enrichment,and mineralization of cesium.In this study,classical molecular dynamics simulations were performed to investigate the properties of Cs-Cl and Cs-F ion pairs in hydrothermal fluids.The association constants(log10K A(m))under a wide range of temperature(i.e.298-1273 K)and fluid density(i.e.0.1-1.0 g/cm3)were derived from the potential of mean force(PMF)curves.The results indicate that Cs-Cl and Cs-F ion pairs have similar stabilities.This is different from other alkali metal cations(e.g.,Li+,Na+,and K+),which prefer binding with F over Cl.The stabilities of Cs-Cl and Cs-F ion pairs increase with increasing temperature(except for the fluid density≤0.1 g/cm3)or decreasing fluid density,which is similar to other alkali halide ion pairs.Comparisons among the stabilities of Cs-Cl/F and other alkali halide ion pairs indicate that the Li-F ion pair has the highest stability in hydrothermal fluids.

Keywords Cesium·Hydrothermal fluids·Association constant·Hydration structure·Molecular dynamics

1 Introduction

In recent years,cesium has been classified as one of the critical metals(U.S.Department of the Interior and U.S.Geological Survey 2017;European Union 2018).Cesium can be widely used in many modern technologies such as atom clocks,nuclear energy,and medical treatment(Takamoto et al.2005;Varala and Rao 2015;Kumar et al.2016).One major cesium source is granite pegmatites(U.S.Geological Survey,2011;Bradley et al.2017),in which cesium is usually enriched in ore minerals such as lepidolite((K,Cs)(Li,Al)3[(Al,Si)4O10](F,OH)2)and pollucite(CsAlSi2O6)(e.g.Cˇerny´2005;Wang et al.2004,2006).The magmatic-hydrothermal processes play an important role in cesium’s mineralization(e.g.London et al.1998;Wang et al.2007,2009),and the speciation of cesium is fundamental for understanding these processes.It has been widely accepted that the lithium-cesium-tantalum(LCT)type of granitic pegmatite,which is related to most kinds of cesium-bearing minerals,is enriched in some volatile components such as F and Cl(Cˇerny´and Ercit 2005;London 2018).Therefore,quantifying the stabilities of Cs-Cl and Cs-F ion pairs in hydrothermal fluids is useful for understanding the transport and enrichment of cesium during the magmatic-hydrothermal processes.

The properties of alkali halide ion pairs have been studied since the 1880s(Marcus and Hefter 2006).There are numerous experimental studies obtaining the association constants of these ion pairs in aqueous solutions.The most commonly used method is conductivity measurement(e.g.Franck 1961;Paterson et al.1971;Pethybridge and Spiers 1974,1977;Fuoss 1980;Hefter and Salomon 1996).For the stabilities of Cs-Cl/F ion pairs at the ambient conditions,Paterson et al.(1971)and Pethybridge and Spiers(1974,1977)reported their association constants in aqueous solutions with different concentrations.Fuoss(1980)further improved the electrolytic conductance theory and obtained more accurate results(KA),i.e.0.62 and 0.49 for Cs-Cl and Cs-F ion pairs,respectively.In the past two decades,dielectric spectroscopy was used to obtain association constants(e.g.Chen et al.2003;Brandes et al.2014).TheKAvalues of CsCl corresponding to the CIP(i.e.,contact ion pair)and SIP(i.e.,solvent shared ion pair)states derived in this way are 59±26 and 10.2±3.5 respectively(Brandes et al.2014),which are much higher than the conductimetric results.Furthermore,the study under elevated T-P conditions is quite rare.To the best of our knowledge,the only one is the conductivity measurements by Franck(1961),which obtained theKAvalues of Cs-Cl ion pair at different temperatures(i.e.823,923,and 1023 K)and fluid densities(i.e.0.3,0.5,and 0.7 g/cm3).

Molecular simulation has become a powerful tool for understanding the physicochemical properties of aqueous solutions at the micro-level(Allen and Tildesley 2017;Stubbs 2016;Brugger et al.2016).Different kinds of simulation methods have been used to calculate the association constants of ion pairs under different T-P conditions,e.g.classical molecular dynamics(CMD)(Chen and Pappu 2007;Fennell et al.2009;Mohoricˇand Bren 2017),first principles molecular dynamics(FPMD)(Wang et al.2021),and Monte Carlo(MC)method(Gujt et al.2014).Specifically,the association constants of Li/Na/K-F/Cl ion pairs in hydrothermal fluids have been obtained by using the CMD simulations(Zhang and Duan 2004;Yui et al.2010;He et al.2017;Zhang et al.2020;Wang et al.2021).However, the molecular dynamics study on the hydrothermal Cs-Cl/F ion pairs is still scarce except for one study of CsCl at 300-600 K(Mohoricˇand Bren 2017).

In this study,we performed classical molecular dynamics(CMD)simulations to calculate the association constants of Cs-Cl and Cs-F ion pairs in a wide range of temperature(i.e.298-1273 K)and fluid density (i.e.0.1-1.0 g/cm3).The microstructures of these ion pairs were also characterized by simulation trajectories.These results can provide a physicochemical basis for understanding the transport and enrichment of cesium in hydrothermal fluids.

2 Methods

2.1 Model

The simulation cells were cubic boxes with three-dimensional periodic boundary conditions.An ion pair(CsCl/CsF)and 250 water molecules were placed in these boxes randomly.The details of simulation systems at different T-ρ conditions are tabulated in Table S1 in theSupportingInformation.The previous study has shown that the finitesize effects can be avoided by using the systems containing 250 water molecules(Zhang et al.2020).For the simulations at the ambient conditions and 573 K,water densities at the corresponding saturated vapor pressures(Psat)were used,which were calculated from the equation of state of water(Wagner et al.2000).For higher temperatures,we firstly performed the simulations of Cs-Cl ion pair at 823 K-0.3 g/cm3and 1023 K-0.3 g/cm3with different force fields(see Sect.2.2),then a series of T-ρconditions(i.e.T=673,873,1073,and 1273 K;ρ=0.1,0.3 and 0.6 g/cm3)were used.The selected T-ρconditions can reflect the property of geological fluids from low-density vapor-like,in which cesium can also be transported(Foustoukos and Seyfried Jr.2007),to moderate-pressure liquid-like,e.g.,the pressures for pure water at 673 K-0.6 g/cm3and 1273 K-0.6 g/cm3are approximately equal to 0.5 kbar and 5.0 kbar,respectively(Zhang and Duan 2005).Besides,the choice of these T-ρconditions contributes to the direct comparisons with the association constants of Cs-Cl/F ion pairs derived from the conductivity measurements(Franck 1961;Fuoss 1980)and the CMD results of other alkali halide ion pairs(Yui et al.2010;He et al.2017;Zhang et al.2020;Wang et al.2021).’’

2.2 Classical molecular dynamics simulation

All CMD simulations were performed by using the DL_POLY 4.09 package(Todorov et al.2006)with the force field parameters tabulated in Table S2 in theSupporting Information.The SPC/E water model(Berendsen et al.1987)was employed.To evaluate the performances of different force fields,we first calculated the association constants of Cs-Cl/F ion pairs at the experimental conditions in Fuoss(1980)(i.e.,at the ambient conditions)and Franck(1961)(i.e.,at 823 K-0.3 g/cm3and 1023 K-0.3 g/cm3)with the parameters from Dang and Smith(1992,1995),Smith and Dang(1994)(denoted as‘‘DS’’)and Joung and Cheatham III(2008)(denoted as‘‘JC’’).The comparisons between the experimental and computational results suggested that the association constants calculated with the JC force field matched the experimental data better(see Sect.3.1),and thus the JC force field was used to calculate the association constants under other T-ρconditions tabulated in Table S1 in theSupporting Information.

Table 1 Comparisons among the log10K A(m)values of Cs-Cl ion pairs in previous experimental studies and this study

The cutoff radius for Columbic and van der Waals interactions was set to be 10.0 A˚,except for the simulations at the ambient conditions,in which a cutoff of 9.5 A˚ was used because of the side length of the simulation box(i.e.19.58 A˚).The previous study has shown that such a small difference in the cutoff radius does not influence the calculated PMF curves or association constants(Zhang et al.2020).Each simulation was performed for 1.0 ns with a time step of 0.5 fs,after an equilibration simulation of 100.0 ps.All simulations were performed in the NVT ensemble and the Nose´-Hoover chain thermostat was employed to control the temperature(Nose´1984a,b;Hoover 1985).

2.3 Potential of mean force(PMF)and association constant(KA)

The potential of mean force,W(r),between the cation and anion as a function of their distancerwas calculated by integrating the mean forcef(r)along withr:

In this study,we used the constrained MD method to derive the force of constraint,which is used to constrain the cation and the anion at a required distance.In each constrained MD simulation,the distance between the cation and anion was fixed,and the constraint force in every configuration was derived with the constraint algorithm(Ryckaert et al.1977;Andersen 1983).The mean forcef(r)can be treated as the average value of the constraint force during the production run.The distance between the two ions changed from 8.0 to 1.6 A˚with an interval of 0.2.r0is the reference distance which was set to 8.0.

The effective potential at the reference distance,W(r0),was calculated from Eq.(2)(Yui et al.2010):

whereεLJandσLJare the Lennard-Jones(LJ)parameters between Cs+and Cl-/F-tabulated in Table S2 in theSupporting Information,εis the relative dielectric constant of the SPC/E water model tabulated in Table S3 in theSupporting Information(Plugatyr and Svishchev 2009),andε0is the vacuum permittivity.

The association constants(KA)were calculated from the PMF(Justice and Justice 1976;Chialvo and Simonson 2003):

Table2Log10K A(m)valuesofCs-Cl/FionpairscalculatedwiththeJCforcefield 1273K 1073K 873K 673K 573K 298K T-ρcondition 0.6g/cm3 0.3g/cm3 0.1g/cm3 0.6g/cm3 0.3g/cm3 0.1g/cm3 0.6g/cm3 0.3g/cm3 0.1g/cm3 0.6g/cm3 0.3g/cm3 0.1g/cm3 0.712g/cm3 0.997g/cm3 1.45±0.02 1.57±0.03 3.67±0.03 4.14±0.04 7.93±0.06 8.79±0.07 1.35±0.03 1.49±0.03 3.64±0.05 4.10±0.05 8.54±0.07 9.36±0.11 1.22±0.03 1.31±0.04 3.43±0.05 3.80±0.06 8.76±0.08 9.48±0.08 1.02±0.03 1.06±0.05 2.83±0.05 3.09±0.04 7.98±0.12 8.40±0.08 0.58±0.04 0.56±0.04 0.02±0.11 0.16-0.02±Cs-Cl Cs-F

whereNAis the Avogadro constant,kBis the Boltzmann constant,Tstands for the temperature in Kelvin,andrcis the cutoff distance dividing the‘‘free ions’’and‘‘associated ion pair’’.In previous studies,both therccorresponding to the CIP state and CIP+SIP state have been taken to calculate the association constants from the ambient to the elevated conditions(e.g.Gao 1994;Cui and Harris 1994;Chialvo and Simonson 2003;Zhang and Duan 2004),and it was found that at the elevated conditions,the differences between the association constants calculated from the CIP state and CIP+SIP state become smaller(Yui et al.2010;He et al.2017).In this study,we firstly used therccorresponding to the CIP state and CIP+SIP state to calculate the association constants at the ambient conditions respectively and found that the results for the CIP+SIP state matched the experimental values better (see Sect.3.1).Therefore,at the elevated conditions,thercfor the CIP+SIP state was adopted.The calculatedKAwas converted into the association constant based on the molal concentration,i.e.KA(m),by multiplying the density of water,i.e.ρ(g/cm3):

3 Results

3.1 Association constants

The association constants for Cs-Cl ion pairs calculated from the CMD simulations with the DS force field and the JC force field are tabulated in Table 1.The values derived from conductivity measurements(Franck 1961;Fuoss 1980)are also included for comparison.At the ambient conditions,the calculated log10KA(m)values for the CIP+SIP state are 0.02 for both force fields,which are closer to the values based on the conductance data[i.e.-0.21 in Fuoss(1980)and-0.18 in Hefter and Salomon(1996),respectively]than the calculated values for the CIP state(i.e.-0.58 and-0.61 calculated with the DS and JC force field,respectively).At higher temperatures,the association constants calculated with the JC force field are about 0.2 logarithm units larger than those calculated with the DS force field and match the experimental data in Franck(1961)better(i.e.within 1.3 logarithm units).

For the Cs-F ion pair at the ambient conditions,the calculated log10KA(m)values for the CIP+SIP state calculated with the DS and JC force field are 0.00 and-0.02 respectively,and those for the CIP state are-1.10 and-0.93 respectively.Again,the results corresponding to the CIP+SIP state are in better agreement with the experimental values of-0.31(Fuoss 1980),-0.52(Hefter and Salomon 1996),and-0.37(Manohar and Atkinson 1993).

In general,there is no significant difference between the results calculated with these two force fields.Considering the slightly better performance of the JC force field at higher temperatures,the association constants at other T-ρ conditions were calculated with the JC force field.Fig.1 shows the calculated log10KA(m)values of Cs-Cl/F ion pairs as a function of T-ρconditions,and the data are tabulated in Table 2.The log10KA(m)values at the fluid density of 0.6 g/cm3are less than 1.6,indicating that the associations of these ion pairs are weak at higher densities.When the fluid densities are higher than 0.1 g/cm3,the log10KA(m)values of both Cs-Cl and Cs-F ion pairs increased with the temperature increasing or fluid density decreasing.At the density of 0.1 g/cm3,the curves of log10KA(m)values reach a maximum at the temperature of 873 K.These data indicate that the Cs-Cl/F ion pairs can be stable at high-temperature low-density fluids.The PMF curves in Fig.S1 in theSupporting Informationshow that at high temperatures and low densities,the CIP positions have the deepest wells,which also indicates that the ion pairs can be stable at these conditions.

Fig.1 Log10K A(m)values of a Cs-Cl and b Cs-F ion pairs(calculated with the JC force field)as a function of temperature and fluid density

3.2 Hydration structures

Figures 2 and 3 show the RDF and CN curves derived from the MD trajectories of hydrothermal Cs-Cl/F ion pairs at the CIP positions under different T-ρconditions.In Fig.2,the first peaks of Cs-O RDFs occur at about 2.9-3.0 A˚at all conditions.The corresponding CN values are tabulated in Table 3 and the snapshots are shown in Fig.4 for example.The results indicate that the hydration numbers of Cs+decrease with temperature and increase with fluid density.The hydration numbers of Cs+in the Cs-Cl ion pair are smaller than those in the Cs-F ion pair at the same condition,e.g.at 873 K-0.3 g/cm3,the numbers of coordinated water molecules around Cs+in Cs-Cl and Cs-F ion pairs are 4.8 and 6.5,respectively.This phenomenon suggests that F-can contribute to more coordinated water molecules around Cs+.The RDF and CN curves for Cl-H and F-H(at the CIP positions)are displayed in Fig.3,in which the first peaks for Cl-H RDF locate at about 2.3 A˚,and those for F-H RDF occur at about 1.7 A˚.Data in Table 3 show that the CNs for Cl-H are slightly smaller than those for F-H at the same condition,e.g.at 873 K-0.3 g/cm3,the CN for Cl-H is 3.1 while that for F-H is 3.6.Besides,the CNs for Cl-H and F-H also become smaller at higher temperatures or lower densities.

Fig.2 Cs-O RDF and CN curves in a,b Cs-Cl and c,d Cs-F ion pairs at the CIP conditions at a,cρ=0.30 g/cm3 and b,d T=873 K

Fig.3 Cl/F-H RDF and CN curves in a,b Cs-Cl and c,d Cs-F ion pairs at the CIP conditions at a,cρ=0.30 g/cm3 and b,d T=873 K

Fig.4 Snapshots of the hydration structures of Cs-Cl(a,c and e)and Cs-F(b,d,and f)ion pairs at the CIP positions.For clarity,only ions and neighboring water molecules are displayed.Aquamarine ball=cesium,chartreuse ball=chlorine,cyan ball=fluorine,red ball=oxygen,pink ball=hydrogen

Table 3 Coordination numbers of the first coordination shell at the CIP positions

4 Discussion

As shown in Fig.1 and Table 2,Cs-Cl and Cs-F ion pairs have similar association constants from the ambient to the elevated conditions.This is different from other alkali halide ion pairs,i.e.,Li/Na/K-Fvs.Li/Na/K-Cl.As shown in previous studies,the association constants of Li/Na/K-F ion pairs are higher than those of Li/Na/K-Cl ion pairs in hydrothermal fluids,especially in high-temperature lowdensity fluids(Wang et al.2021;Yui et al.2010;He et al.2017;Zhang et al.2020),implying that these alkali metal cations prefer forming ion pairs with F over Cl.This difference indicates that the transport of cesium in hydrothermal fluids can be different from that of other alkalis,i.e.,Cs-Cl and Cs-F ion pairs may have similar possibilities.

Figure 5 presents the comparison among the association constants of Cs-F and other alkali fluoride ion pairs(i.e.,Li/Na/K-F)in hydrothermal fluids(Wang et al.2021;Zhang et al.2020).It should be noted that the association constants of Cs-F ion pairs were calculated with the JC force field,whereas those for other alkali fluoride ion pairs were calculated with the DS force field.Compared with the DS force field,the JC force field produced slightly larger(but closer to the experimental data)results for Cs-F ion pairs(see Table 1)and lower results for Li/Na/K-F ion pairs(Wang et al.2021;Zhang et al.2020).It can be seen in Fig.5 that the association constants of hydrothermal alkali halide ion pairs show similar trends as temperature or fluid density changes.Among these alkali halide ion pairs,Li-F ion pairs have the highest association constants at all T-ρconditions,while the stabilities of Na/K/Cs-F and Cs-Cl ion pairs are similar.

Fig.5 Comparison among the calculated log10K A(m)values of Cs-F and Li/Na/K-F ion pairs in hydrothermal fluids.Data for Li/Na/K-F ion pairs were derived from Wang et al.(2021)and Zhang et al.(2020)

The hydration structures of Cs-Cl/F ion pairs(at the CIP positions)shown in Figs.2 and 3 and Table 3 suggest that their hydration numbers increase with temperature decreasing or fluid density increasing.Similarly,the hydration numbers of other alkali halide ion pairs(at the CIP positions)in previous CMD studies also showed such a temperature and fluid density dependence(Wang et al.2021;Zhang et al.2020).These phenomena indicate the weakening of both Cs-OH2and hydrogen bonding interactions in high temperature or low-density fluids,and thus the ion pairs are easier to associate under such conditions.Besides,Table S4 in theSupporting Informationpresents the comparison among the hydration numbers of Cs+and other alkali metal cations(i.e.,Li+,Na+,and K+)at the CIP positions at 873 K.It shows that under the same condition cesium has the highest hydration number while lithium has the lowest,e.g.,the numbers of coordinated water molecules of Li/Na/K/Cs-F ion pairs(at the CIP positions)are 2.9,4.0,5.8 and 6.5 at 873 K and 0.3 g/cm3,respectively.

5 Conclusions

In this study,we investigated the stabilities and microstructures of Cs-Cl and Cs-F ion pairs in aqueous solutions from the ambient to the elevated conditions by using classical molecular dynamics simulations.The calculated association constants show that Cs-Cl and Cs-F ion pairs have close stabilities,suggesting that Cl and F play a similar role in the transport of cesium in hydrothermal fluids.This is different from other alkali halide ion pairs(i.e.Li/Na/K-F/Cl)in which alkali ions prefer binding with fluorine over chlorine.Similar to other alkali halide ion pairs,the association constants of Cs-Cl/F ion pairs increase with temperature increasing or the fluid density decreasing for the density＞0.1 g/cm3.Among these alkali halide ion pairs,the Li-F ion pair has the highest stability in hydrothermal fluids.Structure analyses show that the hydration number of cesium is the highest while lithium has the lowest,and the hydration numbers of all alkali halide ion pairs generally decreased with temperature increasing or fluid density decreasing.

Supplementary InformationThe online version contains supplementary material available at https://doi.org/10.1007/s11631-022-00527-0.

AcknowledgementsThis study was supported by the National Natural Science Foundation of China(Nos 92062213,91855209,42125202 and 41872041).We acknowledge the financial support from the State Key Laboratory for Mineral Deposits Research at Nanjing University.We are grateful to the High Performance Computing Center(HPCC)of Nanjing University for doing the numerical calculations in this paper on its blade cluster system.

Declarations

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