Yttrium speciation in subduction-zone fluids from ab initio molecular dynamics simulations

Abstract

The rare Earth elements (REEs) are important geochemical tracers for geological processes such as high-grade metamorphism. Aqueous fluids are considered important carriers for the REEs in a variety of geological environments including settings associated with subduction zones. The capacity of a fluid to mobilize REEs strongly depends on its chemical composition and on the presence of suitable ligands such as fluoride and chloride. In this study, we present structural and thermodynamic properties of aqueous yttrium–chloride and yttrium–fluoride species at a temperature of 800 ${}^{\circ }$C in a pressure range between 1.3 and 4.5 $\mathrm{GPa}$ derived from ab initio molecular dynamics simulations. The total yttrium coordination by ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ and halide ions changes from seven to eight within the pressure range. For the yttrium–chloride species, a maximum number of three chloride ligands was observed. The derived thermodynamic data show that aqueous yttrium–fluoride complexes are more stable than their yttrium–chloride counterparts in chloride- and fluoride-rich environments at conditions relevant to slab dehydration. Mixed $\mathrm{Y}(\mathrm{Cl},\mathrm{F})$ complexes are found to be unstable even on the molecular dynamics timescale. Furthermore, in contrast to field observations, thermodynamic modeling indicates that yttrium should be mobilized at rather low fluoride concentrations in high-grade metasomatic systems. These results suggest a rather low fluoride activity in the majority of subduction-zone fluids because yttrium is one of the least-mobile REEs. Additionally, the simulations indicate that yttrium drives the self-ionization of hydration water molecules as it was observed for other high-field-strength elements. This might be a general property for highly charged cations in aqueous solutions under high-temperature and high-pressure conditions.

1 Introduction

Subduction zones have been the most important sites for exchange of matter and energy between the Earth's crust and mantle for billions of years until now . Magnetotelluric anomalies (e.g., ) suggest the occurrence of a high proportion of melts and water-rich fluids in the subducted slabs due to partial melting and the dehydration of water-bearing minerals such as serpentine and amphibole . Quite naturally, these fluids do not consist of pure water; much rather, they are brines with high salinity up to several mass percent of Cl and contain Si, Al and alkali cations (Na, K) as major solutes with minor amounts of Ca, Fe and Mg .

Aqueous fluids play an important role in subduction zones. Their interaction with minerals and rocks results in alterations including the dissolution and precipitation of minerals and/or the exchange of chemical elements and isotopes. Furthermore, fluid-mediated transport of trace elements such as high-field-strength elements (HFSEs) and rare Earth elements (REEs) is a major process of the deep element cycles . The distribution of these elements between minerals and fluids or melts is used as petrogenetic indicator for fractionation processes in igneous and metasomatic petrology . It is known that REE patterns of subducted rocks are affected by the chemical composition of the metamorphic fluid (e.g., ) due to their chemical complexation with dissolved anions, e.g., ${\mathrm{F}}^{-}$, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ and ${\mathrm{Cl}}^{-}$ .

The speciation of REEs at moderate pressures (up to few 100 $\mathrm{MPa}$) and high temperatures (250–300 ${}^{\circ }$C) was studied experimentally using in situ X-ray absorption spectroscopy (XAS) and solubility experiments (see a recent review by ) to understand physicochemical properties of hydrothermal fluids related to REE ore deposition focusing on chloride and fluoride complexes. Due to the high stability of fluoride complexes , it is a widely shared notion that fluoride complexes are most important for REE transport in hydrothermal fluids, but suggested that REE fluoride complexes are not the major carrier of REEs due to the low solubility of REE fluoride minerals such as bastnaesite, $(\mathrm{Ce},\mathrm{La},\mathrm{Nd},\mathrm{Y})\left[{\mathrm{FCO}}_{\mathrm{3}}\right]$, and due to the low fluoride activity in low pH environments. According to , REE chloride and sulfate complexes appear to be the main species for REE transport in hydrothermal systems. However, this interpretation is questioned by other authors .

So far, the number of in situ studies that address the complexation or thermodynamic properties of REE aqueous species at pressure ($P$) and temperature ($T$) conditions of subduction zones is very limited due to the challenging experimental setups . Apart from field observations (e.g., fluid inclusion analysis), our main understanding of the behavior of REEs under high $P/T$ conditions is derived from fluid/mineral partitioning and solubility experiments (e.g., ) and from numerical simulations. In a case study, modeled the hydration shell of REEs in solution by static energy calculations of an explicit first hydration shell and an implicit solvent model. Temperature effects were introduced by considering changes of the dielectric constant of the solvent and by calculations of the entropy. In this study, it is concluded that the hydration energies of all lanthanides become more similar at high $P/T$ conditions and that the availability of ligands becomes a controlling factor for the fractionation of light REEs (LREEs) and heavy REEs (HREEs) by subduction-zone fluids. Experiments suggest that LREEs (e.g., La) are more mobile than HREEs in chloride-rich solutions . The presence of fluoride in the system enhances the mobility of HREEs and this leads to fractionation processes .

From a geochemical perspective, yttrium is considered a HREE

International Union of Pure and Applied Chemistry: Nomenclature of Inorganic Recommendations 2005.

, and as such it is very common to use yttrium as a representative of the whole group of HREEs because of their similar chemical properties. Further, the comparable behavior of Y and the majority of the HREEs in high-grade metasomatism processes supports this assumption. The hydration shell of ${\mathrm{Y}}^{\mathrm{3}+}$ in aqueous solutions and possible complexation of ${\mathrm{Y}}^{\mathrm{3}+}$ with chloride has been subject to a number of studies at room temperature (e.g., ). Molecular simulations in conjunction with advanced sampling methods indicate that in the absence of other ligands ${\mathrm{Y}}^{\mathrm{3}+}$ is coordinated by eight hydration water molecules , whereas at high pH $\left[\mathrm{Y}\right(\mathrm{OH}{)}_{\mathrm{3}}({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{3}}{]}_{\mathrm{aq}}$ complexes are formed . The reported Y–O distance of 2.38 Å agrees with extended X-ray absorption fine structure (EXAFS) and X-ray absorption near edge structure (XANES) measurements .

The complexation of ${\mathrm{Y}}^{\mathrm{3}+}$ with ${\mathrm{Cl}}^{-}$ was studied by at hydrothermal conditions up to 340 ${}^{\circ }$C using in situ EXAFS spectroscopy. The authors claim that yttrium is coordinated by eight to nine neighbors at hydrothermal conditions and does not associate with chloride but rather forms polyatomic yttrium species. In another study by , a strong association of Y with chloride up to ${\mathrm{YCl}}_{\mathrm{4}}^{-}$ at 500 ${}^{\circ }$C and a linear reduction of the total number of coordinating atoms from eight to four towards high temperatures are reported. The results of that study indicate that yttrium behaves like 3d transition metal ions rather than a HREE under hydrothermal conditions at high chloride activity. Solubility experiments up to 1 $\mathrm{GPa}$ and 800 ${}^{\circ }$C in a hydrothermal piston-cylinder apparatus performed by indicate that yttrium is transported as $[\mathrm{YClOH}{]}^{-}$ complexes in $\mathrm{NaCl}$ brines and that ${\mathrm{YF}}_{\mathrm{2}}^{+}$ is the major complex in a fluorine-rich environment.

While the stability and distribution of $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ species in aqueous solutions at ambient conditions have been subject to a number of studies (e.g., ), the knowledge of thermodynamic properties of yttrium species in hydrothermal fluids is limited to theoretical predictions based on regressions using the Helgeson–Kirkham–Flowers (HKF) model and one experimental study by . Stability constants of $\mathrm{Y}-\mathrm{Cl}$ and $\mathrm{Y}-\mathrm{F}$ complexes at subcrustal high $P/T$ conditions have been barely investigated.

The capacity of a fluid to mobilize a certain element or to dissolve a certain amount of a component in the fluid depends on the chemical potential of the formed aqueous complexes . For a better understanding of the mobility of Y (as a representative of the HREE group) in subduction-zone fluids at high $P/T$ conditions, knowledge of the relation between the concentration of the molecular species in aqueous fluids and their thermodynamic properties is required. Yttrium in general is one of the most immobile REEs in high-grade metasomatic environments. But the high mobility in certain locations not only associated with hydrothermal ore deposits but also in metamorphic or diagenetic context indicates that the dissolution or transport of Y is constrained to a very narrow range of fluid compositions. Therefore, yttrium could be a potential indicator for certain geological fluid compositions. In this study, we use ab initio molecular dynamics (AIMD) simulations to investigate the atomic-scale structure and probe the free energy of different $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ complexes in the $P/T$ range of subduction zones.

Table 1

Number of atoms in the different simulation cells together with the size of the simulation cell. A and B refer to the system density of 1072 (1.3 $\mathrm{GPa}$) and 1447 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ (4.5 $\mathrm{GPa}$)${}^{\mathrm{a}}$.

Cell	${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$	${\mathrm{Y}}^{\mathrm{3}+}$	${\mathrm{Cl}}^{-}$	${\mathrm{F}}^{-}$	${\mathrm{Na}}^{+}$	No. atoms	$a^{\mathrm{b}}$
A1	84	1	6	0	3	262	14.29
A2	84	1	5	1	3	262	14.25
A3	84	1	4	2	3	262	14.21
A4	84	1	3	3	3	262	14.16
B1	84	1	6	0	3	262	12.93
B2	84	1	3	3	3	262	12.82

${}^{\mathrm{a}}$ Volume estimated using the empirical equation of state from for 2 $\mathrm{molal}$ $\mathrm{NaCl}$ solution. ${}^{\mathrm{b}}$ Edge length of the simulation box (Å).

2 Methods

2.1 Ab initio molecular dynamics

The AIMD simulation approach is based on a quantum-mechanical description of the electronic structure within the density functional theory (DFT) . Here, we used AIMD simulations with the Car–Parrinello method to model the molecular structure of $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ complexes in aqueous solutions. We performed simulations with the widely used Car–Parrinello molecular dynamics (CPMD) code . Within the code, the BLYP exchange correlation functional was employed and the plane-wave expansion of the Kohn–Sham orbitals was truncated at a cutoff energy of 80 $\mathrm{Ry}$. To reduce the computation effort, the core electrons of all atoms in the cubic simulation cell were approximated by Goedecker-type pseudopotentials . To separate the electronic and nuclear motion of the CPMD, a fictive electron mass of 600 $\mathrm{a}.\mathrm{u}.$ with fictitious kinetic energy of 0.24 $\mathrm{a}.\mathrm{u}.$ was used. All results presented here are based on simulations performed with a constant number of atoms $N$ at constant volume $V$ and a temperature $T$ of 800 ${}^{\circ }$C (a so-called $NVT$ ensemble) and with a time step of 0.1 $\mathrm{fs}$. The temperature in the simulation was controlled using a Nosé thermostat . We chose two different pressure conditions of 1.3 and 4.5 $\mathrm{GPa}$ for this study. The volumes of the simulation cell were estimated from the correlation function provided by assuming 2 $\mathrm{molal}$ $\mathrm{NaCl}$ solution for all cells (see Table ). Due to this approximation, the pressures listed in Table  have to be considered estimates.

The initial atomic configurations were derived from AIMD simulations of pure $\mathrm{NaCl}$ solutions (200 ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ and 10 $\mathrm{NaCl}$). This configuration had been equilibrated for a few tens of picoseconds ($\mathrm{ps}$) at 1000 $\mathrm{K}$ and a simulation box size of 13.74 Å. The original $\mathrm{NaCl}$ solutions were generated in classical molecular dynamics simulations using Matsuoka–Clementi–Yoshimine (MCY) potentials . We substituted one of the sodium atoms by yttrium and decreased the number of water molecules and chlorine atoms stepwise until we reached configuration A1 in Table . Figure  shows a snapshot of the simulation cell A1 with a $\left[{\mathrm{YCl}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$ complex. The chloride ions not initially bonded to the yttrium ion are constrained to remain at larger distance (6–7 Å) from the yttrium ion. All other simulation boxes (see Table ) were generated from this initial one and equilibrated for several picoseconds.

For the fluoride-bearing cells, we used different cell compositions due to the strong association of hydrogen and fluoride at low pressures. Only initially bonded ${\mathrm{F}}^{-}$ ions were included to avoid the formation of hydrofluoric acid by the reaction ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}+{\mathrm{F}}^{-}\rightleftharpoons {\mathrm{OH}}^{-}+\mathrm{HF}$.

Figure 1

Snapshot of simulation cell A1 with a $\left[{\mathrm{YCl}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$ complex. The water molecules are indicated by red–white bond sticks, sodium by yellow balls and chlorine by cyan balls. The ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ in the first hydration shell of the yttrium atom (copper colored) is presented as red–white balls and sticks. The constraint distances between the yttrium ion and the constraint ${\mathrm{Cl}}^{-}$ are colored in gray.

[Figure omitted. See PDF]

2.2 Analysis of interaction distances and coordination number

The average atomic structure of disordered systems such as aqueous solutions is commonly described in terms of partial radial distribution functions $g_{ij}\left(r\right)$. These functions describe the probability for finding a pair of atoms of elements $j$ and $i$ at a distance $r$ normalized to the particle number density ${\mathit{\rho }}_N$ of the fluid:

1$$g_{ij}\left(r\right)={\displaystyle \frac{\mathrm{1}}{c_{\mathrm{i}}{c}_j\mathrm{4}\mathit{\pi }{r}^{\mathrm{2}}{\mathit{\rho }}_{N}N}}\sum _{a=\mathrm{1}}^{N_{\mathrm{i}}}\sum _{b=\mathrm{1}}^{N_j}\mathit{\delta }(r-|{\mathit{R}}_a-{\mathit{R}}_b\left|\right),$$ where $c_{\mathrm{i}}(=N_{\mathrm{i}}/N)$ and $c_j$ are the concentrations of elements $i$ and $j$, $N$ is the number of particles in the simulation box, $\mathit{\delta }\left(x\right)$ the Dirac delta function, and ${\mathit{R}}_a$ and ${\mathit{R}}_b$ are the position vectors of particles $a$ and $b$. In the numerical implementation of Eq. (), the distance of each particle in respect to all other particles is calculated at every time step. The evolving list of distances is normalized to the number of particles and the volume of the simulation cell.

The positions of the first maximum of $g_{ij}\left(r\right)$ represent the distances with the highest density of particles around the central position, which are usually interpreted as the nearest-neighbor distances (or bond distances) between elements $i$ and $j$. The average coordination number of one element is derived from counting the number of neighbors for each atom of this element within a given cutoff distance respecting periodic boundary conditions and averaging over all particles of the same kind and over time. The cutoff distance is taken from the first minimum of the respective $g_{ij}\left(r\right)$. The association of ${\mathrm{OH}}^{-}$ groups with cations is evaluated by considering oxygen atoms coordinated by one hydrogen only. To distinguish pure ${\mathrm{OH}}^{-}$ from two ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ sharing one hydrogen, the cutoff between the oxygen and the hydrogen is set to 1.3 Å. Additionally, the distance of the hydrogen within this cutoff distance to the next oxygen is taken as the distinguishing criterion. Only if the oxygen of the next water molecule is located at a distance of $>\mathrm{1.6}$ Å, the ${\mathrm{OH}}^{-}$ is counted as hydroxide. This value represents approximately the hydrogen bond distance between ${\mathrm{OH}}^{-}$ and ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ . To evaluate the formation of a certain species during a simulation run, only complexes with a constant coordination of chloride and fluoride over at least 3 $\mathrm{ps}$ are considered. The average halogen ion hydration number is computed by counting the number of hydrogen oriented towards the ion of the vicinal water molecules.

Figure 2

(a) Potential of mean force of the dissociation reaction of ${\mathrm{YCl}}_{\mathrm{3}}$ to ${\mathrm{YCl}}_{\mathrm{2}}^{+}$ at 1.3 $\mathrm{GPa}$ and 800 ${}^{\circ }$C over a distance between 2.63 Å and 6.0 Å. The evolution of the Helmholtz free energy is shown in panel (b). In panel (c), an example of the progress of the constraint force with simulation time at a Y–Cl distance of 3.0 Å is shown (stage e). Panels (d)–(g) indicate the different stages of the dissociation of the initial complex (see text for details).

[Figure omitted. See PDF]

2.3 Constraint molecular dynamics simulations and thermodynamic integration

A single MD simulation only yields the total internal energy of the system. Thermodynamic integration (TI) is used to derive free energy differences between different states. This approach usually requires a number of intermediate MD simulations along a certain integration variable . Here, we use the constrained molecular dynamics approach and thermodynamic integration in terms of the blue moon sampling . This method has been used already by different groups to investigate the stability of metal complexes in aqueous solutions, not only at ambient pressure and temperature conditions but also at hydrothermal and deep crustal high-density fluid conditions .

Using this method, the Helmholtz free energy difference (${\mathrm{\Delta }}_{\mathrm{r}}A$) of a chemical dissociation reaction is obtained from the average forces $F\left(r\right)$ between two atoms under the constraint of keeping their bond distance $r$ constant. ${\mathrm{\Delta }}_{\mathrm{r}}A$ is derived from the integration of $F\left(r\right)$ between two different distances ($r_{\mathrm{1}}$ and $r_{\mathrm{2}}$), which correspond to the associated and the dissociated states of the aqueous complexes :

2$${\mathrm{\Delta }}_{\mathrm{r}}{A}_{\mathrm{1}⟶\mathrm{2}}=-\underset{\mathrm{1}}{\overset{\mathrm{2}}{\int }}⟨F\left(r\right)⟩\mathrm{d}r.$$ The formal relation between the Helmholtz free energy and the Gibbs free energy is given by 3$${\mathrm{\Delta }}_{\mathrm{r}}G={\mathrm{\Delta }}_{\mathrm{r}}{A}_{\mathrm{1}⟶\mathrm{2}}+V\underset{\mathrm{1}}{\overset{\mathrm{2}}{\int }}\mathrm{d}P.$$ $V$ is the volume of the simulation cell. Here, we use an $NVT$ ensemble, and the change in pressure averaged over the whole trajectory is approximately zero (${\int }_{\mathrm{1}}^{\mathrm{2}}\mathrm{d}P=\mathrm{0}$).

To look into the formation of $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ species in equilibrium reactions, we removed one of the ligands from the Y ion in multiple integration steps by constraining the $\mathrm{Y}-\mathrm{Cl}$ and $\mathrm{Y}-\mathrm{F}$ distances. During the integration, the $\mathrm{Y}-\mathrm{Cl}$ and $\mathrm{Y}-\mathrm{F}$ distances of the Cl or F ions that are not initially bonded to the yttrium ion are fixed at 6.0 to 7.0 Å to avoid disruptions during the integration. Figure  illustrates an example of the dissociation reaction in the simulation box at 1.3 $\mathrm{GPa}$ and 800 ${}^{\circ }$C: R1$$\left[{\mathrm{YCl}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\to \left[{\mathrm{YCl}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}+{\mathrm{Cl}}^{-}.$$

The integration starts at the first distance (Fig. d), which corresponds approximately to the equilibrium distance of the Y–Cl contact ion pair. This interatomic distance (Fig. a) is estimated as the intercept with the zero force line from a linear interpolation between the first and the second integration steps, the first step starting at 2.6 Å for Y–Cl (and at 2.0 Å for Y–F). With increasing displacement of the chloride ion, the constraint force is attractive (Fig. e) until a water molecule takes the place of the ion in the first coordination shell around the ${\mathrm{Y}}^{\mathrm{3}+}$ (Fig. f). At this point, the force becomes repulsive. By integration over the potential of mean force (PMF) (Fig. a), the Helmholtz free energy (${\mathrm{\Delta }}_{\mathrm{r}}A$) difference between the initial complex and the product of the reaction is derived (Fig. b). For the Y–Cl complexes, we assume a ligand in a distance of 6.0 Å  as being dissociated, and for Y–F this distance reduces to 5.0 Å (see Fig. S1 in the Supplement ). In Fig. g, the dissociation is completed and $\left[{\mathrm{YCl}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$ is formed. To estimate the convergence of the constraint force, the standard deviation of the average force is computed. As a convergence criterion, a value of 5 $\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$ over the last 2 $\mathrm{ps}$ is taken. This value is also considered as approximate error of the computed reaction free energies. To satisfy this criterion, the constraint AIMD simulations are performed for values between 4.5 and 40 $\mathrm{ps}$.

2.4 Thermodynamic approaches

It is textbook knowledge that the formation of monomeric complexes in equilibrium reactions develops in steps . During the formation process, a ligand $L$ is added to the metal cation $M$. This formation of a $M{L}_n$ complex can be written as a sequence of stepwise reactions: $$\begin{array}{c}{\displaystyle }M+L\to ML\\ {\displaystyle }ML+L\to M{L}_{\mathrm{2}}\\ {\displaystyle }\mathrm{⋮}\\ {\displaystyle }M{L}_{n-\mathrm{1}}+L\to M{L}_n,\end{array}$$ with the respective logarithmic equilibrium constants:

4$$\log K_n={\displaystyle \frac{-{\mathrm{\Delta }}_{\mathrm{r}}G{{}^{\circ }}_n}{\mathrm{2.303}RT}}.$$ Having determined equilibrium constants for all of those reactions, the stability constant (also referred to as cumulative stability constant or overall stability constant) ${\mathit{\beta }}_n$ of species $M{L}_n$ is defined as 5$$\log {\mathit{\beta }}_n=\log K_{\mathrm{1}}+\log K_{\mathrm{2}}+\log K_{\mathrm{3}}\mathrm{\dots }\log K_n.$$ The standard Gibbs free energy (${\mathrm{\Delta }}_{\mathrm{r}}G{{}^{\circ }}_{\mathrm{i}}$) depends on the reaction Gibbs free energy (${\mathrm{\Delta }}_{\mathrm{r}}{G}_{\mathrm{i}}$) derived from the MD simulation, temperature $T$, gas constant $R$, molality of the ions $m_{\mathrm{i}}$ and the activity coefficient ${\mathit{\gamma }}_{\mathrm{i}}$: 6$${\mathrm{\Delta }}_{\mathrm{r}}G{{}^{\circ }}_{\mathrm{i}}={\mathrm{\Delta }}_{\mathrm{r}}{G}_{\mathrm{i}}-RT\ln {\displaystyle \frac{m_{M{L}_{\mathrm{i}}}{\mathit{\gamma }}_{M{L}_{\mathrm{i}}}}{m_{M{L}_{\mathrm{i}-\mathrm{1}}}{\mathit{\gamma }}_{M{L}_{\mathrm{i}-\mathrm{1}}}\cdot m_L{\mathit{\gamma }}_L}}.$$ The standard state for a solute in aqueous solution is a 1 $\mathrm{molal}$ hypothetical solution with properties of an infinitely diluted solution . The concentration and behavior of the solutes in the simulation cell (see Table ) are quite different from this standard state. Therefore, we computed the activity coefficient corresponding to this hypothetical solution using the B-dot model , which is an empirical extension of the Debye–Hückel theory : 7$$\log {\mathit{\gamma }}_{\mathrm{i}}=-{\displaystyle \frac{z_{\mathrm{i}}^{\mathrm{2}}{A}_{\mathrm{DH}}\sqrt{I}}{\mathrm{1}+{\mathit{å}}_{\mathrm{i}}{B}_{\mathrm{DH}}\sqrt{I}}}+\dot{B}I,$$ where $z_{\mathrm{i}}$ is the charge of ion $i$, å the mean distance of closest approach between the ions and $I$ the ionic strength: 8$$I=\mathrm{0.5}\sum _{i=\mathrm{1}}^n{m}_{\mathrm{i}}{z}_{\mathrm{i}}^{\mathrm{2}}.$$ $A_{\mathrm{DH}}$ and $B_{\mathrm{DH}}$ are the Debye–Hückel parameters: $$\begin{array}{}\text{9}& {\displaystyle }{A}_{\mathrm{DH}}={\displaystyle \frac{\mathrm{1.8248}\cdot {\mathrm{10}}^{\mathrm{6}}\sqrt{{\mathit{\rho }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}}{(T\mathit{ϵ}{)}^{\frac{\mathrm{3}}{\mathrm{2}}}}}\text{10}& {\displaystyle }{B}_{\mathrm{DH}}={\displaystyle \frac{\mathrm{50.292}\sqrt{{\mathit{\rho }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}}{(T\mathit{ϵ}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}}}.\end{array}$$ They depend on temperature, density (${\mathit{\rho }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$) and dielectric constant ($\mathit{ϵ}$) of the solvent. $\mathit{ϵ}$ is computed using the equation provided by and for pure ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ assuming a fluid density of the simulation. The value of $\dot{B}$ is calculated by the CHNOSZ software package applying the extrapolation suggested by . The $\stackrel{\mathrm{˚}}{a}$ parameter for the different complexes and ions (for all $\mathrm{Y}-\mathrm{Cl}/\mathrm{F}$ complexes, a constant value of 4.5 Å is applied) are taken from and Eq. () was solved in a Python implementation of the EQBRM program .

Table 2

The listed atomic distances and coordination numbers are averaged over the whole-simulation runs for all unbiased simulations. “AIMD time” corresponds to the total simulation time, whereas lifetime (“LT”) refers to persistence of the initial yttrium–halide complex. Under “formed”, the most abundant complexes${}^{\mathrm{a}}$ of the last 10 ps of the simulation are listed.

		Distances (Å)			Coordination numbers						Complex
ID	Cell	Y–O${}^{\mathrm{b}}$	Y–Cl	Y–F	Y–O	$\mathrm{Y}-{\mathrm{OH}}^{-}$	Y–Cl	Y–F	$\mathrm{Y}(\mathrm{Cl},\mathrm{F})-\mathrm{Na}$	$\mathit{\varnothing }$ Y	Initial	Formed	AIMD time/LT (ps)
# 1	A1	2.32	2.58	–	5.9	0.6	1.0	–	0.1	6.9	$\left[\mathrm{YCl}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{\mathrm{2}+}$	$\left[\mathrm{YClOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	25/25
# 2	A1	2.36	2.59	–	4.8	0.3	2.0	–	0.5	6.8	$\left[{\mathrm{YCl}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}^{+}$	$\left[{\mathrm{YCl}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	23/23
# 3	A1	2.37	2.60	–	3.6	0.0	3.0	–	0.4	6.6	$\left[{\mathrm{YCl}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{3}}{]}_{\mathrm{aq}}$	$\left[{\mathrm{YCl}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$	24/24
# 4	A1	2.35	2.58	–	2.9	0.0	3.4	–	0.8	6.3	$\left[{\mathrm{YCl}}_{\mathrm{4}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{2}}{]}^{-}$	$\left[{\mathrm{YCl}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{3}}{]}_{\mathrm{aq}}\cdot {\mathrm{Na}}^{+}$	26/14
# 5	A1	2.41	2.63	–	1.2	0.0	4.8	–	1.2	6.0	$\left[{\mathrm{YCl}}_{\mathrm{5}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}){]}^{\mathrm{2}-}$	$\left[{\mathrm{YCl}}_{\mathrm{4}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{2}}{]}^{-}\cdot {\mathrm{Na}}^{+}$	26/22
# 6	A2	2.37	–	2.08	5.6	1.0	–	1.0	–	6.6	$\left[\mathrm{YF}\right(\mathrm{OH}\left)\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	$\left[\mathrm{YFOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	25/25
# 7	A3	2.39	–	2.08	4.8	0.0	–	2.0	–	6.8	$\left[{\mathrm{YF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$,	$\left[{\mathrm{YF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	29/29
# 8	A4	2.43	–	2.11	3.5	0.0	–	3.0	–	6.5	$\left[{\mathrm{YF}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$	$\left[{\mathrm{YF}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$	29/29
# 9	A2	2.39	2.63	2.08	4.8	0.1	1.0	1.0	0.2	6.7	$\left[\mathrm{YClF}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	$\left[\mathrm{YClF}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	29/29
# 10	A2	2.39	2.62	2.10	3.9	0.1	1.6	1.0	0.3	6.5	$\left[{\mathrm{YCl}}_{\mathrm{2}}\mathrm{F}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$	$\left[\mathrm{YClF}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	29/ 17
# 11	A3	2.40	2.62	2.07	4.3	0.0	0.4	2.0	0.0	6.8	$\left[{\mathrm{YClF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$	$\left[{\mathrm{YF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}\cdot {\mathrm{Na}}^{+}$	27/12
# 12	A1	2.35	–	–	7.2	0.5	–	–	–	7.2	$\left[\mathrm{Y}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{3}+}$	$\left[\mathrm{YOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{\mathrm{2}+}$	29/29
# 13	B1	2.36	2.64	–	7.4	0.4	0.5	–	0.0	7.8	$\left[\mathrm{YCl}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{\mathrm{2}+}$	$\left[\mathrm{YCl}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$	27/13
# 14	B1	2.34	2.65	–	5.9	0.1	2.0	–	0.3	7.9	$\left[{\mathrm{YCl}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	$\left[{\mathrm{YCl}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{+}$	27/27
# 15	B1	2.36	2.61	–	7.4	0.4	0.5	–	0.0	7.9	$\left[{\mathrm{YCl}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$	$\left[\mathrm{YOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$	27/0
# 16	B2	2.34	–	2.06	6.7	0.2	–	1.0	0.1	7.7	$\left[\mathrm{YF}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$	$\left[\mathrm{YF}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$	24/24
# 17	B2	2.38	–	2.08	5.5	0.2	–	2.0	0.4	7.5	$\left[{\mathrm{YF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	$\left[{\mathrm{YF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{+}$	27/27
# 18	B2	2.39	–	2.15	5.6	0.1	–	2.6	1.0	8.1	$\left[{\mathrm{YF}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$	$\left[{\mathrm{YF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{+}\cdot {\mathrm{Na}}^{+}$	25/14
# 19	B2	2.36	2.70	2.10	6.5	0.1	0.4	1.0	1.0	8.0	$\left[\mathrm{YClF}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$	$\left[\mathrm{YF}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}\cdot {\mathrm{Na}}^{+}$	25/10
# 20	B2	2.36	2.71	2.11	5.9	0.1	0.8	1.0	0.6	7.6	$\left[{\mathrm{YCl}}_{\mathrm{2}}\mathrm{F}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}_{\mathrm{aq}}$	$\left[\mathrm{YFCl}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{+}$	25/3
# 21	B2	2.34	2.78	2.13	5.2	0.0	0.4	2.0	0.4	7.6	$\left[{\mathrm{YClF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}_{\mathrm{aq}}$	$\left[{\mathrm{YF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{+}\cdot {\mathrm{Na}}^{+}$	27/9
# 22	B1	2.32	–	–	7.8	0.8	–	–	0.0	7.8	$\left[\mathrm{Y}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{8}}{]}^{\mathrm{3}+}$	$\left[\mathrm{YOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$	27/27

${}^{\mathrm{a}}$ A list of all formed complexes observed for at least 3 $\mathrm{ps}$ over the AIMD is given in the Supplement in Table S3. ${}^{\mathrm{b}}$ The decomposed $\mathrm{Y}-\mathrm{OH}{}_{\mathrm{2}}$ and $\mathrm{Y}-{\mathrm{OH}}^{-}$ distances are listed in Table S1.

Note that in the simulations we modeled a dissociation reaction but to be in line with the modern nomenclature of aqueous geochemistry all the derived thermodynamic data presented below correspond to the formation reaction. In this study, we investigated two kinds of reactions: R2$${\mathrm{Y}}^{\mathrm{3}+}+n{\mathrm{Cl}}^{-}\to {\mathrm{YCl}}_n^{\mathrm{3}-n}(n=\mathrm{1}-\mathrm{3})$$ and R3$${\mathrm{Y}}^{\mathrm{3}+}+n{\mathrm{F}}^{-}\to {\mathrm{YF}}_n^{\mathrm{3}-n}(n=\mathrm{1}-\mathrm{3}).$$

In the following, all results referring to one of those reactions are indexed by TI $-n$. Furthermore, for reasons of clarity, we do not include the number of hydration water molecules in the formula of the aqueous complexes in some of the presented figures and tables.

3 Results

3.1 Yttrium coordination in high-density aqueous fluid

AIMD simulations were performed for the hydrated ${\mathrm{Y}}^{\mathrm{3}+}$ and for 11 different yttrium–halogen complexes: five $\mathrm{YCl}$${}_n^{\mathrm{3}-n}$, $n=\mathrm{1}$–5; three $\mathrm{YF}$${}_n^{\mathrm{3}-n}$, $n=\mathrm{1}$, 2, 3; and three mixed $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ complexes. Simulation conditions and obtained structural data are compiled in Table . Moreover, the formed aqueous species are listed in Table . Note that the composition of the simulation cells varies slightly as cells A1–A4, B1 and B2 contain different amounts of F and Cl.

Figure 3

Radial distribution functions of $\mathrm{Y}-(\mathrm{Cl},\mathrm{O},\mathrm{F})$ scaled to the maximum of the $g_{ij}\left(r\right)$ from run # 9 and 12 together with a snapshot of a $\left[\mathrm{YClF}\right(\mathrm{OH}\left)\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}]$ complex. In the snapshot (from a run at 1.3 $\mathrm{GPa}$ and 800 ${}^{\circ }$C), the central yttrium atom is surrounded by chlorine (cyan), fluorine (green), a hydroxyl group and water molecules (red – O and white – H balls). Water molecules of the second shell are shown as blue sticks. This visualization illustrates the relation between $g_{ij}\left(r\right)$ and the atomic structure of the aqueous species. The colors of the ligands in the snapshot are equivalent to those in the $g_{ij}\left(r\right)$ functions.

[Figure omitted. See PDF]

Different partial radial distribution functions for Y–(${\mathrm{Cl}}^{-}$, ${\mathrm{F}}^{-}$, ${\mathrm{OH}}^{-}$ and ${\mathrm{OH}}_{\mathrm{2}}$) are shown in Fig. . To facilitate the comparison of the first peaks, the $g_{ij}\left(r\right)$ values are scaled to equal maximum intensity. The observed sequence of atomic distances between the central metal ion and its ligands holds for all the complexes. The shortest distance is found between ${\mathrm{Y}}^{\mathrm{3}+}$ and ${\mathrm{F}}^{-}$ ($\sim \mathrm{2.1}$ Å), followed by $\mathrm{Y}-{\mathrm{OH}}^{-}$ ($\sim \mathrm{2.1}$–2.2 Å) and Y–${\mathrm{OH}}_{\mathrm{2}}$ ($\sim \mathrm{2.3}$–2.4 Å). The largest distance is observed for Y–Cl pairs ($\sim \mathrm{2.6}$–2.8 Å) (for an overview, see Tables  and S1). Further, in $g_{\mathrm{Y}-\mathrm{O}}\left(r\right)$, a second maximum is observed. It corresponds to the second hydration shell, which is formed around all complexes at all studied $P/T$ conditions. For all fluid compositions, association of ${\mathrm{NaCl}}_n$ ($n=\mathrm{1}$–3) and of ${\mathrm{Na}}^{+}$ with the yttrium complex is observed. The mean Na–Cl coordination, NaCl species distribution and second hydration shell positions are listed in Table S1.

Figure 4

Average yttrium coordination by chloride, fluoride and oxygen for run # 1–22.

[Figure omitted. See PDF]

Figure 5

Presence of selected ion pairs or complexes during the different AIMD simulations. This includes the formation of ${\mathrm{OH}}^{-}$ within the first hydration shell of yttrium, the association of sodium with the coordinating halogens and the stability of the initial Y–halide complexes over AIMD time.

[Figure omitted. See PDF]

At 1.3 $\mathrm{GPa}$, the average Y coordination by O, Cl and/or F is about seven (see Fig. ) with two exceptions, run # 4 and 5, even if the initial coordination is lower (run # 1–3). In run # 1–3, the initial Y–Cl coordination is retained over the entire AIMD runtime, whereas for # 4 and 5, the 4-fold and 5-fold Cl coordinations do not persist and the time average yttrium coordination is below seven. The initial complex of # 5 seems to be unstable and transforms to the initial complex of # 4 after 26 ps. In # 4, after the 4-fold coordinated chloride complex is dissociated, a total coordination number of seven is reached at the end of the simulation run. Frequently, the formation of ${\mathrm{OH}}^{-}$ in the first hydration shell of yttrium is observed. The major hydroxide formation mechanism will be discussed below.

Figure  provides an overview over the formation and dissolution of selected structural units in the course of the simulations, i.e., the stability of the initial complexes, Y–hydroxide association and bonding of the coordinating halogen to sodium. All five chloride complexes associate with sodium. The strongest association is observed in # 2–5, where the sodium is connected to one or two chloride ligands of the Y complex. This association increases with the number of halide ligands in these complexes. Furthermore, in # 4 and 5, this association initiates the dissociation of $\left[{\mathrm{YCl}}_{\mathrm{4}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{2}}{]}^{-}$ and $\left[{\mathrm{YCl}}_{\mathrm{5}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}){]}^{\mathrm{2}-}$. Even larger clusters of sodium, constraint chlorides and the Y–chloride complex appear over timescales of less than 3 $\mathrm{ps}$. Moreover, from # 1–5, the average Y coordination decreases with the increasing number of initial chloride ligands (Fig. ). The Y–O distances do not change significantly with the increasing number of chloride ligands from # 1 to 4. Only in run # 5 is a significantly longer Y–O distance of 2.41 Å observed. The Y–chloride distances range from 2.58 Å in # 1 to 2.63 Å in # 5.

For pure Y–fluoride complexes under the same conditions (# 6–8), a slight increase of the Y–O distance with increasing number of fluoride ligands from 2.37 Å in $\left[\mathrm{YF}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{\mathrm{2}+}$ to 2.43 Å in $\left[{\mathrm{YF}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$ is observed. In all three runs, the initial complex persists over the whole simulation time. As in the case of Y–chloride solutions, the association of fluoride with sodium is observed but it is less pronounced. In # 6, where only one fluoride is initially bonded to the central ion, ${\mathrm{OH}}^{-}$ is formed within the first hydration shell of yttrium. The $\mathrm{Y}-\mathrm{F}$ distances within the complexes are approximately 0.5 Å shorter than those of the $\mathrm{Y}-\mathrm{Cl}$ species.

For the $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ mixed complexes (run # 9–11), only the run starting initially from $\left[\mathrm{YClF}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$ does not show the formation of multiple complexes over time. Here, only short separations of the ${\mathrm{Cl}}^{-}$ over $\sim \mathrm{1}$ $\mathrm{ps}$ from the complex occur. In # 10, $\left[{\mathrm{YCl}}_{\mathrm{2}}\mathrm{F}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$ dissociates to $\left[\mathrm{YClF}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{+}$ after 11.5 $\mathrm{ps}$. This complex is present over approximately 10 $\mathrm{ps}$, in conjunction with the formation of $\left[\mathrm{YClFOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$, followed by the reassociation of the initial complex. In # 11, starting from $\left[{\mathrm{YClF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{4}}{]}_{\mathrm{aq}}$, the initially bonded chloride is released after $\sim \mathrm{12}$ $\mathrm{ps}$ and $\left[{\mathrm{YF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$ is formed. The $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ distances of the mixed complexes are comparable to those of the pure ones.

In run # 12, starting from $\left[\mathrm{Y}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{3}+}$, hydroxide is formed within the first 8 $\mathrm{ps}$ (see Fig. ), which results in the formation of $\left[\mathrm{YOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{\mathrm{2}+}$ that is present over 14 $\mathrm{ps}$ of the AIMD time, followed by a reassociation and a redissociation, which suggests a dynamic change between these two species.

In the high-pressure runs at 4.5 $\mathrm{GPa}$ (# 13–22), the average Y coordination is about eight (see Fig. ). In the case of the Y–chloride complexes, the dissociation of the 1-fold and 3-fold coordinated complexes is observed in run # 13 and 15 (see Fig. ). Only in run # 14 does the initial Y–chloride complex $\left[{\mathrm{YCl}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$ persist over the whole 27 $\mathrm{ps}$ trajectory. The higher-coordinated Y–chloride complexes break apart within the equilibration run and the results are not further analyzed. This breakdown is partly driven by the association of the coordinating chloride with sodium. For instance, in run # 13, one sodium chloride unit associates with the Y complex before the chloride dissociates from the Y complex and $[{\mathrm{NaCl}}_{\mathrm{3}}{]}^{\mathrm{2}-}$ is formed for $\sim \mathrm{3}$ $\mathrm{ps}$ (Fig. ). The resulting $\left[\mathrm{Y}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{8}}{]}^{\mathrm{3}+}$ associates with ${\mathrm{OH}}^{-}$ shortly afterwards. In all high $P$ runs, the formation of ${\mathrm{OH}}^{-}$ by self-dissociation of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ close to the yttrium ion can be seen as in the low $P$ runs.

Figure 6

Formation of $\left[\mathrm{YOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$ and reassociation of initial complex $\left[\mathrm{YCl}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$ in run # 13. The blue-colored water molecules are located in the second hydration shell, red–white ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ and ${\mathrm{OH}}^{-}$ are bonded in the first coordination shell, and sodium is colored yellow. Chloride ions are in cyan and the ${\mathrm{Y}}^{\mathrm{3}+}$ ion is copper colored. The red numbers indicate the time progress in the AIMD simulation. In the center, the proton transfer state is highlighted with a green sphere.

[Figure omitted. See PDF]

Figure  illustrates the ${\mathrm{OH}}^{-}$ formation mechanism as it evolves for Y–chloride and Y–fluoride complexes at low- and high-pressure conditions for the example of $\left[\mathrm{YCl}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$ in run # 13. After the initial complex dissociates within the first 7 $\mathrm{ps}$ into $\left[\mathrm{Y}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{8}}{]}^{\mathrm{3}+}$, a proton is transferred between one ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ in the first hydration shell and a water molecule of the second shell after additional 2–3 $\mathrm{ps}$. The resulting $\left[\mathrm{YOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$ complex is present over 14 $\mathrm{ps}$, followed by reassociation with chloride. The thus-formed $\left[\mathrm{YClOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{+}$ persists during the remaining simulation time with some interruptions due to short-lived proton transfers.

For the Y–fluoride complexes, under the high-pressure conditions, the 1-fold or 2-fold Y by F coordination persists over the whole-simulation runs (# 16, 17). All complexes show association with one or two sodium ions over several picoseconds but this interaction does not lead to a dissociation of fluoride from yttrium. In # 18, the initially 3-fold coordinated complex dissociates after 14.5 $\mathrm{ps}$ and $\left[{\mathrm{YF}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{+}$ is formed. None of the mixed $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ complexes at 4.5 $\mathrm{GPa}$ persist over the entire simulation run. In each of those runs, all chloride ions are dissociated from the yttrium after at most 10 $\mathrm{ps}$, and pure $\mathrm{Y}-\mathrm{F}$ or $\mathrm{Y}-(\mathrm{F},\mathrm{OH})$ complexes remain. As in the low-pressure run for the pure hydrated ${\mathrm{Y}}^{\mathrm{3}+}$, a hydroxide ion is observed in the first hydration shell for a significant amount of simulation time.

The nearest Y-(O,Cl,F) distances show only small variations between both pressure conditions, typically in the range of 0.01–0.06 Å for the stable complexes (Table ). A closer look at the distances between oxygen of the second hydration shell and the yttrium ion (see $\mathrm{Y}-\mathrm{O}\left(\mathrm{2}\mathrm{nd}\right)$ in Table S1) reveals a continuous increase from the purely hydrated ion with increasing Cl coordination at low pressures from 4.7 to 5.2 Å. In all other cases, these distances are rather similar in a range between 4.3 and 4.6 Å.

Comparing the average halide ion coordination by ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ molecules, differences between purely hydrated halide ions or halide ions associated with the yttrium ion as well as pressure-induced changes are observed (see Table S2). For the chloride ion at 1.3 $\mathrm{GPa}$, the number of hydrating water molecules increases from two to four between ${\mathrm{YCl}}_{\mathrm{2}}^{+}-{\mathrm{YCl}}_{\mathrm{3}\mathrm{aq}}$ and dissociated ${\mathrm{Cl}}^{-}$, whereas in the mixed complexes ${\mathrm{YClF}}_{\mathrm{2}\mathrm{aq}}$ the chloride hydration number is close to three. For ${\mathrm{YCl}}^{\mathrm{2}+}$, ${\mathrm{YCl}}_{\mathrm{4}}^{-}$ and ${\mathrm{YCl}}_{\mathrm{5}}^{\mathrm{2}-}$, this number lies between two and three. The fluoride ion is coordinated by one water molecule in the Y–fluoride and mixed complexes. At 4.5 $\mathrm{GPa}$, ${\mathrm{F}}^{-}$ is hydrated by four ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ and by two when associated with the metal ion, whereas the dissociated ${\mathrm{Cl}}^{-}$ coordination increases to five solvent molecules. For those Y–chloride complexes that persist for at least 10 $\mathrm{ps}$ in the AIMD run, the chloride ion is coordinated by three water molecules.

To conclude this section, the main findings from the AIMD simulations are summarized. Firstly, the pure Y–chloride complexes ${\mathrm{YCl}}^{\mathrm{2}+}$, ${\mathrm{YCl}}_{\mathrm{2}}^{+}$ and ${\mathrm{YCl}}_{\mathrm{3}\mathrm{aq}}$ do not dissociate within the course of the simulation of at least 23 $\mathrm{ps}$ at a pressure of 1.3 $\mathrm{GPa}$ ($\mathit{\rho }=\mathrm{1072}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) but they do at higher pressure (4.5 $\mathrm{GPa}$, 1447 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) except for YCl${}_{\mathrm{2}}^{+}$. Furthermore, ${\mathrm{YCl}}_{\mathrm{4}}^{-}$ and ${\mathrm{YCl}}_{\mathrm{5}}^{\mathrm{2}-}$ decompose in the unbiased AIMD simulations. Secondly, the pure Y–fluoride complexes persist over the simulation time at 1.3 $\mathrm{GPa}$. It is important to state that at these conditions the formation of hydrofluoric acid is very strong when the fluoride is not associated with the metal ion. At 4.5 $\mathrm{GPa}$, the neutral complex $\left[{\mathrm{YF}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}_{\mathrm{aq}}$ dissociates and a lower coordinated species forms. Thirdly, ${\mathrm{OH}}^{-}$ is formed in the first hydration shell of Y–chloride and Y–fluoride complexes due to the self-dissociation of water. Its abundance increases with decreasing number halide ligands. Furthermore, chloride as well as fluoride form mixed complexes with yttrium and hydroxide at both $P/T$ conditions, whereas mixed $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ complexes are rather unstable. The overall coordination of yttrium changes from $\sim \mathrm{7}$ at 1.3 $\mathrm{GPa}$ to $\sim \mathrm{8}$ at 4.5 $\mathrm{GPa}$.

3.2 Free energy exploration

Although several complexes are observed in some of the runs described above, it is not feasible to derive the corresponding formation constants of the aqueous species directly from the AIMD simulations. This would require much longer simulation times to ensure the statistical significance of the relative species abundance. On the tens of picoseconds timescale, we did not observe a complete exchange of the halogen ligands of the Y–chloride or Y–fluoride complexes except for run # 15. To overcome this problem, we apply the TI approach within a constraint AIMD simulation to compute the reaction free energies of aqueous complex dissociation. As described in the previous section, the mixed $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ complexes have a tendency to dissociate already during the conventional AIMD runs, which indicates their low stability. Therefore, no TI runs are performed for those complexes. When the constraint halide ion associates with hydrogen during the constraint MD, the simulation is stopped and the integration step is repeated from a different starting configuration. Thus, we confirm that all results are reproducible within a series of simulations. The derived Helmholtz free energies (${\mathrm{\Delta }}_{\mathrm{r}}A$) of Reactions () and () are listed in Table .

Table 3

List of the formation Gibbs free energies ($\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$) derived from ab initio thermodynamic integration. The error of the free energies is estimated with 5 $\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$.

#	Reaction	Cell	Temperature (${}^{\circ }$C)	Pressure ($\mathrm{GPa}$)	$\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ (Å)	${\mathrm{\Delta }}_{\mathrm{r}}G$	${\mathrm{\Delta }}_{\mathrm{r}}G{}^{\circ }$	$\log K$
TI-1	$[\mathrm{YOH}{]}^{\mathrm{2}+}+{\mathrm{Cl}}^{-}\to [\mathrm{YClOH}{]}^{+}$	A1	800	1.3	2.57	$-$36.1	$-$57.4	2.79
TI-2	$[\mathrm{YClOH}{]}^{+}+{\mathrm{Cl}}^{-}+{\mathrm{H}}_{\mathrm{3}}{\mathrm{O}}^{+}\to [{\mathrm{YCl}}_{\mathrm{2}}{]}^{+}+\mathrm{2}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$	A1	800	1.3	2.57	$-$38.8	$-$57.6	2.80
TI-3	$[{\mathrm{YCl}}_{\mathrm{2}}{]}^{+}+{\mathrm{Cl}}^{-}\to [{\mathrm{YCl}}_{\mathrm{3}}{]}_{\mathrm{aq}}$	A1	800	1.3	2.64	$-$26.4	$-$41.8	2.03
TI-4	$[\mathrm{YOH}{]}^{\mathrm{2}+}+{\mathrm{Cl}}^{-}\to [\mathrm{YClOH}{]}^{+}$	B1	800	4.5	2.62	$-$29.6	$-$57.4	2.79
TI-5	$[\mathrm{YClOH}{]}^{+}+{\mathrm{Cl}}^{-}+{\mathrm{H}}_{\mathrm{3}}{\mathrm{O}}^{+}\to [{\mathrm{YCl}}_{\mathrm{2}}{]}^{+}+\mathrm{2}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$	B1	800	4.5	2.64	$-$8.5	$-$34.4	1.68
TI-6	$[\mathrm{YOH}{]}^{\mathrm{2}+}+{\mathrm{F}}^{-}\to [\mathrm{YFOH}{]}^{+}$	B5	800	4.5	2.12	$-$45.9	$-$85.6	4.17
TI-7	$[\mathrm{YF}{]}^{\mathrm{2}+}+{\mathrm{F}}^{-}\to [{\mathrm{YF}}_{\mathrm{2}}{]}^{+}$	B5	800	4.5	2.13	$-$38.5	$-$74.7	3.64
TI-8	$[{\mathrm{YF}}_{\mathrm{2}}{]}^{+}+{\mathrm{F}}^{-}\to [{\mathrm{YF}}_{\mathrm{3}}{]}_{\mathrm{aq}}$	B5	800	4.5	2.14	$-$36.3	$-$72.5	3.53

At 1.3 $\mathrm{GPa}$, in TI-1 starting from $\left[\mathrm{YCl}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{\mathrm{2}+}$, in nearly all TI steps close to the yttrium ${\mathrm{OH}}^{-}$ formed by the hydrolysis of water molecules within the first 2–4 $\mathrm{ps}$. Therefore, the obtained dissociation energy of 36.1 $\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$ does not distinguish between $\left[\mathrm{YCl}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{\mathrm{2}+}$ and $\left[\mathrm{YOHCl}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}^{+}$ complexes. For TI-2, we observe a similar dissociation energy of 38.8 $\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$ but much fewer hydroxide ions are formed such that in average over all integration steps $\mathrm{Y}-{\mathrm{OH}}^{-}$ appears in only 14 % of the total simulation time. The lowest dissociation energy of the Y–chloride complexes at 1.3 $\mathrm{GPa}$ occurs in TI-3 with 26.4 $\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$. In TI-3, only little amounts of ${\mathrm{OH}}^{-}$ are formed. While the integration proceeds, the yttrium–oxygen coordination changes (including ${\mathrm{OH}}^{-}$ and ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$) for all complexes at a constraint $\mathrm{Y}-\mathrm{Cl}$ distance between 3.6 and 4.0 Å. The removed chloride as well as the remaining Y–chloride complex associate with sodium during the dissociation process for a few picoseconds. Dissociation energies of the Y–fluoride complexes at 1.3 $\mathrm{GPa}$ could not be obtained due to the strong association of ${\mathrm{H}}^{+}$ and ${\mathrm{F}}^{-}$. The formation of hydrofluoric acid prevents the required long constraint Y–fluoride bond distances for the integration of the PMF.

Figure 7

(a) Potential of mean force of TI-4. The green circle indicates the constraint distance of 3.6 Å for which the force and cumulative mean force is shown in panel (b). Dots in panel (c) indicate the presence of $\mathrm{Y}-{\mathrm{OH}}^{-}$ bonds during the simulation run.

[Figure omitted. See PDF]

As in the lower-pressure runs, the integration at high pressure does not distinguish between complexes in which ${\mathrm{OH}}^{-}$ is present or absent. In TI-4, starting from $\left[\mathrm{YClOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{+}$ and forming $\left[\mathrm{YOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$, a free energy difference of 29.6 $\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$ is obtained. As illustrated in Fig. , the reassociation of ${\mathrm{OH}}^{-}$ with an excess proton leads to a change of the constraint force due to the decreasing attraction of the metal cation to the constraint chloride ligand. The most frequent $\mathrm{Y}-{\mathrm{OH}}^{-}$ association is observed at constraint Y–Cl distances above 3.0 Å.

The dissociation energies of the Y–chloride complexes significantly decrease between $\left[\mathrm{YCl}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$ and $\left[{\mathrm{YCl}}_{\mathrm{2}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}}{]}^{+}$ at 4.5 GPa. TI-5 yields a low dissociation energy of 8.5 $\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$ where in several of the integration steps the second chloride ligand also dissociates and $\left[\mathrm{YOH}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{7}}{]}^{\mathrm{2}+}$ is formed. The $\mathrm{Y}-\mathrm{Cl}$ dissociation is preceded by an $\mathrm{Y}-{\mathrm{OH}}^{-}$ association and results in the formation of $\left[{\mathrm{YOH}}_{\mathrm{0}-\mathrm{1}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{6}-\mathrm{7}}{]}^{\mathrm{3}-\mathrm{2}+}$. In those cases, the initial complex was reset and the integration step was restarted. For $\left[{\mathrm{YCl}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}_{\mathrm{aq}}$, the dissociation energy was not derived because the initial complex dissociated at short $\mathrm{Y}-\mathrm{Cl}$ constraint distances within the first picosecond of each simulation. Therefore, this complex is not considered as important at such high-pressure conditions. An approach to derive a dissociation energy for such unstable complexes was shown by by constraining the remaining ligands at the equilibrium distance.

Table 4

List of the logarithmic stability constants $\log {\mathit{\beta }}^{*}$ derived from ab initio thermodynamic integration in comparison to theoretical predictions based on HKF regression.

Density ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$)	$P$ ($\mathrm{GPa}$)${}^{\mathrm{b}}$	Reaction	$\log {\mathit{\beta }}^{*{}^{\mathrm{a}}}$	$\log {\mathit{\beta }}^{\mathrm{c}}$
1072	1.3	${\mathrm{Y}}^{\mathrm{3}+}+{\mathrm{Cl}}^{-}\to {\mathrm{YCl}}^{\mathrm{2}+}$	2.8	3.4${}^{\mathrm{e}}$
1072	1.3	${\mathrm{Y}}^{\mathrm{3}+}+\mathrm{2}{\mathrm{Cl}}^{-}\to {\mathrm{YCl}}_{\mathrm{2}}^{+}$	5.6	4.6${}^{\mathrm{e}}$
1072	1.3	${\mathrm{Y}}^{\mathrm{3}+}+\mathrm{3}{\mathrm{Cl}}^{-}\to {\mathrm{YCl}}_{\mathrm{3}\mathrm{aq}}$	8.4	–
1447	4.5	${\mathrm{Y}}^{\mathrm{3}+}+{\mathrm{Cl}}^{-}\to {\mathrm{YCl}}^{\mathrm{2}+}$	2.8	$-$1.7${}^{\mathrm{e}}$
1447	4.5	${\mathrm{Y}}^{\mathrm{3}+}+\mathrm{2}{\mathrm{Cl}}^{-}\to {\mathrm{YCl}}_{\mathrm{2}}^{+}$	4.5	0.04${}^{\mathrm{e}}$
1447	4.5	${\mathrm{Y}}^{\mathrm{3}+}+{\mathrm{F}}^{-}\to {\mathrm{YF}}^{\mathrm{2}+}$	4.2	3.4${}^{\mathrm{d}}$
1447	4.5	${\mathrm{Y}}^{\mathrm{3}+}+\mathrm{2}{\mathrm{F}}^{-}\to {\mathrm{YF}}_{\mathrm{2}}^{+}$	7.8	10.9${}^{\mathrm{d}}$
1447	4.5	${\mathrm{Y}}^{\mathrm{3}+}+\mathrm{3}{\mathrm{F}}^{-}\to {\mathrm{YF}}_{\mathrm{3}\mathrm{aq}}$	11.2	–

${}^{\mathrm{a}}$ Listed values involve further transition states that cannot be separated during the TI. ${}^{\mathrm{b}}$ Pressure estimates assuming a 2 $\mathrm{molal}$ $\mathrm{NaCl}$ solution using the empirical equation of state by . ${}^{\mathrm{c}}$ Values computed using the DEW model with HKF parameters reported by ${}^{\mathrm{d}}$ for Y and ${}^{\mathrm{e}}$ for Ho complexes.

In TI-6 to TI-8, the dissociation energies of Y–fluoride complexes following Reaction () are investigated. The equilibrium distance between the yttrium ion and its fluoride ligands is smaller than the Y–Cl bond distance. Due to this shorter distance, the integration range to reach the dissociated state was reduced to 5.0 Å. As shown in Fig.  (and Fig. S1), the convergence of the free energy is still reached. In each run of this simulation series, we observe the temporary formation of hydrofluoric acid by the protolysis from ${\mathrm{H}}_{\mathrm{3}}{\mathrm{O}}^{+}$ to one of the constraint ${\mathrm{F}}^{-}$and, compared to the low-pressure runs, to a lesser extent by the protolysis of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$.

The dissociation energies (Fig. ) are quite similar between the Y–fluoride complexes. During the integration, the Y complexes as well as the removing fluoride ions interact with sodium. No self-dissociation of the complexes is observed except for $\left[{\mathrm{YF}}_{\mathrm{3}}\right({\mathrm{H}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{5}}{]}_{\mathrm{aq}}$. In the latter case, at a constant distance of 2.6 Å, one of unconstrained fluorides separates from the initial complex. However, this behavior is not reproducible.

During all integration runs at both pressures, the $\mathrm{Y}-\mathrm{O}$ coordination number increases in average by one during the transformation from the associated with the dissociated state. For the Y–chloride complexes, the increase in hydration happens at a $\mathrm{Y}-\mathrm{Cl}$ distance between 3.6 and 4.1 Å, which is in the range of the repulsive part of the constraint force. For the dissociation of Y–fluoride complexes, this distance decreases to 3.4–3.6 Å. A significant influence of the second hydration shell on the constraint force of the reaction coordinate is not observed. All complexes and the removing halide ions interact with sodium. This interaction cannot be quantified by the PMF but the association of the yttrium complex with sodium decreases with the number of halide ligands.

Figure 8

Evolution of the Helmholtz free energy derived from thermodynamic integration for $\mathrm{Y}-\mathrm{Cl}/\mathrm{F}$ complexes at 800 ${}^{\circ }$C and $\mathrm{1.3}/\mathrm{4.5}$ $\mathrm{GPa}$ Inset: (a) potential of mean force of the dissociation reaction of $[\mathrm{YFOH}{]}^{+}$ to ${\mathrm{YOH}}_{\mathrm{2}}^{+}$ for a integration length of 5.0 and 6.0 Å. (b) Resulting Helmholtz free energy along the integration pathway (for higher magnification, see Fig. S1).

[Figure omitted. See PDF]

3.3 Thermodynamic data

Finally, the reaction free energies derived from thermodynamic integration are transformed into standard state properties by applying the activity corrections described in the methods section above. The standard state correction yields ${\mathrm{\Delta }}_{\mathrm{r}}G{}^{\circ }$ values that are significantly smaller compared to ${\mathrm{\Delta }}_{\mathrm{r}}G$ (see Table ). This treatment does not consider explicitly the formation of $\mathrm{HCl}$ and $\mathrm{HF}$ as well as the association of yttrium with hydroxide because their formation during TI is not systematic and it is not possible to quantify their contributions to the reaction free energies. As these contributions seem not to be negligible as shown in Fig. , the logarithmic stability constants include not only the reactions listed in Table  and are therefore indicated with a star ($\log {\mathit{\beta }}^{*}$). Figure  shows the ${\mathrm{\Delta }}_{\mathrm{r}}G{}^{\circ }$ and $\log {\mathit{\beta }}^{*}$ for the different species. While pressure does not affect the formation reaction of ${\mathrm{YCl}}^{\mathrm{2}+}$, the stability of ${\mathrm{YCl}}_{\mathrm{2}}^{+}$ decreases substantially with increasing pressure. Comparing the stability constants of chloride and fluoride species at 4.5 $\mathrm{GPa}$, it is found that those of the fluoride species are higher by 1.4 ($\log {\mathit{\beta }}_{\mathrm{1}}^{*}$) to 3.3 ($\log {\mathit{\beta }}_{\mathrm{2}}^{*}$) log units.

4 Discussion

4.1 Molecular structure of the aqueous complexes

As mentioned in the introduction, the number of studies focusing on the hydration environment of yttrium or other HREEs in aqueous fluids at high $T$ and $P$ conditions is very limited. The average coordination of seven nearest neighbors that is observed in the simulations at 1.3 $\mathrm{GPa}$ fits in the range of experimental results. observed eight to nine nearest neighbors at lower $T$ of 250 ${}^{\circ }$C and vapor pressures in highly concentrated chloride solution. But they did not find an association of Y–chloride, whereas reported a strong association with an average coordination of four at 500 ${}^{\circ }$C. The present simulations predict that ${\mathrm{YCl}}_{\mathrm{4}}^{-}$ is not stable at high $P/T$ conditions. The reason for this could be the too-low average Y coordination in ${\mathrm{YCl}}_{\mathrm{4}}^{-}$ and ${\mathrm{YCl}}_{\mathrm{5}}^{\mathrm{2}-}$. An average coordination of seven as present in the stable Y–chloride complexes of the simulations cannot be achieved in these highly chlorinated species due to steric constraints. The Y–O distances derived from EXAFS spectra by in the range of 2.36–2.39 Å are in good agreement with the atomic distances from the presented simulations, while the conference abstract by does not comprise quantitative data. Experiments and simulations are only partly comparable, as the simulations were performed at higher $T$ (800 ${}^{\circ }$C) than the experiments by or .

Figure 9

(a) Reaction Gibbs free energy ${\mathrm{\Delta }}_{\mathrm{r}}G{}^{\circ }$ of the different formation reactions; (b) change of the logarithmic stability constant of the different $\mathrm{Y}-(\mathrm{Cl},\mathrm{F})$ complexes.

[Figure omitted. See PDF]

The formation of stable $\mathrm{Y}-\mathrm{Cl}$ species at 1.3 and 4.5 $\mathrm{GPa}$ is consistent with observations by, e.g., , that the mobility of yttrium increases with increasing availability of ${\mathrm{Cl}}^{-}$ in aqueous systems. The Y–Cl complexes become less stable with increasing $P$. The destabilization of metal–halide species with rising pressure is known from a variety of systems (see overview by ). This is commonly explained by the increase of the dielectric constant with increasing density at constant $T$. The increase of $\mathit{ϵ}$ ($\mathit{ϵ}=\mathrm{17.1}$ at 1072 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ and $\mathit{ϵ}=\mathrm{26.2}$ at 1447 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ according to ) leads to a stronger hydration of the metal ion by ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ and the stabilization of charged species. Therefore, ${\mathrm{YCl}}_{\mathrm{3}\mathrm{aq}}$ decomposes at 4.5 $\mathrm{GPa}$. The direct competition of both halide ligands in the mixed complexes shows a clear preference of yttrium to form stable bonds with ${\mathrm{F}}^{-}$ rather than with ${\mathrm{Cl}}^{-}$. Moreover, the lower reaction Gibbs free energies of the Y–fluoride species in Fig. a strongly support this observation.

Figure  shows a comparison of the $\mathrm{Y}-{\mathrm{OH}}^{-}$ formation between both pressure conditions in those aqueous complexes that persist over at least 10 $\mathrm{ps}$ in the unconstrained AIMD simulations. A higher abundance of hydroxide groups is observed for $\mathrm{Y}-\mathrm{Cl}$ and $\mathrm{Y}-\mathrm{F}$ complexes at lower $P$. Furthermore, the amount of formed ${\mathrm{OH}}^{-}$ decreases with increasing number of ligands and is particularly high for the purely hydrated ${\mathrm{Y}}^{\mathrm{3}+}$ at both pressures. This observation cannot be explained by the increased self-dissociation of water, which increases with pressure (see, e.g., ). According to , the ${\mathrm{OH}}^{-}$ molality is in the range of 10${}^{-\mathrm{6}}$ to 10${}^{-\mathrm{5}}$ under the investigated $P/T$ conditions. Therefore, the hydroxide formation in the simulation must be driven by charge compensation in the absence of other ligands around the yttrium. The low abundance of $\mathrm{Y}-{\mathrm{OH}}^{-}$ in the neutral species (e.g., ${\mathrm{YF}}_{\mathrm{3}\mathrm{aq}}$) supports this interpretation. It has to be stressed that the initial simulation cells do not contain excess hydroxide ions. Therefore, the $\mathrm{Y}-{\mathrm{OH}}^{-}$ association is expected to be much higher if the ${\mathrm{OH}}^{-}$ concentration increases. But it has to be underlined that the values in Fig.  are rather imprecise because an equilibrium distribution of Y–hydroxide bonds cannot be achieved in the rather short simulation time.

Figure 10

Comparison of the average $\mathrm{Y}-{\mathrm{OH}}^{-}$ coordination number between the different complexes that exist for at least over 10 $\mathrm{ps}$ in the unconstrained AIMD at 1.3 (a) and 4.5 GPa (b).

[Figure omitted. See PDF]

The association of yttrium with hydroxide as observed in the simulations was also noticed in high $P/T$ solubility experiments by in $\mathrm{NaCl}$ brines but not in the ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}-\mathrm{NaF}$ system . The authors explained the association by a geometrical effect. The smaller HREEs (in comparison to a LREE) have less attraction to a large chloride ion due to the so-called “steric hindrance”, as discussed by . However, in the case of the Y–fluoride complexes, especially for ${\mathrm{YF}}^{\mathrm{2}+}$, the ${\mathrm{OH}}^{-}$ formation is also controlled by the protolysis of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ close to the metal ion. Therefore, the geometrical explanation does not hold to explain the Y–hydroxide bonding. The fact that this process was not detected in the experiments by might be caused by the high fluoride content in the experiments. The majority of ${\mathrm{YF}}_{\mathrm{2}}^{+}$ in the experiments underlines this conclusion. The fact that the formation of $\mathrm{Y}-(\mathrm{F}/\mathrm{OH})$ species was not detected in solubility experiments by up to 250 ${}^{\circ }$C indicates that the protolysis of vicinal water by yttrium is a high temperature process. Protolysis at high temperatures was also reported by for ${\mathrm{Ti}}^{\mathrm{4}+}$ and might be a general property of high-field-strength elements.

Entropy is an additional driving force for the ion association due to changes in the local solvent structure near the critical point and above as discussed, e.g., by and based on AIMD simulations. The dominant effect arises from the translational entropy through hydration changes of the ions during the aqueous reaction. Such a concept was already discussed by . They proposed that the change in the electrostriction volume of the solvent controls the association of aqueous species, e.g., ${\mathrm{NH}}_{\mathrm{3}\mathrm{aq}}$, ${\mathrm{HCl}}_{\mathrm{aq}}$ and ${\mathrm{NaCl}}_{\mathrm{aq}}$. In the present simulation study, a relation between the change of hydration of chloride and the stability of certain complexes is observed at 1.3 $\mathrm{GPa}$. For example, the formation of ${\mathrm{YCl}}_{\mathrm{2}}^{+}$ and ${\mathrm{YCl}}_{\mathrm{3}\mathrm{aq}}$ (TI-2, TI-3) releases two inner-sphere solvation water molecules because the hydration of ${\mathrm{Cl}}^{-}$ decreases from four to two. This may explain the very similar reaction Gibbs energy of $-$38.8 $\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$ in TI-2 compared to TI-1 ($-$36.1 $\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$), where only one ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ is released (besides the effect of ${\mathrm{OH}}^{-}$ formation). Normally, one would expect a decrease in released reaction free energy with increasing number of ligands. Due to the spontaneous formation of $\mathrm{HF}$, the number of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ in the hydration shell of an ${\mathrm{F}}^{-}$ion cannot be derived directly from the 1.3 $\mathrm{GPa}$ simulations. At 4.5 $\mathrm{GPa}$, two hydration ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ are exchanged for ${\mathrm{F}}^{-}$and ${\mathrm{Cl}}^{-}$ during the dissociation of the initial complexes. This convergence of the hydration change for the halogens supports the assumptions by and that the entropic effect decreases with increasing density at constant temperature ($\ge \mathrm{200}$ ${}^{\circ }$C).

4.2 Comparison of the thermodynamic data with HKF regression

Experimental high $P$ and high $T$ thermodynamic properties of Y–chloride and Y–fluoride species are not available to compare and evaluate the present simulation results. Therefore, the Deep Earth Water (DEW) model is used to compute stability constants of ${\mathrm{YF}}^{\mathrm{2}+}$ and ${\mathrm{YF}}_{\mathrm{2}}^{+}$ derived from solubility experiments up to 250 ${}^{\circ }$C by using a HKF regression to 800 ${}^{\circ }$C and a fluid density equal to that of the simulations. There are no stability data of Y–chloride complexes available but due to the similarities of yttrium– and holmium (Ho)–chloride complexes at room temperature , the behavior of Y–chloride complexes is assumed to be similar to Ho–chloride complexes . The Ho–chloride HKF parameters are taken from , . In addition, $\mathrm{Y}/\mathrm{Ho}-{\mathrm{OH}}^{-}$ stabilities are derived from data of , . The results are shown in Fig. .

Comparing the results in Fig. a, it can be observed that the stabilities of ${\mathrm{YCl}}^{\mathrm{2}+}$ and ${\mathrm{YCl}}_{\mathrm{2}}^{+}$ are similar to those of the Ho–chloride species within approximately one logarithmic unit. For ${\mathit{\beta }}_{\mathrm{3}}$ ($\mathrm{Y}-\mathrm{Cl}$), Y and Ho show opposite behavior. In this case, an increase in the stability of ${\mathrm{YCl}}_{\mathrm{3}\mathrm{aq}}$ is observed, whereas the Ho ${\mathit{\beta }}_{\mathrm{3}}$ ($\mathrm{Ho}-\mathrm{Cl}$) decreases . do not report any ${\mathit{\beta }}_{\mathrm{3}}$ (Cl) values. For Y–hydroxide and Ho–hydroxide species, the stabilities are in the range of ${\mathrm{YCl}}^{\mathrm{2}+}$. At higher density (Fig. b), the HKF model yields significantly lower stability constants of Ho–/Y–chloride species compared to the AIMD for ${\mathrm{YCl}}^{\mathrm{2}+}$ and ${\mathrm{YCl}}_{\mathrm{2}}^{+}$. The formation of $\mathrm{Y}-{\mathrm{OH}}^{-}$ in the AIMD runs suggests that $\mathrm{Y}-{\mathrm{OH}}^{-}$ association may occur in high-density brines.

Due to the lack of AIMD $\log {\mathit{\beta }}_n$ ($\mathrm{Y}-\mathrm{F}$) data, Fig. c only shows values derived from the HKF parameters. Here, the Y–fluoride species have the highest stability compared to Ho–fluoride , Y–hydroxide and Ho–hydroxide species. It should be mentioned here that the model from is suspected to overestimate the HREE-F stability . In Fig. d, it is shown that for the higher fluid density the ${\mathit{\beta }}_{\mathrm{1}}^{*}$ ($\mathrm{Y}-\mathrm{F}$) values from the simulations are consistent with regression based on experimental data within one log unit. The formation constant of ${\mathrm{YF}}_{\mathrm{2}}^{-}$ from is higher compared to the AIMD results. On the other hand, the Ho–fluoride complexes and the Y–/Ho–hydroxide complexes have a much lower stability. Those differences indicate a different behavior of the geochemical twins (Y and Ho) in fluoride-rich environments, which could explain the fractionation of these elements even at high $P/T$ conditions as it was observed, e.g., in hydrothermal systems and discussed by . Furthermore, the comparable stability of both Y–chloride and Ho–chloride complexes confirms the comparable geochemical behavior of the ions in chloride-rich solutions as assumed by .

In general, at lower fluid densities, we note similar stabilities of Y– and Ho–chloride complexes. The strong divergence for the neutral complexes might be explained by the origin of the HKF parameters. The data from were derived from thermodynamic predictions based on measurements at 25 ${}^{\circ }$C and 1 $\mathrm{bar}$ where only a very small amount of a neutral species is formed and the uncertainty of the extrapolation is large. This interpretation is supported by reported experiments of and , who do not observe the formation of neutral complexes up to 250 ${}^{\circ }$C. Only in situ measurements by show higher Cl coordination of yttrium, whereas for Ho no data are available.

Figure 11

Comparison of the aqueous species stability derived from the AIMD simulation and HKF regression using the DEW model at 800 ${}^{\circ }$C and a fluid density of 1072 and 1447 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$.

[Figure omitted. See PDF]

4.3

The role of chloride and fluoride for the mobilization of ${\mathrm{Y}}^{\mathrm{3}+}$ under subduction-zone conditions

Stable Y–chloride and Y–fluoride complexes are found over the investigated $P/T$ range. Fluoride-bearing species are predominant in chloride- and fluoride-rich metamorphic fluids. This outcome is consistent with a variety of studies from the field (see, e.g., ) as well as from experiments . Figure a–c show the distribution of the Y–fluoride and Y–chloride complexes as a function of the dissolved salt concentration. In Fig. d, the speciation change as function of the ${\mathrm{F}}^{-}$concentration is shown. Besides the presented thermodynamic data from this study, competing aqueous reactions as the formation of ${\mathrm{HCl}}_{\mathrm{aq}}$, ${\mathrm{NaCl}}_{\mathrm{aq}}$, ${\mathrm{NaF}}_{\mathrm{aq}}$, ${\mathrm{HF}}_{\mathrm{aq}}$ and ${\mathrm{NaOH}}_{\mathrm{aq}}$ as expected in high-grade metamorphic fluids and observed in the simulations are included. Possible complexation with silica components due to the lack of thermodynamic data as well as mineral reactions are not included. The neutral p$\mathrm{H}/$pOH is fixed by the ion product of water . The thermodynamic properties of the competing aqueous species at high $P/T$ conditions are computed using the HKF model parameters reported by , , and using the DEW model . For this simple model, a Y concentration of 23 $\mathrm{ppm}$ in solution is assumed. This amount corresponds to measurements in subducted material such as mid-ocean-ridge basalt and is in the range of Y solubility in ${\mathrm{Cl}}^{-}$-bearing brines reported by .

Figure 12

Y–chloride and Y–fluoride species distribution computed from the AIMD data assuming 23 $\mathrm{ppm}$ dissolved yttrium in aqueous solution at 800 ${}^{\circ }$C. (a, b) Y–chloride species distribution versus the logarithmic sodium chloride concentration in low- and high-density fluid; (c) Y–fluoride species versus the amount of dissolved $\mathrm{NaF}$ and (d) $\mathrm{Y}-\mathrm{Cl}/\mathrm{F}$ species in a 1 $\mathrm{molal}$ $\mathrm{NaCl}$ solution with increasing ${\mathrm{F}}^{-}$ concentration. The pH of the solutions is changing in a range of $\pm \mathrm{1.7}$ over the concentration ranges.

[Figure omitted. See PDF]

Comparing the Y–chloride speciation in Fig. a–b, it can be noted that yttrium is mainly associated with ${\mathrm{Cl}}^{-}$ at an $\mathrm{NaCl}$ concentration of $\sim \mathrm{0.005}$ $\mathrm{molal}$ in the lower-density fluid, whereas under higher-density conditions this state is reached at $\sim \mathrm{0.01}$ $\mathrm{molal}$ $\mathrm{NaCl}$ in solution. For a fluid density of 1072 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, ${\mathrm{YCl}}_{\mathrm{3}}$ is the main species above an $\mathrm{NaCl}$ concentration of 0.01 $\mathrm{molal}$. Under high-density conditions, ${\mathrm{YCl}}_{\mathrm{2}}^{+}$ is the major species at a higher $\mathrm{NaCl}$ concentration of 0.14 $\mathrm{molal}$. This shows that rather small amounts of dissolved $\mathrm{NaCl}$ are needed to form Y–chloride species under subduction-zone $P/T$ conditions if the yttrium is in solution. The required amounts of chloride presumably occur in subduction fluids as demonstrated in a variety of studies . However, it has been discussed by that fluoride aqueous complexes are much more capable to mobilize significant amounts of yttrium. As shown in Fig. c for a fluid-bearing $\mathrm{Y}+\mathrm{NaF}$ (in the absence of other ligands) at an $\mathrm{NaF}$ concentration of about 0.0014 $\mathrm{molal}$, Y–fluoride complexes are the majority species. Deep aqueous fluids in the Earth's crust are presumed to be mixtures of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ plus salt (mainly $\mathrm{NaCl}$, ${\mathrm{CaCl}}_{\mathrm{2}}$ and $\mathrm{KCl}$) . Therefore, in Fig. d, the speciation in a 1 $\mathrm{molal}$ $\mathrm{NaCl}$ brine with increasing amount of fluoride in the solution at a density of 1447 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ is computed. Here, ${\mathrm{YCl}}_{\mathrm{2}}^{+}$ and ${\mathrm{YF}}_{\mathrm{3}\mathrm{aq}}$ are the dominant species. Below a ${\mathrm{F}}^{-}$ concentration of approximately 0.01 $\mathrm{molal}$, Y–chloride complexes are formed, and above this concentration, Y–fluoride becomes more important, at least in high-density brines. At lower densities, the formation of Y–fluoride evolves at lower fluoride concentration due to the strong increase of complex stability as shown in Fig. c.

Estimates of the fluoride content in aqueous phases evolving in subducting slabs are given in the range of $\sim \mathrm{0.003}$–0.18 $\mathrm{molal}$ derived from analysis of melt inclusions and metamorphic rocks , thermodynamic modeling and the fugacity ratio $f_{\mathrm{HF}}/f_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ based on partitioning data between mineral and fluid . The simple model outlined above shows that only a low fluoride concentration in the metamorphic aqueous fluids is needed to change the yttrium speciation and therefore its mobility. This outcome is in line with recent thermodynamic modeling by that shows that a small amount of dissolved ${\mathrm{F}}^{-}$ (tens of ppm) in ore-forming solutions in equilibrium with a granite assemblage is able to mobilize significant amount of REEs. But the high immobility of yttrium in crustal as well as subduction-zone metasomatism might reflect a low fluorine activity in the majority of metasomatic fluids due to formation of fluoride species, e.g., ${\mathrm{HF}}_{\mathrm{aq}}$, ${\mathrm{HF}}_{\mathrm{2}}^{-}$, ${\mathrm{SiF}}_{\mathrm{6}}^{\mathrm{2}-}$ or $\mathrm{Si}(\mathrm{OH}{)}_{\mathrm{2}}{\mathrm{F}}_{\mathrm{2}\mathrm{aq}}$. Therefore, the included thermodynamic properties of $\mathrm{HF}$ from might underestimate its stability. A very low solubility of yttrium-bearing minerals as suggested by and the retrograde solubility of REE phosphates could play an important role in the fixation of yttrium and other HREEs.

At high densities, the stability of Y–hydroxide complexes might be higher than the HKF regression indicates. As shown by , at room temperature, the $\mathrm{Y}-{\mathrm{OH}}^{-}$ complexation is sensitive to pH changes and could evolve under more alkaline conditions as presumed for deep metasomatism . But with the applied methods in this study, no further statement can be made due to the limitations of the PMF thermodynamic integration approach to extract reaction free energies for all relevant individual reactions including different intermediate states and dynamic changes such as proton transfer. To overcome this problem, multiple constraints could be applied as suggested by but this would lead to even longer simulation times and might require more than the available computation resources for reaching sufficient convergence of structures and energies . Therefore, other free energy sampling methods could be promising alternatives, e.g., the metadynamics approach . This method has been successfully applied to compute acid constant (pK${}_{\mathrm{a}}$) in combination with quantum mechanical molecular dynamics (e.g., ) and to find new reaction pathways (see ). To probe aqueous systems or mineral–fluid interactions of more realistic size and composition, classical molecular dynamics simulations in combination with force field interaction potentials or/and machine learning potentials certainly have potential to provide significant progress in this field.

5 Conclusions

The results of the ab initio molecular dynamics simulations provide new insight into the $\mathrm{Y}-\mathrm{Cl}$, $\mathrm{Y}-\mathrm{F}$ and $\mathrm{Y}-{\mathrm{OH}}^{-}$ complexation in highly saline solutions, as they occur in geological settings, e.g., of subduction zones. Firstly, Y–chloride aqueous complexes are observed at 800 ${}^{\circ }$C and 1.3 or 4.5 $\mathrm{GPa}$ where ${\mathrm{YCl}}_{\mathrm{3}}$ and ${\mathrm{YCl}}_{\mathrm{2}}^{+}$ are the major species. Moreover, the destabilization of ${\mathrm{YCl}}_{\mathrm{4}}^{-}$ and ${\mathrm{YCl}}_{\mathrm{5}}^{\mathrm{2}-}$ indicates that there are no other Y–chloride species that have to be considered at least in high-grade metamorphic processes.

The extracted thermodynamic properties of Y–chloride species presented in this study are the first published data set to our knowledge. The stability of Ho–chloride complexes derived from thermodynamic calculations based on an HKF regression at a solution density of 1073 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ suggests that Y and Ho behave very similar in ${\mathrm{Cl}}^{-}$ rich solutions but with increasing solution density Y–chloride complexes are more favorable than Ho–chloride complexes. On the contrary, Y shows a much stronger association with fluoride compared to Ho, which could explain their different behavior in F-rich aqueous environments. A different association behavior of both elements with respect to ${\mathrm{OH}}^{-}$ would have an even higher impact on the $\mathrm{Y}/\mathrm{Ho}$ fractionation because mineral solubilities and mineral surface adsorption are much more controlled by the pH and pOH values than the halide content of the aqueous fluid.

The formation of Y–fluoride complexes in high-density aqueous fluids happens even at very low ${\mathrm{F}}^{-}$concentrations and should lead to a high mobility of Y (HREE) as observed in natural samples. Only in very fluorine-rich environments (e.g., ) significant amounts of HREEs are mobilized. This finding may indicate that the fluoride activity in the majority of metamorphic aqueous fluids is rather low. The reason for that could be the high affinity of fluorine to silicate , which is one of the main components of aqueous phases in metamorphic processes . In high-grade metamorphic fluids, it is also necessary to consider that more cations are competing for fluoride to form, e.g., $\mathrm{MgF}{}^{+}$ and $\mathrm{CaF}{}^{+}$. This leads to a decrease of the F${}^{+}$ activity as well.

As discussed by , the HREE mass change by fluid–rock interaction is much more determined by the mineral assemblages and phosphate mobility and therefore the halide content of the fluid phase might not be the only controlling factor for the HREE mobility. Nevertheless, the thermodynamic data reported here are consistent with the results of the HKF regression. Furthermore, the stability constants are affected by the formation of hydroxide mixed complexes and $\mathrm{HCl}/\mathrm{HF}$ formation during the thermodynamic integration. Therefore, the presented thermodynamic quantities can only be considered as semi-quantitative. Furthermore, it must be emphasized that the applied activity correction could be a source of huge uncertainty. As demonstrated by , the Ewald summation, that is used to build up the periodic boundary condition within the simulation, disrupts the solvation free energy of highly charged ions. But this kind of perturbation is not accounted in the Debye–Hückel approach. Therefore, a more systematic evaluation of the impact of artificial periodic electrostatics and neutralizing background charge on the computed thermodynamic properties derived from AIMD simulations especially at high temperatures is required in future studies.

Code and data availability

All simulations were performed using the CPMD code that is distributed free of charge to non-profit organizations under the CPMD Free Licence (https://www.cpmd.org, last access: 27 April 2020). The molecular structures were visualized using the VMD code (; https://www.ks.uiuc.edu/Research/vmd/, last access: 27 April 2020), which is distributed under an Open Source Licence.

The supplement related to this article is available online at: https://doi.org/10.5194/se-11-767-2020-supplement.

Author contributions

JS and SJ designed the study. JS performed the numerical simulations and SJ supervised the analysis and interpretation of the results. Both authors discussed the results of the study and wrote the manuscript.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Exploring new frontiers in fluids processes in subduction zones”. It is a result of the EGU Galileo conference “Exploring new frontiers in fluids processes in subduction zones”, Leibnitz, Austria, 24–29 June 2018.

Acknowledgements

This study was supported by the Deutsche Forschungsgemeinschaft in the framework of project JA 1469/10-1. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (https://www.gauss-centre.eu, last access: 27 April 2020) for funding this project (CHPO15) by providing computing time through the John von Neumann Institute for Computing (NIC) on the JUWELS supercomputer and JURECA at Jülich Supercomputing Centre (JSC). We appreciate helpful discussions with Thomas Wagner. Yuan Mei and an anonymous reviewer are thanked for their very useful comments and suggestions to improve the manuscript.

Financial support

This research has been supported by the DFG (grant no. JA 1469/10-1) and the Gauss Center for Supercomputing (grant no. CHPO15).

Review statement

This paper was edited by Simone Tumiati and reviewed by Yuan Mei and one anonymous referee.