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

The rare earth elements (REE) are important geochemical tracers for geological processes such as high grade metamorphism. Aqueous fluids are considered important carriers for the REE in a variety of geological environments including settings associated with subduction zones. The capacity of a fluid to mobilize REE 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 fluoride species at a temperature of 800 ◦C in a pressure range between 1.3 5 and 4.5 GPa derived from ab initio molecular dynamics simulations. The total yttrium coordination by H2O 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 at conditions relevant to slab dehydration. Mixed Y-(Cl,F) complexes are found to be unstable. Furthermore, in contrast to field observations thermodynamic modeling indicates that yttrium should be mobilized at 10 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 REE. 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. 15


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 (Tang et al., 2016). Magnetotelluric anomalies (e.g. Worzewski et al. (2011);McGary et al. (2014)) suggest the occurrence of a high proportion of melts and water-rich fluids in the subducted slabs due to partial melting 20 (Zheng et al., 2016) and due to the dehydration of water-bearing minerals such as serpentine (Ulmer and Trommsdorff, 1995) and amphibole (Schmidt and Poli, 1998). Quite naturally these fluids do not consist of pure water much rather they are brines

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) (Hohenberg and Kohn, 1964;Kohn and Sham, 1965). Here, we used AIMD simulations with the  Parrinello (Car and Parrinello, 1985) method to model the molecular structure of Y-(Cl,F) complexes in aqueous solutions. We performed simulations with the widely used CPMD code (CPMD, 1990;Marx and Hutter, 2000). Within the code the BLYP exchange correlation functional (Becke, 1988) was employed and the plane-wave expansion of the Kohn-Sham orbitals was truncated at a cutoff energy of 80 Ry. To reduce the computation afford the core electrons of all atoms in the cubic simulation cell were approximated by Goedecker-type pseudopotentials (Goedecker et al., 1996;Hartwigsen et al., 1998;Krack, 2005). 95 To separate the electronic and nuclear motion of the Car-Parrinello molecular dynamics, a fictive electron mass of 600 a.u. with fictitious kinetic energy of 0.24 a.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 • C (so called N V T ensemble) and with a time step of 0.1 fs. The temperature in the simulation was controlled using a Nosé thermostat (Nosé, 1984;Hoover, 1985). We chose two different pressure conditions of 1.3 GPa and 4.5 GPa for this study. The volumes of the simulation cell were estimated from 100 the correlation function provided by Mantegazzi et al. (2013) assuming 2 molal NaCl solution for all cells (see Tab. 1). Due to this approximation the pressures listed in Tab. 1 have to be considered estimates.
The initial atomic configurations were derived from AIMD simulations of pure NaCl solutions (200 H 2 O and 10 NaCl).
This configuration had been equilibrated for a few tens of picoseconds (ps) at 1000 K and a simulation box size of 13.74 Å.
The original NaCl solutions were generated in classical molecular dynamics simulations using MCY potentials (Matsuoka 105 et al., 1976). 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 Tab. 1. Figure 1 shows a snapshot of the simulation cell A1 with a [YCl 3 (H 2 O) 4 ] 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 Tab. 1) were generated from this initial one and equilibrated for several picoseconds. 110 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 F − ions were included to avoid the formation of hydrofluoric acid by the reaction H 2 O + F − OH − + HF.

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 115 distribution functions g ij (r). These functions describe the probability for finding a pair of atoms of elements j and i at a 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 kg m −3 (1.3 GPa) and 1447 kg m −3 (4.5 GPa).  distance r normalized to the particle number density ρ N of the fluid: where c i (= N i /N ) and c j are the concentrations of elements i and j, N is the number of particles in the simulation box, δ(x) the Dirac delta function, and R a and R b are the position vectors of particles a and b. In the numerical implementation of Eq. 1 120 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 (r) 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 125 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 (r). The association of OH − groups with cations is evaluated by considering oxygen atoms coordinated by one hydrogen only. To distinguish pure OH − from two H 2 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 distinguishing criterion. Only if the oxygen of the next water molecule 130 is located in a distance of > 1.6 Å the OH − is accounted as hydroxide. This value represents approximately the hydrogen bond distance between OH − and H 2 O (Stefanski et al., 2018). 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 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.

135
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 (Resat and Mezei, 1993). Here, we use the constrained molecular dynamics approach and thermodynamic integration in terms of the blue moon sampling (Ciccotti et al., 2005). 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 140 temperature conditions (Bühl and Golubnychiy, 2007;Bühl and Grenthe, 2011) but also at hydrothermal (Mei et al., 2013(Mei et al., , 2015(Mei et al., , 2016 and deep crustal high density fluid conditions (Mei et al., 2018).
Using this method, the Helmholtz free energy difference (∆ r A) of a chemical dissociation reaction is obtained from the average forces F (r) between two atoms under the constraint of keeping their bond distance r constant. ∆ r A is derived from the integration of F (r) between two different distances r 1 and r 2 , which correspond to the associated and the dissociated state 145 of the aqueous complexes (Sprik and Ciccotti, 1998): The formal relation between the Helmholtz free energy and the Gibbs free energy is given by: V is the volume of the simulation cell. Here we use an N V T ensemble and the change in pressure averaged over the whole 150 trajectory is approximately zero ( 2 1 dP = 0). To look into the formation of Y-(Cl,F) species in equilibrium reactions we removed one of the ligands from the Y ion in multiple integration steps by constraining the Y-(Cl/F) distance. During the integration the Y-Cl/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 2 illustrates an example of the dissociation reaction in the simulation box at 1.3 GPa and 800 • C: The integration starts at the first distance ( Fig. 2 (a)), which corresponds approximately to the equilibrium distance of the Y-Cl contact ion pair. This interatomic distance (a) is estimated as the intercept with the zero force line from a linear interpolation between the first and the second integration step, 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. 2 (b)) until a water molecule takes the place 160 of the ion in the first coordination shell around the Y 3+ (Fig. 2 (c)). At this point the force becomes repulsive. By integration over the potential of mean force (PMF) (Fig. 2 (I)) the Helmholtz free energy (∆ r A) difference between the initial complex and the product of the reaction is derived in (Fig. 2 (II)). 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 Supporting Information Fig. SI1). In Fig. 2  average force is computed. As convergence criterion a value of 5 kJ mol −1 over the last 2 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 between 4.5 and 40 ps.

Thermodynamic approaches
It is textbook knowledge that the formation of monomeric complexes in equilibrium reactions develops in steps (Atkins and 170 De Paula, J., 2006;Brown and Ekberg, 2016). During the formation process a ligand L is added to the metal cation M. This formation of a ML n complex can be written as a sequence of stepwise reactions:

175
ML n−1 + L → ML n with the respective logarithmic equilibrium constants Having determined equilibrium constants for all of those reactions, the stability constant (also referred to as cumulative stability constant or overall stability constant) β n of species ML n is defined as 180 log β n = log K 1 + log K 2 + log K 3 · · · log K n (5) The standard Gibbs free energy (∆ r G • i ) depends on the reaction Gibbs free energy (∆ r G i ) derived from the MD simulation, temperature T , gas constant R, molality of the ions m i and the activity coefficient γ i : The standard state for a solute in aqueous solution is a 1 molal hypothetical solution with properties of a infinitely diluted 185 solution (IUPAC, 1982). The concentration and behavior of the solutes in the simulation cell (see Tab. 1) is quite different from this standard state. Therefore, we computed the activity coefficient corresponding to this hypothetical solution using the B-dot model (Helgeson et al., 1981;Helgeson, 1969), which is an empirical extension of the Debye-Hückel theory (Hückel and Debye, 1923): where z i is the charge of ion i,å the mean distance of closest approach between the ions and I the ionic strength: They depend on temperature, density (ρ H2O ) and dielectric constant ( ) of the solvent. is computed using the equation provided by Pan et al. (2013) and Sverjensky et al. (2014) for pure H 2 O assuming a fluid density of the simulation. The value ofḂ is calculated by the CHNOSZ software package (Dick, 2008) applying the extrapolation suggested by Manning et al. (2013). Theå parameter for the different complexes and ions (for all Y-Cl/F complexes a constant value of 4.5 Å is applied) are taken from Kielland (1937) and Eq. 7 was solved in a Python implementation of the EQBRM program (Anderson and Crerar, 200 1993).
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: 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. Different partial radial distribution functions for Y-(Cl − , F − , OH − , OH 2 ) are shown in Fig. 3. To facilitate the comparison of the first peaks the g ij (r) 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 Y 3+ and F − (∼2.1 Å) followed by Y-OH − (∼2.1-2.2 Å) and Y-OH 2 (∼2.3-2.4 Å). The largest distance is observed for Y-Cl pairs (∼2.6-2.8 Å) (for an overview 220 see Tab. 2 and Tab. SI1). Further, in g Y-O (r) a second maximum is observed. It corresponds to the 2 nd hydration shell, which 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 abundance complexes 1 of the last 10 ps of the simulation are listed. is formed around all complexes at all studied P/T conditions. For all fluid compositions, association of NaCl n (n=1-3) and of Na + with the yttrium complex is observed. The mean Na-Cl coordination, NaCl species distribution and second hydration shell positions are listed in Tab. SI1. #5, even if the initial coordination is lower (runs #1-3). In runs #1-3, the initial Y-Cl coordination is retained over the whole AIMD run time, whereas for #4 and #5 the fourfold and fivefold Cl-coordination 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 fourfold coordinated chloride complex is dissociated a total coordination number of seven is reached at the end of the simulation run. Frequently, the formation of OH − in the first hydration shell of yttrium is observed. The major 230 hydroxide formation mechanism will be discussed below. 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 5.5 6.0 6.5 7.0 7.5 #1 to #5 the average Y coordination decreases with the increasing number of initial chloride ligands (Fig. 4) distances of the mixed complexes are comparable to those of the pure ones. In run #12 starting from [Y(H 2 O) 7 ] 3+ hydroxide is formed within the first 8 ps (see Fig. 5), which results in the formation of that is present over 14 ps of the AIMD time followed by a reassociation and a redissociation, which suggests a dynamic change between these two species.

255
In the high pressure runs at 4.5 GPa (#13-22) the average Y coordination is about eight (see Fig. 4). In case of the Y chloride complexes, the dissociation of the onefold and threefold coordinated complexes is observed in runs #13 and #15 (see Fig. 5). Only in run #14 the initial Y chloride complex [YCl 2 (H 2 O) 5 ] + persists over the whole 27 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 260 chloride unit associates with the Y complex before the chloride dissociates from the Y complex and [NaCl 3 ] 2− is formed for ∼3 ps (Fig. 6). The resulting [Y(H 2 O) 8 ] 3+ associates with OH − shortly afterwards. In all high P runs, the formation of OH − by self-dissociation of H 2 O close to the yttrium ion can be seen as in the low P runs. Comparing the average halide ion coordination by H 2 O molecules, differences between purely hydrated halide ions or halide ions associated to the yttrium ion as well as pressure-induced changes are observed (see Supporting Information Tab. SI2). For the chloride ion at 1.3 GPa the number of hydrating water molecules increases from two to four between YCl + 2 -YCl 3 aq and dissociated Cl − , whereas in the mixed complexes YClF 2 aq the chloride hydration number is close to three. For YCl 2+ , YCl − 4 285 and YCl 2− 5 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 GPa, F − is hydrated by four H 2 O and by two when associated with the metal ion, whereas the dissociated Cl − coordination increases to five solvent molecules. For those Y chloride complexes that persist for at least 10 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 com-290 plexes YCl 2+ , YCl + 2 and YCl 3 aq do not dissociate within the course of the simulation of at least 23 ps at a pressure of 1.3 GPa (ρ =1072 kg m −3 ) but they do at higher pressure (4.  forms. Thirdly, 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 fluoride form mixed complexes with yttrium and hydroxide at both P/T conditions whereas mixed Y-(Cl,F) complexes are rather unstable. The overall coordination of yttrium changes from ∼7 at 1.3 GPa to ∼8 at 4.5 GPa.

Free energy exploration 300
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 thermodynamic integration (TI) approach within a constraint AIMD simulation to compute the 305 reaction free energies of aqueous complex dissociation. As described in the previous section the mixed Y-(Cl,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.  In those cases, the initial complex was reset and the integration step was 330  Table 3.
List of the formation Gibbs free energies (kJ mol −1 ) derived from ab initio thermodynamic integration. The error of the free energies is estimated with  In TI-6 to TI-8 the dissociation energies of Y fluoride complexes following reaction R3 are investigated. The equilibrium 335 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. 8 (and Fig. SI1) 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 H 3 O + to one of the constraint F − and, compared to the low pressure runs, to a lesser extent by the protolysis of H 2 O.

340
The dissociation energies (Fig. 8) Figure 8. Evolution of the Helmholtz free energy derived from thermodynamic integration for Y-Cl/F complexes at 800 • C and 1.3/4.5 GPa Inset: potential of mean force of the dissociation reaction of [YFOH] + to YOH + 2 for a integration length of 5.0 Å and 6.0 Å. (II) resulting Helmholtz free energy along the integration pathway (for higher magnification see Fig. SI1). cannot be quantified by the PMF but the association of the yttrium complex with sodium decreases with the number of halide 350 ligands.

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 ∆ r G • that are significantly smaller compared to ∆ r G (see Tab. 3). This treatment does not consider explicitly the formation of HCl and HF 355 as well 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. 7 the logarithmic stability constants include not only the reactions listed in Tab. 4 and are therefore indicated with a star (log β * ). Figure 9 shows the ∆ r G • and log β * for the different species. While pressure does not affect the formation reaction of YCl 2+ , the stability of YCl + 2 decreases substantially with increasing pressure. Comparing the stability constants of chloride 360 and fluoride species at 4.5 GPa those of the fluoride species are higher by 1.4 (log β * 1 ) to 3.3 (log β * 2 ) log units.   (Mantegazzi et al., 2013), † values computed using the DEW model (Sverjensky et al., 2014)  in aqueous fluids at high T and P conditions is very limited. The average coordination of seven nearest neighbors that is nearest neighbors at lower T of 250 • C and vapour pressures in highly concentrated chloride solution. But they did not find an association of Y chloride, whereas Mayanovic et al. (2002) reported a strong association with an average coordination of four at 500 • C. The present simulations predict that YCl − 4 is not stable at high P/T conditions. The reason for this could be the too low average Y coordination in YCl − 4 and YCl 2− 5 . An average coordination of seven as present in the stable Y chloride 370 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 Vala Ragnarsdottir et al. (1998)  shows a clear preference of yttrium to form stable bonds with F − rather than with Cl − . Moreover, the lower reaction Gibbs free energies of the Y fluoride species in Fig. 9 (a) strongly support this observation. Figure 10 shows a comparison of the Y-OH − formation between both pressure conditions in those aqueous complexes that 385 persist over at least 10 ps in the unconstrained AIMD simulations. A higher abundance of hydroxide groups is observed for Y-(Cl/F) complexes at lower P . Furthermore, the amount of formed OH − decreases with increasing number of ligands and is particularly high for the purely hydrated Y 3+ at both pressures. This observation cannot be explained by the increased selfdissociation of water, which increases with pressure (see e.g. Rozsa et al. (2018); Goncharov et al. (2005)). According to Marshall and Franck (1981) the OH − molality is in the range of 10 -6 to 10 -5 under the investigated P/T conditions. Therefore, 390 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 Y-OH − in the neutral species (e.g. YF 3 aq ) supports this interpretation. It has to be stressed that the initial simulation cells do not contain excess hydroxide ions. Therefore, the Y-OH − association is expected to be much higher if the OH − concentration increases. But it has to be underlined that the values in Fig. 10 are rather imprecise because an equilibrium distribution of Y hydroxide bonds cannot be achieved in the rather short simulation time.

395
The association of yttrium with hydroxide as observed in the simulations was also noticed in high P/T solubility experiments by Tropper et al. (2011) in NaCl brines but not in the H 2 O−NaF system (Tropper et al., 2013). The authors explained the association by a geometrical effect. The smaller HREE (in comparison to a LREE) have less attraction to a large chloride ion due to the so called "steric hindrance" as discussed by Mayanovic et al. (2009). But in case of the Y fluoride complexes especially for YF 2+ the OH − formation is also controlled by the protolysis of H 2 O close to the metal ion. Therefore, the geo-400 metrical explanation does not hold to explain the Y hydroxide bonding. That this process was not detected in the experiments  by Tropper et al. (2013) might be caused by the high fluoride content in the experiments. The majority of YF + 2 in the experiments underline this conclusion. That the formation of Y-(F/OH) species was not detected in solubility experiments by Loges et al. (2013) up to 250 • C indicates that the protolysis of vicinal water by yttrium is a high temperature process. Protolysis at high temperatures was also reported by van Sijl et al. (2010) for Ti 4+ and might be a general property of high-field-strength 405 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 Sherman (2010) and Mei et al. (2014) 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 Mesmer et al. (1988). They proposed that the change in the electrostriction volume of the solvent controls 410 the association of aqueous species, e.g. NH 3 aq , HCl aq and NaCl 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 GPa. For example the formation of YCl + 2 and YCl 3 aq (TI-2, TI-3) releases two inner-sphere solvation water molecules because the hydration of Cl − decreases from four to two. This may explain the very similar reaction Gibbs energy of -38.8 kJ mol −1 in TI-2 compared to TI-1 (-36.1 kJ mol −1 ) where only one H 2 O is released (besides the effect of OH − formation). Normally, one would expect a decrease 415 in released reaction free energy with increasing number of ligands. Due to the spontaneous formation of HF the number of H 2 O in the hydration shell of an F − ion cannot be derived directly from the 1.3 GPa simulations. At 4.5 GPa two hydration H 2 O are exchanged for F − and Cl − during the dissociation of the initial complexes. This convergence of the hydration change for the halogens supports the assumptions by Mesmer et al. (1988) and Mei et al. (2015) that the entropic effect decreases with increasing density at constant temperature (≥ 200 • C).

Comparison of the thermodynamic data with HKF regression
Experimental high P and high T thermodynamic properties of Y chloride and fluoride species are not available to compare and evaluate the present simulation results. Therefore, the Deep Earth Water (DEW) model (Sverjensky et al., 2014) is used to compute stability constants of YF 2+ and YF + 2 derived from solubility experiments up to 250 • C by Loges et al. (2013) using a HKF regression to 800 • C and a fluid density equal to that of the simulations. There are no stability data of Y chloride 425 complexes available but due to the similarities of yttrium and holmium (Ho) chloride complexes at room temperature (Luo and Byrne, 2001) the behavior of Y chloride complexes is assumed to be similar to Ho chloride complexes (Migdisov et al., 2019).
Comparing the results in Fig. 11 (a) it can be observed that the stabilities of YCl 2+ and YCl + 2 are similar to those of the Ho 430 chloride species within approximately one logarithmic unit. For β 3 (Y-Cl), Y and Ho show opposite behavior. In this case an increase in the stability of YCl 3 aq is observed whereas the Ho β 3 (Ho-Cl) decreases (Haas et al., 1995). Migdisov et al. (2009) do not report any β 3 (Cl) values. For Y and Ho hydroxide species the stabilities are in the range of YCl 2+ . At higher density ( Fig. 11 (b)) the HKF model yields significantly lower stability constants of Ho/Y chloride species compared to the AIMD for YCl 2+ and YCl + 2 . The formation of Y-OH − in the AIMD runs suggests that Y-OH − association may occur in high density 435 brines.
Due to the lack of AIMD log β n (Y−F) data, Fig. 11 (c) only shows values derived from the HKF parameters. Here, the Y fluoride species (Loges et al., 2013) have the highest stability compared to Ho fluoride (Migdisov et al., 2009;Haas et al., 1995), Y hydroxide (Shock et al., 1997) and Ho hydroxide (Haas et al., 1995) species. It should be mentioned here that the model from (Haas et al., 1995) is suspected to overestimate the HREE-F stability (Migdisov and Williams-Jones, 2014). In  (Migdisov et al., 2009) complexes and the Y/Hohydroxide 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 445 it was observed e.g. in hydrothermal systems (Bau and Dulski, 1995) and discussed by Loges et al. (2013). 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 Migdisov et al. (2019).
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 Haas et al. (1995) were derived 450 from thermodynamic predictions based on measurements at 25 • C and 1 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 Migdisov et al. (2009) andLoges et al. (2013) who do not observe the formation of neutral complexes up to 250 • C. Only in situ measurements by (Mayanovic et al., 2002) show higher Cl coordination of yttrium whereas for Ho no data is available. YCl 3 n n (this study) HoCl 3 n n (Haas et. al 1995) YOH 2 + (Shock et al., 1997) HoOH 2 + (Haas et al. 1995) HoCl 3 n n ( Migdisov et al., 2009)  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. Newton et al. (1998)) as well as from experiments (Hetherington et al., 2010;Tsay et al., 2014;Tropper et al., 2013). Fig-ure12(a-c) shows the distribution of the Y fluoride and chloride complexes as a function of the dissolved salt concentration. In Fig. 12 (d) the speciation change as function of the F − concentration is shown. Beside the presented thermodynamic data from 460 this study competing aqueous reactions as the formation of HCl aq , NaCl aq , NaF aq , HF aq , NaOH aq as expected in high grade metamorphic fluids (Aranovich and Safonov, 2018;Galvez et al., 2016;Manning, 2018)) 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 pH/pOH is fixed by the ion product of water (Marshall and Franck, 1981). The thermodynamic properties of the competing aqueous species at high P/T conditions are computed using the HKF model parameters reported 465 by Mei et al. (2018); Shock and Helgeson (1988); Shock et al. (1997) and Shock et al. (1989) using the DEW model (Sver-jensky et al., 2014). For this simple model, a Y concentration of 23 ppm in solution is assumed. This amount corresponds to measurements in subducted material such as mid-ocean-ridge basalt (Sun and McDonough, 1989) and is in the range of Y solubility in Cl − -bearing brines reported by Tropper et al. (2011);Schmidt et al. (2007a). The pH of the solutions is changing in a range of ±1.7 over the concentration ranges.
Comparing the Y chloride speciation in Fig. 12 (a,b) it can be noted that yttrium is mainly associated with Cl − at an NaCl 470 concentration of ∼0.005 molal in the lower density fluid whereas at the higher density conditions this state is reached at ∼0.01 molal NaCl in solution. For a fluid density of 1072 kg m −3 YCl 3 is the main species above an NaCl concentration of 0.01 molal. At high density conditions YCl + 2 is the major species at a higher NaCl concentration of 0.14 molal. This shows that rather small amounts of dissolved NaCl are needed to form Y chloride species at subduction zone P/T conditions if the yttrium is in solution. The required amounts of chloride presumably occurs in subduction fluids as demonstrated in a variety 475 of studies (Barnes et al., 2018;Aranovich and Safonov, 2018). However, it has been discussed by Tropper et al. (2013) that fluoride aqueous complexes are much more capable to mobilize significant amounts of yttrium. As shown in Fig. 12 (c) for a fluid bearing Y + NaF (in the absence of other ligands) at an NaF concentration of about 0.0014 molal, Y fluoride complexes are the majority species. Deep aqueous fluids in the Earth crust are presumed to be mixtures of H 2 O+salt (mainly NaCl, fluoride in the solution at a density of 1447 kg m −3 is computed. Here, YCl + 2 and YF 3aq are the dominant species. Below a F − concentration of approximately 0.01 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 11 (c).
Estimates of the fluoride content in aqueous phases evolving in subducting slabs are given in the range of ∼0.003-0.18 molal 485 derived from analysis of melt inclusions and metamorphic rocks (Portnyagin et al., 2007;Hughes et al., 2018;Aranovich and Safonov, 2018), thermodynamic modeling (Zhu and Sverjensky, 1991) and f HF/f H 2 O ratio based on partitioning data between mineral and fluid (Sallet, 2000;Yardley, 1985). 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 Xing et al. (2018) that shows that a small amount of dissolved 490 F − (tens of ppm) in ore forming solutions in equilibrium with a granite assemblage is able to mobilize significant amount of REE. But the high immobility of yttrium in crustal as well subduction zone metasomatism (Ague, 2017) might reflect a low fluorine activity in the majority of metasomatic fluids due to formation of fluoride species, e.g. HF aq , HF − 2 , SiF 2− 6 or Si(OH) 2 F 2 aq . Therefore the included thermodynamic properties of HF from Shock et al. (1989) might underestimate its stability. A very low solubility of yttrium-bearing minerals as suggested by Migdisov et al. (2016) and the retrograde solubility 495 of REE phosphates (Schmidt et al., 2007a;Cetiner et al., 2005) could play an important role in the fixation of yttrium and other HREE.
At high densities, the stability of Y hydroxide complexes might be higher than the HKF regression indicates. As shown by Liu et al. (2012), at room temperature the Y-OH − complexation is sensitive to pH changes and could evolve under more alkaline conditions as presumed for deep metasomatism (Galvez et al., 2016). But with the applied methods in this study no further 500 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 Ivanov et al. (2006) 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 (Mark et al., 1994). Therefore, other free energies sampling methods could be promising alterna-505 tives, e.g. the metadynamics approach (Laio and Parrinello, 2002). This method has been successfully applied to compute acid constant (pKs) in combination with quantum mechanical molecular dynamics (e.g. Sakti et al. (2018); Tummanapelli and Vasudevan (2015)) and to find new reaction pathways (see Pérez de Alba Ortíz et al. (2018); Pietrucci et al. (2018)). 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 (Behler and Parrinello, 2007;Cheng 510 et al., 2019) certainly have potential to provide significant progress in this field.
The results of the ab initio molecular dynamics simulations provide new insight into the Y-Cl, Y-F and Y-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 • C and 1.3 or 4.5 GPa where YCl 3 and YCl + 2 are the major species. Moreover, the destabilization of YCl 4 515 and YCl 5 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 (Migdisov et al., 2009) at a solution density of 1073 kg m −3 suggests that Y and Ho behave very similar in Cl − rich solutions 520 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 OH − would have an even higher impact on the Y/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.

525
The formation of Y fluoride complexes in high density aqueous fluids happens even at very low 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 Pan and Fleet (1996); Harlov et al. (2006)) significant amounts of HREE 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 (Dolejš and Zajacz, 2018;Dalou et al., 2015) which is one of the main components of aqueous phases in metamorphic 530 processes (Manning, 2018;Hermann et al., 2006). In high grade metamorphic fluids, it is also necessary to consider that more cations are competing for fluoride to form, e.g., MgF + and CaF + . This leads to a decrease of the F + activity as well.
As discussed by Ague (2017), 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 535 regression. Furthermore, the stability constants are affected by the formation of hydroxide mixed complexes and HCl/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 Hünenberger and McCammon (1999) does the Ewald summation, that is used to build up the periodic boundary condition within the simulation, disrupt the solvation free energy of highly charged ions. But this kind of perturbation 540 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. Loges, A., Migdisov, A. A., Wagner, T., Williams-Jones, A. E., and Markl, G.: An experimental study of the aqueous solubility and speciation of Y(III) fluoride at temperatures up to 250 • C, Geochim. Cosmochim. Acta, 123, 403-415, https://doi.org/10.1016Acta, 123, 403-415, https://doi.org/10. /j.gca.2013Acta, 123, 403-415, https://doi.org/10. .07.031, 2013.